Имеющие дело с прикладными вычислениями знают, какие неприятности может преподносить конечная точность представления вещественных чисел в ЭВМ. Наиболее известные в этом плане проблемы — это решение чувствительных к возмущениям (так называемых, плохо обусловленных) систем линейных уравнений и нахождение собственных значений несимметричных матриц.
Когда речь идет о повседневных арифметических операциях, проблемы с конечной точностью вычислений не выглядят столь пугающими. И наилучшей проверкой того, что результат получен правильно, является сравнение значений полученных на различных точностях.
Если, например, вычисления, полученные на одинарной и удвоенной точностях совпадают, то создается чувство уверенности в результате, по крайней мере с точностью сопоставимой с одинарной. Здесь, я бы хотел привести один интересный пример, демонстрирующий, что даже в сравнительно несложной арифметической задаче подобная устойчивость при переменной точности представления чисел не может служить основанием для такой уверенности.
Необходимо вычислить значение функции двух переменных при определенных значениях (даны ниже) ее аргументов:
Этот пример был мною замечен, когда я разбирался с библиотекой 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) точность может быть недостаточной. Он также показывает, что нельзя доверять и устойчивости результата, полученного на различных точностях. Применение аппарата интервальных вычислений может быть одним из вариантов решения в этом случае, которое позволяет оценить необходимую точность для получения адекватного результата.
Автор: scidam
