Пользователь yruslan опубликовал хорошую статью про арифметику с плавающей запятой: «Что нужно знать про арифметику с плавающей запятой».
Хочу добавить к ней пару поучительных примеров. Ситуации, которые я опишу, встречались несколько раз в моей практике. Ошибки, которые при этом порождались были очень редкими, трудно воспроизводимыми и сложными в поиске.
Возможно, я сейчас буду рассказывать прописные истины, но вызывает удивление частота, с которой люди забывают, что числа которые описывает стандарт IEEE754 это не то же самое, что вещественные числа.
Перед получением обратного числа неплохо бы проверить делитель на ноль. Некоторые программисты делают это приблизительно так:
Вот небольшой пример:
А вот результат его работы:
Избежать появления бесконечности можно было сравнением абсолютного значения делителя с константой FLT_MIN.
Возьмём к примеру функцию арккосинуса acos(), которая принимает значения в промежутке [-1; 1]. Сейчас мы будем при помощи неё получать угол между двумя векторами. Для этого нужно скалярное произведение нормализованных векторов передать в acos() и тогда на выходе получим угол в радианах. Всё просто. Итак, пишем код:
А теперь проверим, как этот код поведёт себя при почти параллельных векторах:
В этом синтетическом примере на 100 итераций у меня получилось 9 значений угла, которые содержали NaN. В реальном проекте эта ошибка может проявляться крайне редко.
В данном случае необходимо заменить acos(dotProduct) на более безопасный вызов acos(clamp(dotProduct, -1.f, 1.f)).
Если вы видите в коде деление, или следующие функции: asin, acos, sqrt, fmod, log, log10, tan, в которые передаётся вычисленный параметр, подумайте, а может ли параметр попасть на границы области определения? Если да, то будьте уверены, что когда-нибудь он выйдет за границы.
Хочу добавить к ней пару поучительных примеров. Ситуации, которые я опишу, встречались несколько раз в моей практике. Ошибки, которые при этом порождались были очень редкими, трудно воспроизводимыми и сложными в поиске.
Возможно, я сейчас буду рассказывать прописные истины, но вызывает удивление частота, с которой люди забывают, что числа которые описывает стандарт IEEE754 это не то же самое, что вещественные числа.
Деление
Перед получением обратного числа неплохо бы проверить делитель на ноль. Некоторые программисты делают это приблизительно так:
if (d != 0.f) result = 1.f / d;
Вот небольшой пример:
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);
}
А вот результат его работы:
...
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
...
Избежать появления бесконечности можно было сравнением абсолютного значения делителя с константой FLT_MIN.
Выход за область определения функции
Возьмём к примеру функцию арккосинуса acos(), которая принимает значения в промежутке [-1; 1]. Сейчас мы будем при помощи неё получать угол между двумя векторами. Для этого нужно скалярное произведение нормализованных векторов передать в acos() и тогда на выходе получим угол в радианах. Всё просто. Итак, пишем код:
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); // вывод угла между векторами
А теперь проверим, как этот код поведёт себя при почти параллельных векторах:
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
...
В этом синтетическом примере на 100 итераций у меня получилось 9 значений угла, которые содержали NaN. В реальном проекте эта ошибка может проявляться крайне редко.
В данном случае необходимо заменить acos(dotProduct) на более безопасный вызов acos(clamp(dotProduct, -1.f, 1.f)).
Заключение
Если вы видите в коде деление, или следующие функции: asin, acos, sqrt, fmod, log, log10, tan, в которые передаётся вычисленный параметр, подумайте, а может ли параметр попасть на границы области определения? Если да, то будьте уверены, что когда-нибудь он выйдет за границы.
