Введение

Арифметические операции над вещественными числами на компьютере не точны. Но давайте подробнее с этим разберемся.

Ведь в критически важных применениях полагаться на какой-то взятый с потолка \(epsilon=10^{-6}\) звучит очень опасно и можно либо увести беспилотный автомобиль в кювет, либо не туда и не так приземлить самолет, либо проложить маршрут роботу насквозь через стену, либо ошибиться с тем где проложена граница государства.

Причем такая проблема может всплыть не на стадии тестирования программы, а через продолжительное время эксплуатации.

Представление вещественных чисел

Вещественное 16-битное число (т.н. half) представлено знаком, порядком и мантисcой:

Представление 16-битного вещественного числа

  • \(S\) - знак. Если бит знака равен нулю, то число положительное. Если бит знака равен единице, то число отрицательное (иначе говоря перед ним есть знак минуса).

  • \(E\) - порядок по некоторому фиксированному основанию \(B\) (так же называется экспонентой, обычно основание \(B=2\)).

  • \(M\) - мантисса.

Число \(X\) раскладывается на знак, порядок и мантиссу таким образом что:

\[X = (-1)^{S} \cdot M \cdot B^E\]

Например если основание фиксировано и традиционно равно десяти \(B=10\) то вот несколько примеров:

1) \(-234500 = (-1) \cdot 2345 \cdot 10^2\)

2) \(67.89 = (+1) \cdot 6789 \cdot 10^{-2}\)

Т.е. если число \(X = -234500\), то мантисса хранит первые несколько значащих цифр (столько сколько влезает в биты выделенные под мантиссу), а показатель (экспонента) кодирует информацию о том где находится запятая, т.е. иначе говоря о том какой порядок у данного числа.

Подробнее - представление вещественных чисел.

Стандартные типы в программирование

Самые популярные типы для хранения вещественных чисел в таком виде:

  • half (half-precision) - 16 бит (1 бит знак, 5 бит порядок, 10 бит мантисса)

  • float/single (single-precision) - 32 бита (1 бит знак, 8 бит порядок, 23 бита мантисса)

  • double (double-precision) - 64 бита (1 бит знак, 11 бит порядок, 52 бита мантисса)

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

Неточность представления

Заметим что с десятичным основанием \(B=10\) в таком виде невозможно представить абсолютно точно число \(1.3333(3)\), зато можно представить \(1.1\).

Если же основание троичное \(B=3\), то в таком виде легко представить абсолютно точно число \(1.3333(3)\), но теперь нельзя точно представить \(1.1\)!

У компьютеров в культуре и традиции другое основание - бинарное, поэтому в программировании основание равно не десяти а двум: \(B=2\). И таким образом число раскладывается ввиде:

\[X = (-1)^{S} \cdot M \cdot 2^E\]

Из этого вытекают забавные спецэффекты, например если вы попробуете сложить числа \(1.1\) и

Попробуйте посмотреть как представляются такие числа как \(1\), \(2\), \(1024\), \(3\), \(0.5\), \(0.25\), \(0.125\), \(0.375\).

Запустите данный пример:

double a = 0.7;
double b = 0.9;
double x = a + 0.1;
double y = b - 0.1;

System.out.println("x: " + x);
System.out.println("y: " + y);
System.out.println("(x==y): " + (x == y));

Он выведет:

x: 0.7999999999999999
y: 0.8
(x==y): false

Упражнение Измените константы в данном примере так чтобы несмотря на оперирование вещественными числами все результат был эквивалентен математическому.

Подсказка Вспомните что компьютер опреирует степенями двойки:

double a = 0.75;
double b = 0.25;
double x = a - 0.125;
double y = b + 0.375;

System.out.println("x: " + x);
System.out.println("y: " + y);
System.out.println("(x==y): " + (x == y));

И результат:

x: 0.625
y: 0.625
(x==y): true

Ошибки при арифметических операциях

Представьте что у нас десятичное основание и на мантиссу 3 цифры (т.е. можем хранить от \(0\) до \(999\)), и в показателе можем хранить степень от \(-10\) до \(+10\).

Пусть есть следующие числа:

  • \(2.3 = 23 \cdot 10^{-1}\).

  • \(0.00004 = 4 \cdot 10^{-5}\).

  • \(900000000 = 90 \cdot 10^8\).

Чему будут равны их попарные произведения/деления? А суммы/разницы?

Оказывается что ошибка произведения и деления не так коварна как у суммы и разницы!

Ведь прибавляя к очень большому число \(a\) очень маленькое \(b\) мы рискуем получить \(a+b=a\)! В случае когда после приведения мантисс к общему показателю значащие цифры суммы мантисс не влезают в то количество бит которое у нас выделено для представления мантиссы!

Давайте попробуем это проверить:

double powOfTwo = 1.0;
int n = 10;
for (int i = 1; i <= n; ++i) {
    powOfTwo *= 2.0;
}
System.out.println("2^" + n + "=" + powOfTwo);
double sum = powOfTwo + 1.0;
System.out.println("2^" + n + "+1.0=" + sum);
System.out.println("2^" + n + "+1.0 > 2^" + n + "? " + (sum > powOfTwo));

Данный код выводит:

2^10=1024.0
2^10+1.0=1025.0
2^10+1.0 > 2^10? true

Упражнение 1 Попробуйте увеличивать \(n\) до тех пор пока вместо true не будет получено обсурдное false.

Упражнение 2 Подумайте почему вплоть до \(n=52\) выводится true, а начиная с \(n=53\) уже false? (подсказка - вспомните сколько бит выделяется мантиссе в double)

Упражнение 3 Убедитесь что \(2^{53}+1 == 2^{53}\) с точки зрения double.

Трюк 1

Что если мы оценим погрешность вычисления и будем говорить не больше/меньше/равно нулю, а больше/меньше/не уверены как соотносится с нулем?

Но оценка погрешности ведется тоже через вещественные операции! Поэтому эту ошибку надо заложить в оценку ошибки!

Подробнее см. тут и лекцию про устойчивую реализацию алгоритмов вычислительной геометрии.

Трюк 2

Когда предикат не уверен - давайте посчитаем оценку вверх и оценку вниз (сначала сказав процессору округлять вверх, а затем посчитав тот же предикат с округлением вниз).

Если диапазон содержит ноль - мы все еще не уверен.

Трюк 3

Если мы все еще не уверены - что поделать, прибегаем к последней мере - представляем наши числа ввиде рациональных чисел и оперируем ими. Т.к. это потребует долгих вычислений и длинной арифметики - работает медленно. Но т.к. до такого наш код будет очень редко доходить, то в среднем замедленее не заметно!

Рекомендуемые источники