[Из песочницы] Вычисления с плавающей точкой: можно ли доверять результатам?
Когда речь идет о повседневных арифметических операциях, проблемы с конечной точностью вычислений не выглядят столь пугающими. И наилучшей проверкой того, что результат получен правильно, является сравнение значений полученных на различных точностях.
Если, например, вычисления, полученные на одинарной и удвоенной точностях совпадают, то создается чувство уверенности в результате, по крайней мере с точностью сопоставимой с одинарной. Здесь, я бы хотел привести один интересный пример, демонстрирующий, что даже в сравнительно несложной арифметической задаче подобная устойчивость при переменной точности представления чисел не может служить основанием для такой уверенности.
Необходимо вычислить значение функции двух переменных при определенных значениях (даны ниже) ее аргументов:
Этот пример был мною замечен, когда я разбирался с библиотекой C-XSC (система классов языка C для (преимущественно) интервальных научных расчетов с произвольной точностью). Данная библиотека отлично подходит для практического исследования вычислительной устойчивости различных численных алгоритмов.
Для эмуляции вычислений с плавающей точкой на Python установим пакет mpmath. В принципе, если ограничиться точностями IEEE 754, можно было бы использовать и NumPy, но нам нужно показать, что результаты, получаемые в рамках IEEE 754 неверны.
Вычислим значение функции f(a, b) при a = 77617 и b = 33096.
# coding: utf-8
from mpmath import *
def f(a,b):
'''
From: Ramp S.M. Algorithms for verified inclusions -- theory and practice, USA, NY, 1988.
'''
return (333.75 - a * a) * b ** 6 + a * a * (11 * a * a * b * b - 121 * b ** 4 - 2) + 5.5 * b ** 8 + a/(2.0*b)
# Одинарная точность
mp.dps = 8
a = mpf(77617)
b = mpf(33096)
print 'Значение f(a, b) с одинарной точностью:', f(a, b)
# Удвоенная точность
mp.dps = 16
a = mpf(77617)
b = mpf(33096)
print 'Значение f(a, b) с удвоенной точностью:', f(a, b)
Читатель может заметить, что задание точности через dps в mpmath это не совсем правильный подход при эмуляции реальной удвоенной точности. Если речь идет об удвоенной и\или одинарной точности в рамках IEEE то, пожалуй, более корректно описывать их характеристики через двоичную систему. Здесь, однако, это не так важно; зато использование mp.dps более проще интерпретируется — т.е. как количество десятичных значащих цифр в представлении числа.
Выполняя код, получим:
Значение f(a, b) с одинарной точностью: 1.1726039
Значение f(a, b) с удвоенной точностью: 1.172603940053179
Вполне убедительно, не так ли? Только значение-то неправильное! Правильное значение при данных a и b вообще меньше единицы!
Дополним пример следующими вычислениями:
for i in range(8, 40):
mp.dps = i
a = mpf(77617)
b = mpf(33096)
print 'Точность dps={0}, f(a,b)={1}'.format(mp.dps, f(a,b))
Получим (некоторые строки опущены):
Точность dps=8, f(a,b)=1.1726039
Точность dps=9, f(a,b)=1.17260394
Точность dps=10, f(a,b)=-7.737125246e+25
...
Точность dps=13, f(a,b)=1.172603940053
Точность dps=14, f(a,b)=1.1726039400532
Точность dps=15, f(a,b)=1.17260394005318
Точность dps=16, f(a,b)=1.172603940053179
Точность dps=17, f(a,b)=-9.2233720368547758e+18
...
Точность dps=28, f(a,b)=1.172603940053178631858834905
Точность dps=29, f(a,b)=1.1726039400531786318588349045
Точность dps=30, f(a,b)=1.17260394005317863185883490452
...
Точность dps=36, f(a,b)=-0.827396059946821368141165095479816292
Точность dps=37, f(a,b)=-0.827396059946821368141165095479816292
Точность dps=38, f(a,b)=-0.827396059946821368141165095479816292
Точность dps=39, f(a,b)=-0.827396059946821368141165095479816291999
Видно, что несмотря на устойчивость, некоторые значения все-таки существенно отличаются от обычных, что уже заставляет задуматься.
Скажу сразу, значения, получаемые при точности dps=36 и выше, правильные. Но откуда знать, что при дальнейшем увеличении точности не произойдет еще какого-либо скачка, ведь, как было видно со значением 1.17260…, даже устойчивость результата при различных точностях не может гарантировать его правильность.
К счастью, и на этот вопрос можно ответить, используя пакет mpmath. Данный пакет позволяет производить интервальные вычисления, что дает возможность оценить область возможных значений функции при представлении ее аргументов с различной точностью.
Выполним расчеты, используя аппарат интервальной арифметики mpmath:
for j in range(6, 40):
iv.dps = j
a = iv.mpf(77617)
b = iv.mpf(33096)
print 'Точность={0}: f(a, b)={1}'.format(mp.dps, f(a,b))
Получим следующее:
Точность dps=6: f(a, b)=[-5.0706024009e+30, 6.3382542101e+30]
Точность dps=7: f(a, b)=[-3.16912650057e+29, 3.16912654779e+29]
Точность dps=8: f(a, b)=[-5.94211218857e+28, 2.971056097974e+28]
Точность dps=9: f(a, b)=[-3.7138201178561e+27, 4.9517601582944e+27]
Точность dps=10: f(a, b)=[-1.54742504910673e+26, 3.09485009825849e+26]
Точность dps=11: f(a, b)=[-1.934281311383407e+25, 3.86856262277385e+25]
Точность dps=12: f(a, b)=[-4.8357032784585167e+24, 3.6267774588444373e+24]
Точность dps=13: f(a, b)=[-3.02231454903657294e+23, 2.26673591177745118e+23]
Точность dps=14: f(a, b)=[-2.833419889721787128e+22, 2.833419889721790484e+22]
Точность dps=15: f(a, b)=[-3.5417748621522339103e+21, 3.5417748621522344346e+21]
Точность dps=16: f(a, b)=[-442721857769029238784.0, 442721857769029246976.0]
Точность dps=17: f(a, b)=[-27670116110564327424.0, 27670116110564327456.0]
Точность dps=18: f(a, b)=[-3458764513820540927.0, 2305843009213693953.5]
Точность dps=19: f(a, b)=[-432345564227567614.828125, 288230376151711745.179688]
Точность dps=20: f(a, b)=[-36028797018963966.8274231, 18014398509481985.17260742]
Точность dps=21: f(a, b)=[-3377699720527870.8273963928, 1125899906842625.172604084]
Точность dps=22: f(a, b)=[-140737488355326.827396061271, 422212465065985.172603942454]
Точность dps=23: f(a, b)=[-17592186044414.8273960599472, 17592186044417.17260394006735]
Точность dps=24: f(a, b)=[-1099511627774.8273960599468637, 3298534883329.17260394005325]
Точность dps=25: f(a, b)=[-137438953470.827396059946822859, 412316860417.172603940053178917]
Точность dps=26: f(a, b)=[-17179869182.82739605994682137446, 17179869185.17260394005317863941]
Точность dps=27: f(a, b)=[-2147483646.8273960599468213681761, 2147483649.1726039400531786320407]
Точность dps=28: f(a, b)=[-268435454.827396059946821368142245, 268435457.172603940053178631864532]
Точность dps=29: f(a, b)=[-8388606.827396059946821368141165871, 8388609.172603940053178631858840746]
Точность dps=30: f(a, b)=[-1048574.8273960599468213681411651475, 1048577.1726039400531786318588349559]
Точность dps=31: f(a, b)=[-131070.827396059946821368141165095792, 131073.172603940053178631858834907439]
Точность dps=32: f(a, b)=[-8190.827396059946821368141165095483, 1.172603940053178631858834904520184761]
Точность dps=33: f(a, b)=[-1022.8273960599468213681411650954798445, 1.1726039400531786318588349045201837979]
Точность dps=34: f(a, b)=[-126.82739605994682136814116509547981678, 1.17260394005317863185883490452018372567]
Точность dps=35: f(a, b)=[-6.827396059946821368141165095479816292382, 1.172603940053178631858834904520183709123]
Точность dps=36: f(a, b)=[-0.82739605994682136814116509547981629200549, -0.82739605994682136814116509547981629181741]
Точность dps=37: f(a, b)=[-0.827396059946821368141165095479816292005489, -0.827396059946821368141165095479816291981979]
Точность dps=38: f(a, b)=[-0.8273960599468213681411650954798162919996113, -0.827396059946821368141165095479816291998142]
Точность dps=39: f(a, b)=[-0.82739605994682136814116509547981629199906033, -0.82739605994682136814116509547981629199887666]
Интересно проследить эволюцию интервала возможных значений функции: он стабилизируется только при использовании точностей от 36 и более значимых десятичных цифр, хотя и постепенно сужается.
Из интервальных расчетов становится вполне ясным, что доверять следует только результатам, получаемым при точностях от dps=36 и более десятичных цифр.
Этот пример — наглядная демонстрация того, насколько нужно быть осторожным осуществляя вычисления с плавающей точкой, и что даже 128-битная (например, numpy.float128 если говорить о Python и NumPy) точность может быть недостаточной. Он также показывает, что нельзя доверять и устойчивости результата, полученного на различных точностях. Применение аппарата интервальных вычислений может быть одним из вариантов решения в этом случае, которое позволяет оценить необходимую точность для получения адекватного результата.
Комментарии (8)
napa3um
8 июля 2016 в 12:13 (комментарий был изменён)
+1↑
↓
Сайт, посвящённый проблеме точности вычислений в разных языках: http://0.30000000000000004.com/
gbg
8 июля 2016 в 12:30
+20↑
↓
Делай раз — смотрим на формулу: у вас 3×10^4 возводится в степень 8.
Делай два — результат этого возведения — 3×10^32 — упс, вышли за приделы точности double
Делай три — открываем тетрадочку с лекциями по численным методам и читаем «Для обеспечения устойчивого счета, требуется заложить минимум 4 порядка точности арифметики». Суммируем 32+4=36. Ой, что-то теория работает.Это я к чему. Вычислительная математика работает строго по принципу GIGO — мусор на входе — мусор на выходе. Вы заранее засунули компьютеру ерунду, которую он обсчитать не может, и удивляетесь, что он выдал вам ерунду. Так математический софт не пишут.
8 июля 2016 в 16:20
+1↑
↓
Плюсанусь, однако замечу, что не столь важен максимальный порядок, сколь разница между максимальным и минимальным. Собственно, вангую, что не будь «а/2 бе», который таки порядка 10^0, то и устойчивость достигалась бы на меньшей точности, чем 36.-
gbg
8 июля 2016 в 16:29
0↑
↓
«Устойчивость» вообще термин здесь малоуместный.«Устойчивость» означает, что малые изменения параметров системы (обычно — коэффициентов диффура, или граничных условий), ведут к малым изменениям решения.
И этот термин относится к алгоритму счета, там и критерий специальный есть, и теория, которая этот критерий строго доказывает, и так далее.
Здесь же у нас чисто арифметическая задача «сколько брать знаков на расчет» и ответ к ней — давным давно эмпирически предложен — запас должен быть 4 (это как самый минимум!) знака.
Вы совершенно правы, указывая на то, что нужно брать разницу между самой большой и самой меньшей степенями десятки.
Однако, при численном решении реальных задач (в матфизике, например), задачу сначала нормируют (чтобы самое большое число по модулю было равно 1), а затем уже с нормированной и безразмерной задачей работают. Тогда критерий становится простым — «хотим N знаков точности — берем 4+N знака арифметики (а лучше — больше)».
dbagaev
8 июля 2016 в 16:34
0↑
↓
Вычислительная математика работает и по другим принципам, иногде получая вполне корректные данные на входе, но выдавая мусор на выходе, и дело часто отнюдь не только в том, что где-то мы вылазим за пределы допустимых значений или теряем разряды. Проблема в том, что некоторые вычисления являются вычислительно нестабильными, например элементарное вычитание двух чисел либо суммирование ряда.Я не буду сейчас приводить и объяснять конкретные примеры, но в одном из них выполнение серии безобидных операций без выходов за пределы чисел на выходе дает число, в котором все разряды кроме старшего неправильные. Есть и контр-примеры, в которых изменение порядка вычислений сохраняет вычислительную стабильность даже при переполнениях.
-
gbg
8 июля 2016 в 16:44
0↑
↓
В данном случае, под входным мусором понимается как алгоритм, так и данные, которые ему скормили. Так что GIGO и здесь продолжает работать.
8 июля 2016 в 14:07
0↑
↓
Wolfram Alpha считает что ответ примерно равен -4.27229×10^208 июля 2016 в 15:27
+3↑
↓
Ну прям секрет открыли. А еще unsigned int нельзя использовать для вычислений с результатом, который выходит за пределы 2^32.
Если посчитать первый множитель в обычном калькуляторе, то получается результат -7,9171109033773850490791882372801e+36, откуда видно, что для точных вычислений надо как минимум 37 разрядов.
