27 марта 2011 в 19:37

Откуда берутся NaNы?

C++*
Пользователь yruslan опубликовал хорошую статью про арифметику с плавающей запятой: «Что нужно знать про арифметику с плавающей запятой».

Хочу добавить к ней пару поучительных примеров. Ситуации, которые я опишу, встречались несколько раз в моей практике. Ошибки, которые при этом порождались были очень редкими, трудно воспроизводимыми и сложными в поиске.

Возможно, я сейчас буду рассказывать прописные истины, но вызывает удивление частота, с которой люди забывают, что числа которые описывает стандарт 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, в которые передаётся вычисленный параметр, подумайте, а может ли параметр попасть на границы области определения? Если да, то будьте уверены, что когда-нибудь он выйдет за границы.
Алексей Борисов @Imp5
карма
135,0
рейтинг 0,1
Самое читаемое Разработка

Комментарии (22)

  • –5
    Многие программисты об этом забывают по тому, что арифметика с плавающей точкой встречается на порядок реже, чем целочисленная. Да и ошибки здесь на порядок веселее. Вспомнился такой код:
    double d = 0;
    for(int i=0;i<1000;i++) d += 0.2;
    bool eq = (d == 200);
    • 0
      Упс, ошибся, там d+=0,1 и d==100
      • +2
        если не видно разницы, зачем платить больше?
    • +2
      Дело в том, что такие ошибки допускают люди, которые по несколько лет почти ежедневно имеют дело с плавающей точкой. Вектора, матрицы, кватернионы каждый день, а потом раз и такая ошибка в коде. А им ведь деньги за этот код платят.
      • +4
        Большинство программистов даже зная как устроен дабл будут допускать подобные ошибки. И проблема совсем в другом — в понимании. А дабл нужно понимать как value + eps, где value то что вы хотели представить как double, а eps — относительная погрешность округления (машинное Эпсилон). Если взять пример автора и провести анализ относительных погрешностей:
        r (x1) = r (y1) = r (x2) = r (y2) = eps,
        r (x1 * x1 + y1 * y1) = r (x2 * x2 + y2 * y2) = 4 * eps,
        r (len1) = r (len2) = 3 * eps,
        r (x1 / len1) = r (y1 / len1) = r (x2 / len2) = r (y2 / len2) = 5 * eps,
        r (dotProduct) = 12 * eps, если брать x1 * x2 и y1 * y2 одного знака как в примере автора.
        Вот теперь понятно откуда берутся значения больше 1.0. Вы думаете каждый программист без соответствующего мат. образования может выполнить такой анализ?
        • +2
          Надо считать точнее: r(x1*x1+y1*y1)=2*x1*r(x1)+2*y1*r(y1)+(x1^2+y1^2+(x1^2+y1^2))*eps = 2*(x1+y1+x1^2+y1^2)*eps. В итоге результат может оказаться совершенно другим. Например при извлечении корня могут появляться члены вроде eps^{1/2}. А это уже вполне чувствительная величина.
          • 0
            Минуточку x1 * x1 и y1 * y1 строго больше нуля. Поетому r(x1 * x1 + y1 * y1) = min(r(x1*x1), r(y1*y1)) + eps = min(r(x1) + r(x1) + eps, r(y1) + r(y1) + eps) + eps = min(3*eps, 3*eps) + eps = 4*eps.
            • 0
              Поправка вместо min должен быть max, но суть таже.
            • 0
              Дело в том, что погрешность в x1 увеличивается сама по себе даже при точном возведении в квадрат — (x1+eps)*(x1+eps)=x1^2+2*x*eps+eps*eps. К этому мы добавляем погрешность произведения (она относительна в виду природы строения действительных чисел) — x1*x1*eps. Членами eps^2 пренебрегаем.
              PS На самом деле r(x1)=0, и погрешность получается меньше — 2*(x1^2+y1^2)*eps.
        • +1
          Простите, но используйте тег без него сложно читать...
  • 0
    В таком случае я не понимаю чем мой ответ отличается от ваших рассуждений. Относительная ошибка
    r(x1 * x1) у вас тоже 3 * eps.
  • +1
    Как вовремя статейка подвернулась! Я только что мучился с векторным 2d движком и как раз эта проблема была. Я топорно решил: умножал каждое число на 100 000, потом округлял, потом делил на 100 000. Ну тупо рубил количество знаков после запятой т.к. максимальная точность не нужна для этого игродвижка, как выяснилось.
  • +2
    Классика жанра, читать всем кто использует вычисления с плавающей точкой:
    What Every Computer Scientist Should Know About Floating-Point Arithmetic
    ЗЫ pdf гуглится
  • 0
    В первом примере надо сравнивать, делитель меньше FLT_MIN или нет? Прост смущает как что-то может быть меньше чем самое маленькое вещественное число.
    • +1
      FLT_MIN — это не самое маленькое вещественное число, это самое маленькое положительное денормализованное число, которое может хранить float.
      • 0
        На самом деле, самое маленькое число представимое в float. Которое обычно будет subnormal, но может и не быть, если платформа не поддерживает subnormal числа.
  • +2
    Да, таким вещам надо учить с младенчества!
    image
  • –7
    Это же язык для извращенцев, в нормальном языке программирования не должно быть запутанностей и нелогичностей, и штука вроде if (b != 0) { result = a /b } должна работать для абсолютно любых значений переменных, потому что программистам больше делать наверно нечего, как читать многостраничные мануалы по IEEE сколько-то там.

    Если бы Марк Цукенберг занимался такой ерундой, фейсбук бы не появился.
    • +2
      double и в Африке double. И если в вашем коде есть нечто вроде if (b! = 0) {result = a / b} это свидетельствует не о проблеме языка программирования, а о проблеме реализации.
      Для b == 0, a / b == INF, что соответствует реальности, поэтому не вижу никаких неточностей со стороны языка программирования. Что вы называете тогда нормальным языком программирования?
    • +2
      Толсто, очень толсто.
    • НЛО прилетело и опубликовало эту надпись здесь
  • +1
    Ждем статью «Откуда берутся infы?» :)

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.