1.9.1: improbable but not impossible divide-by-zero

[whydoubt]

Inside random_unit_vector(), if random_in_unit_sphere() returns a vector very close to
zero, passing that to unit_vector() could result in a divide-by-zero fault.

One way to avoid this would be if the rejection method in random_in_unit_sphere() also
rejected near-zero values. There is a near_zero() function introduced later (in 10.3), that
would be useful here.

[hollasch]

(I love diving into floating point.) Yes, the problem arises when all three components have values less than $2^{-537}$ (about $10^{-162}$), as squaring such a number yields a value less than $2^{-1074}$ (about $10^{-324}$), the smallest denormalized value, and thus underflowing to zero.

As you've mentioned, we could reject all points in the disk/sphere where the norm is less than some small value. In this case, a good lower limit would be $2^{-535} \approx 10^{-162}$. However, when selecting random points from a unit sphere or disk, it's often the case that $[0,0,0]$ points/vectors should be fine.

I think it might make sense to create a dedicated function for random_unit_vector() that doesn't make use of random_in_unit_sphere() for this reason. In addition, it would get rid of the annoying need to compute the vector magnitude twice (once to test for acceptance, and a second time for normalization). Something like this:

inline vec3 random_unit_vector() {
    while (true) {
        auto p = vec3::random(-1,1);
        auto lensq = p.length_squared();
        if (1e-160 < lensq && lensq < 1)
            return p / sqrt(lensq);
    }
}

[MaliusArth]

• std::uniform_real_distribution returns a value on the interval [a, b) --> b is excluded
  in our case 1 is being excluded which messes with the uniformity of the distribution and moves the value space a bit to the negative.
• lensq is already the squared and possibly underflowed value, so all there is left to do is to ensure it is > 0 which includes any underflowing vectors, cutting off values smaller than 1e-160 is unnecessary and just messes with the distribution further.

So the correct implementation would be

inline vec3 random_unit_vector() {
    while (true) {
        auto p = vec3::random(-1,1+std::numeric_limits<double>::epsilon); // or add epsilon inside vec3::random
        auto lensq = p.length_squared();
        if (0 < lensq && lensq <= 1)
            return p / sqrt(lensq);
    }
}

I believe :)

maybe this could be fixed as part of #1637

[hollasch]

I think you're missing the explanation above. When $\text{lensq} < 2^{-537}$, $\sqrt{\text{lensq}} = 0$. The smallest denorm will underflow the mantissa. The subsequent division by zero will yield the vector $<\infty,\infty,\infty>$. Your suggestion above still contains the same floating-point vulnerability.

Regarding your comment about excluding 1.0 from the random range, this doesn't matter. It changes the range, yes, but not the uniformity. Indeed, we can select uniformly-distributed points from a square of any width, as the resulting accepted points are then converted to normalized vectors. The disk radius could be $[0, 0.90]$, $[0, 0.5)$, or $[0, 100)$. All that matters is that the disk be sampled uniformly, and the result normalized.

[MaliusArth]

When $lensq<2^{−537}$, $\sqrt{lensq} = 0$

That statement is incorrect.
The root of a number between (0, 1) results in a larger number, not a smaller one. The smallest subnormal value of a double is 5e-324 and $\sqrt{5e-324} = 2.2227587494850775e-162$ (a bigger number), there is no danger of underflowing at this stage of the function.

In your original statement you correctly identified that the squaring is the potentially problematic part:

Yes, the problem arises when all three components have values less than
$2^{−537}$ (about $10^{−162}$), as squaring such a number yields a value less than $2^{−1074}$ (about $10^{−324}$), the smallest denormalized value, and thus underflowing to zero.

That said, while squaring has the potential of underflowing to zero, we are talking about numbers that are already virtually zero and the underflow on single vector components is not a problem.

A problem which does arise however is when all 3 components end up being 0 resulting in a squared length of 0 and a division by 0 in the normalization of the vector.
To catch this case ensure that lensq > 0