Имеющие дело с прикладными вычислениями знают, какие неприятности может преподносить конечная точность представления вещественных чисел в ЭВМ. Наиболее известные в этом плане проблемы — это решение чувствительных к возмущениям (так называемых, плохо обусловленных) систем линейных уравнений и нахождение собственных значений несимметричных матриц.
Когда речь идет о повседневных арифметических операциях, проблемы с конечной точностью вычислений не выглядят столь пугающими. И наилучшей проверкой того, что результат получен правильно, является сравнение значений полученных на различных точностях.
Если, например, вычисления, полученные на одинарной и удвоенной точностях совпадают, то создается чувство уверенности в результате, по крайней мере с точностью сопоставимой с одинарной. Здесь, я бы хотел привести один интересный пример, демонстрирующий, что даже в сравнительно несложной арифметической задаче подобная устойчивость при переменной точности представления чисел не может служить основанием для такой уверенности.
Необходимо вычислить значение функции д��ух переменных при определенных значениях (даны ниже) ее аргументов:
Этот пример был мною замечен, когда я разбирался с библиотекой C-XSC (система классов языка C для (преимущественно) интервальных научных расчетов с произвольной точностью). Данная библиотека отлично подходит для практического исследования вычислительной устойчивости различных численных алгоритмов.
Для эмуляции вычислений с плавающей точкой на Python установим пакет mpmath. В принципе, если ограничиться точностями IEEE 754, можно было бы использовать и NumPy, но нам нужно показать, что результаты, получаемые в рамках IEEE 754 неверны.
Вычислим значение функции f(a, b) при a = 77617 и b = 33096.
Читатель может заметить, что задание точности через dps в mpmath это не совсем правильный подход при эмуляции реальной удвоенной точности. Если речь идет об удвоенной и\или одинарной точности в рамках IEEE то, пожалуй, более корректно описывать их характеристики через двоичную систему. Здесь, однако, это не так важно; зато использование mp.dps более проще интерпретируется — т.е. как количество десятичных значащих цифр в представлении числа.
Выполняя код, получим:
Вполне убедительно, не так ли? Только значение-то неправильное! Правильное значение при данных a и b вообще меньше единицы!
Дополним пример следующими вычислениями:
Получим (некоторые строки опущены):
Видно, что несмотря на устойчивость, некоторые значения все-таки существенно отличаются от обычных, что уже заставляет задуматься.
Скажу сразу, значения, получаемые при точности dps=36 и выше, правильные. Но откуда знать, что при дальнейшем увеличении точности не произойдет еще какого-либо скачка, ведь, как было видно со значением 1.17260..., даже устойчивость результата при различных точностях не может гарантировать его правильность.
К счастью, и на этот вопрос можно ответить, используя пакет mpmath. Данный пакет позволяет производить интервальные вычисления, что дает возможность оценить область возможных значений функции при представлении ее аргументов с различной точностью.
Выполним расчеты, используя аппарат интервальной арифметики mpmath:
Получим следующее:
Интересно проследить эволюцию интервала возможных значений функции: он стабилизируется только при использовании точностей от 36 и более значимых десятичных цифр, хотя и постепенно сужается.
Из интервальных расчетов становится вполне ясным, что доверять следует только результатам, получаемым при точностях от dps=36 и более десятичных цифр.
Этот пример — наглядная демонстрация того, насколько нужно быть осторожным осуществляя вычисления с плавающей точкой, и что даже 128-битная (например, numpy.float128 если говорить о Python и NumPy) точность может быть недостаточной. Он также показывает, что нельзя доверять и устойчивости результата, полученного на различных точностях. Применение аппарата интервальных вычислений может быть одним из вариантов решения в этом случае, которое позволяет оценить необходимую точность для получения адекватного результата.
Когда речь идет о повседневных арифметических операциях, проблемы с конечной точностью вычислений не выглядят столь пугающими. И наилучшей проверкой того, что результат получен правильно, является сравнение значений полученных на различных точностях.
Если, например, вычисления, полученные на одинарной и удвоенной точностях совпадают, то создается чувство уверенности в результате, по крайней мере с точностью сопоставимой с одинарной. Здесь, я бы хотел привести один интересный пример, демонстрирующий, что даже в сравнительно несложной арифметической задаче подобная устойчивость при переменной точности представления чисел не может служить основанием для такой уверенности.
Необходимо вычислить значение функции д��ух переменных при определенных значениях (даны ниже) ее аргументов:
Этот пример был мною замечен, когда я разбирался с библиотекой 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) точность может быть недостаточной. Он также показывает, что нельзя доверять и устойчивости результата, полученного на различных точностях. Применение аппарата интервальных вычислений может быть одним из вариантов решения в этом случае, которое позволяет оценить необходимую точность для получения адекватного результата.