Where do NaNs come from?

User yruslan has published a good article about floating-point arithmetic: “What you need to know about floating-point arithmetic” .

I want to add to it a couple of instructive examples. The situations that I will describe have met several times in my practice. The errors that were caused by this were very rare, difficult to reproduce and difficult to find.

Perhaps I’ll be telling common truths now, but the frequency with which people forget that the numbers that the IEEE754 standard describes is not the same as real numbers is surprising.

Division

Before getting the reverse number, it would be nice to check the divisor by zero. Some programmers do it like this:
')

if (d != 0.f) result = 1.f / d; 

Completely forgetting the denormalized numbers.
Here is a small example:

  float var = 1.f; for (int i = 0; i < 50; i++) { var /= 10.f; float result = 1.f / var; printf("1 / %e = %e\n", var, result); } 

But the result of his work:
...
1 / 9.999999e-035 = 1.000000e+034
1 / 9.999999e-036 = 1.000000e+035
1 / 9.999999e-037 = 1.000000e+036
1 / 1.000000e-037 = 1.000000e+037
1 / 9.999999e-039 = 1.000000e+038
1 / 1.000000e-039 = 1.#INF00e+000 << —
1 / 9.999946e-041 = 1.#INF00e+000
1 / 9.999666e-042 = 1.#INF00e+000
1 / 1.000527e-042 = 1.#INF00e+000
1 / 9.949219e-044 = 1.#INF00e+000
1 / 9.809089e-045 = 1.#INF00e+000
1 / 1.401298e-045 = 1.#INF00e+000
...

Infinity could be avoided by comparing the absolute value of the divider with the constant FLT_MIN.

Going beyond the scope of the function

Take for example the function of arccosine acos (), which takes values ​​in the interval [-1; one]. Now we will use it to get the angle between two vectors. To do this, we need to transfer the scalar product of normalized vectors to acos () and then at the output we get an angle in radians. It's simple. So, we write the code:

  float x1 = 10.f; //    float y1 = 1.f; float length1 = sqrt(x1 * x1 + y1 * y1); x1 /= length1; y1 /= length1; //   float x2 = 2.f; //    float y2 = 10.f; float length2 = sqrt(x2 * x2 + y2 * y2); x2 /= length2; y2 /= length2; //   float dotProduct = x1 * x2 + y1 * y2; //    float angle = acos(dotProduct); printf("acos(dotProduct) = %e\n", angle); //     

And now let's check how this code behaves with almost parallel vectors:

  for (int i = 0; i < 100; i++) { float x1 = 10.f; float y1 = 5.e-3f * (rand() % 10000) / 10000; //      float length1 = sqrt(x1 * x1 + y1 * y1); x1 /= length1; y1 /= length1; float x2 = 10.f; float y2 = 5.e-3f * (rand() % 10000) / 10000; float length2 = sqrt(x2 * x2 + y2 * y2); x2 /= length2; y2 /= length2; float dotProduct = x1 * x2 + y1 * y2; float angle = acos(dotProduct); printf("dotProduct = %1.8f acos(dotProduct) = %e\n", dotProduct, angle); } 

...
dotProduct = 1.00000000 acos(dotProduct) = 0.000000e+000
dotProduct = 1.00000000 acos(dotProduct) = 0.000000e+000
dotProduct = 0.99999994 acos(dotProduct) = 3.452670e-004
dotProduct = 1.00000012 acos(dotProduct) = -1.#IND00e+000 << NaN
dotProduct = 1.00000000 acos(dotProduct) = 0.000000e+000
dotProduct = 1.00000012 acos(dotProduct) = -1.#IND00e+000 << NaN
...

In this synthetic example, at 100 iterations, I got 9 angle values ​​that contained NaN. In a real project, this error can occur very rarely.
In this case, it is necessary to replace acos (dotProduct) with a safer acos call (clamp (dotProduct, -1.f, 1.f)).

Conclusion

If you see division in the code, or the following functions: asin, acos, sqrt, fmod, log, log10, tan, to which the calculated parameter is transferred, consider whether the parameter can fall on the boundaries of the definition domain? If so, be sure that someday he will go beyond the boundaries.