March 27, 2011 at 19:37Where do NaNs come from?
User yruslan published a good article about floating point arithmetic: "What you need to know about floating point arithmetic . "
I want to add a couple of instructive examples to it. The situations that I will describe have met several times in my practice. Errors that were generated at the same time were very rare, difficult to reproduce, and difficult to find.
I may be going to tell uppercase truths now, but it is surprising that the frequency with which people forget that the numbers described by the IEEE754 standard is not the same as real numbers.
Before getting the return number, it would be nice to check the divisor by zero. Some programmers do this roughly like this:
Here is a small example:
And here is the result of his work: It was possible to avoid the appearance of infinity by comparing the absolute value of the divisor with the constant FLT_MIN.
Take for example the arccosine function acos (), which takes values in the interval [-1; 1]. Now we will use it to get the angle between two vectors. To do this, you need to pass the scalar product of normalized vectors to acos () and then at the output we get the angle in radians. Everything is simple. So, we write the code:
Now let's check how this code behaves with almost parallel vectors:
In this synthetic example, for 100 iterations, I got 9 angle values that contained NaN. In a real project, this error can occur extremely rarely.
In this case, it is necessary to replace acos (dotProduct) with a more secure call to acos (clamp (dotProduct, -1.f, 1.f)).
If you see division in the code, or the following functions: asin, acos, sqrt, fmod, log, log10, tan, into which the calculated parameter is passed, think about whether the parameter can fall within the boundaries of the definition domain? If so, be sure that one day he will go beyond the borders.
I want to add a couple of instructive examples to it. The situations that I will describe have met several times in my practice. Errors that were generated at the same time were very rare, difficult to reproduce, and difficult to find.
I may be going to tell uppercase truths now, but it is surprising that the frequency with which people forget that the numbers described by the IEEE754 standard is not the same as real numbers.
Division
Before getting the return number, it would be nice to check the divisor by zero. Some programmers do this roughly like this:
if (d != 0.f) result = 1.f / d;
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);
}
And here is the result of his work: It was possible to avoid the appearance of infinity by comparing the absolute value of the divisor with the constant FLT_MIN.
...
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
...
Going beyond the function definition domain
Take for example the arccosine function acos (), which takes values in the interval [-1; 1]. Now we will use it to get the angle between two vectors. To do this, you need to pass the scalar product of normalized vectors to acos () and then at the output we get the angle in radians. Everything is 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); // вывод угла между векторами
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, for 100 iterations, I got 9 angle values that contained NaN. In a real project, this error can occur extremely rarely.
In this case, it is necessary to replace acos (dotProduct) with a more secure call to acos (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, into which the calculated parameter is passed, think about whether the parameter can fall within the boundaries of the definition domain? If so, be sure that one day he will go beyond the borders.