Геометрия: абсолютная точность в предикате левого поворота
Введение
Арифметические операции над вещественными числами на компьютере не точны. Но давайте подробнее с этим разберемся.
Ведь в критически важных применениях полагаться на какой-то взятый с потолка \(epsilon=10^{-6}\) звучит очень опасно и можно либо увести беспилотный автомобиль в кювет, либо не туда и не так приземлить самолет, либо проложить маршрут роботу насквозь через стену, либо ошибиться с тем где проложена граница государства.
Причем такая проблема может всплыть не на стадии тестирования программы, а через продолжительное время эксплуатации.
Представление вещественных чисел
Вещественное 16-битное число (т.н. half
) представлено знаком, порядком и мантисcой:
-
\(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
Если мы все еще не уверены - что поделать, прибегаем к последней мере - представляем наши числа ввиде рациональных чисел и оперируем ими. Т.к. это потребует долгих вычислений и длинной арифметики - работает медленно. Но т.к. до такого наш код будет очень редко доходить, то в среднем замедленее не заметно!