В отличие от того, что было продемонстрировано на Chux's и
aka.nice свои ответы и упоминал в разных комментариях я
полагаю, что мы сделали поспешный вывод, сосредоточившись только на одном вычислении,
то есть sin(5). Даже если ответ ОП ближе к истинному значению,
Надо сказать, что как алгоритм ОП и факториал
алгоритм столь же плох для вычисления sin(5), но в конце концов,
Факториал алгоритм лучше.
Этот ответ не будет вдаваться в подробности о Floating-Point
Арифметика. Превосходную монографию можно найти по этой теме в Что
Каждый ученый должен знать о с плавающей точкой
Арифметика . Также я не открою банку с червями Средний
Точность с плавающей точкой .
отказ от ответственности: Я не хочу утверждать, что то, что я здесь пишу, является на 100% правильным, и я точно не авторитет в этой области. Я просто нашел вопрос чрезвычайно интересным и попытался извлечь из него как можно больше. Я, конечно, приветствую любой комментарий!
некоторые интересные чтения, которые приводят к этому:
Некоторые исходные понятия
представление с плавающей точкой (FPR): в конечной плавающей точке
представление, число записывается как ± 0.d 1 d 2
... d p x b e . Он имеет точность p и четность
base b где все цифры d i являются целыми числами с 0 & le;
d i . Это просто, что такое представление
не может представлять все действительные числа. Например, если b = 10 и p =
3 , число 1/3 приблизительно равно 0,333 x 10 0 или если
b = 2 и p = 24 нужно приблизительно 0,1 как
0,110011001100110011001101 x 2 5 .
абсолютная ошибка (AE): , если действительное число x имеет AE, равное
& Delta; x , тогда мы знаем, что ожидаемое истинное значение x̃ из x
лежит между x - & Delta; x & le; x̃ & le; х + & Delta; х . Если x
ближайший FPR к истинному значению x̃ , то мы знаем, что его AE
is ((b / 2) b -p - 1 ) x b e . Например, если b
= 10 и p = 3 , число 1/3 приблизительно равно 0,333 x
10 0 и имеет абсолютную ошибку 0,0005 x 10 0
указывая, что 0,3325 x 10 0 & le; 1/3 & le; 0,3335 х
10 0
Что все это значит для серии Тейлор-Маклаурин?
Для стандартного IEEE binary64 (двойная точность) вычисления
( b = 2 и p = 53 ) из серии Тейлора-Маклаурина sin(5) один
быстро видит, что лошадь уже сбежала на 2-й и 3-й
срок. Здесь FPR является точным только до 2^-49, как можно видеть в
таблица ниже (представляющая ближайшее двоичное представление
истинные дроби):
j term FPR e AE
0.12345678901234567
1 5^1/1! 5.00000000000000000E+00 3 2^-51 4.4E-16
3 -5^3/3! -2.08333333333333321E+01 5 2^-49 1.8E-15
5 5^5/5! 2.60416666666666679E+01 5 2^-49 1.8E-15
7 -5^7/7! -1.55009920634920633E+01 4 2^-50 8.9E-16
j sum FPR e AE
0.12345678901234567
1 5 5.00000000000000000E+00 3 2^-51 4.4E-16
3 -95/6 -1.58333333333333339E+01 4 2^-50 8.9E-16
5 245/24 1.02083333333333339E+01 4 2^-50 8.9E-16
7 -5335/1008 -5.29265873015873023E+00 3 2^-51 4.4E-16
Точность до 2^-49 можно понять следующим образом. Если
вы смотрите на термин 5^3/3!, то ближайший FPR этого числа
фракция (5864062014805333 / 9007199254740992) * 2^5. Как и ты
видите, нам здесь не хватает части, а именно
5^3/3! - 5864062014805333 / 9007199254740992) * 2^5
= 1/844424930131968
~ 1.1842378929335002E-15 < 2^-49
Так что независимо от того, что мы пытаемся сделать, мы всегда пропускаем эту часть. так что это не так
можно вычислить sin(5) с точностью лучше 2^-49. Как
это точность наибольшего термина (абсолютное значение)
прибавил к сумме. Тем не менее, другие термины также вносят ошибки и все
эти ошибки накапливаются. Итак, после первых 30 семестров, мы знаем, что
AE для sin(5) накапливается так:
2^-51 + 2^-49 + 2^-49 + 2^50 + ... = 5.45594...E-15
Это число, однако, слишком идеалистично, поскольку больше потерь
писправление из-за потери значимости .
Какой алгоритм лучше?
Два представленных алгоритма (при условии, что все переменные IEEE binary64 (двойная точность) переменные, за исключением i):
алгоритм 1: Это слегка принятая версия алгоритма, представленная здесь
fact = 1.0_dp; pow = x; sum = x; xx = x*x
do i=2,30
j=2*i-1
fact = fact*j*(j-1)
pow = -pow * xx
term = pow/fact
sum = sum + term
end do
алгоритм 2: Это слегка принятый вариант алгоритма, представленный здесь
term = x; sum = x; xx = x*x
do i=2,30
j=2*i-1
term = (-term*xx)/(j*(j-1))
sum = sum + term
end do
Хотя оба алгоритма математически одинаковы, существует небольшая разница в численности.В операциях с плавающей запятой деления и умножения считаются безопасными, то есть они всегда приводят к ближайшему FPR истинного результата.То есть с заданным входом.Это, однако, не означает, что множественные деления и умножение приведут к ближайшему FPR полного вычисления.
Именно поэтому алгоритм 2 немного хуже, чем алгоритм 1 .Вычисление term = (-term*xx)/(j*(j-1)) содержит несколько умножений и делений и уже использует приблизительную версию term.Это в отличие от алгоритма 1 , где term =
pow/fact - это одиночная операция.
В таблице ниже показаны различия в term, начиная с j=5:
j term(alg-1) term(alg-2)
0.12345678901234567 0.12345678901234567
1 5.00000000000000000E+00 5.00000000000000000E+00
3 -2.08333333333333321E+01 -2.08333333333333321E+01
5 2.60416666666666679E+01 2.60416666666666643E+01
7 -1.55009920634920633E+01 -1.55009920634920633E+01
9 5.38228891093474449E+00 5.38228891093474360E+00
11 -1.22324747975789649E+00 -1.22324747975789627E+00
13 1.96033249961201361E-01 1.96033249961201306E-01
15 -2.33372916620477808E-02 -2.33372916620477738E-02
17 2.14497166011468543E-03 2.14497166011468499E-03
Размер создаваемой ошибки имеет величину:
AE(term(j-2)) * xx / (j*(j-1))
Можем ли мы еще улучшить?
Ответ однозначно да, но нам нужноиспользовать некоторые «трюки».Самый простой способ - использовать более высокую точность, а затем преобразовать все в двойную точность, но это обман!
Ранее упоминалось, что потеря значимости приведет к большему количеству ошибок ив этом комментарии он упоминается начиная с 23! и далее, вы теряете точность в факториале, так как он больше не может быть корректно представлен в виде двоичного числа64.
Мы попробуемчтобы улучшить ответ, отслеживая термин ошибки и используя так называемое компенсированное суммирование или суммирование Кахана .Как указывалось в начале этого поста, обратный факториал по существу приводит к потере точности.
вычисляя обратный факториал: вместо вычисления факториала мы вычислимобратный факториал следующим образом.Представьте, что мы вычислили f ~ 1 / (n-1)! с соответствующим значением ошибки q таким, что f + q ~ 1 / (n-1)! Мы знаем, что «FPR» f можно записать следующим образом: f = a * b ep с a, b, e, p целые числа.Таким образом, можно использовать целочисленное деление для вычисления 1 / n! и соответствующего ему члена ошибки, используя следующий набор операций:
f = (a/n) * b^(e-p)
q = (q + mod(a,n) * b^(e-p))/n
в Фортране, что приводит к следующему коду:
f = 1.0_dp; q = 0.0_dp
do i = 2,10
! retrieving the integer from 1/(i-1)!
a = int(scale(fraction(f),digits(f)),kind=INT64)
! computing the error on f while keeping track of the
! previous error
q = (q + scale(real(mod(a,i),kind=dp),exponent(f)-digits(f)))/i
! computing f/i resembling 1/i!
f = scale(real(a/i ,kind=dp),exponent(f)-digits(f))
! rebalancing the terms
t = f; f = f + q; q = q - (f - t)
end do
Здесь последняя строка перебалансирует f с ошибкой q.Трюк, вытекающий из суммирования Кахана .
Суммирование Кахана с обратным факториалом ... новый sin:
Теперь это возможнообъединить этот трюк с алгоритмом 2 , чтобы сделать его немного лучше:
pow = x; sum = x; xx = x*x; err = 0.0_dp;
f = 1.0_dp; q = 0.0_dp
do i=2,30
j=2*i-1
! The factorial part
a = int(scale(fraction(f),digits(f)),kind=INT64)
q = (q + scale(real(mod(a,j*(j-1)),kind=dp),exponent(f)-digits(f)))/(j*(j-1))
f = scale(real(a/(j*(j-1)) ,kind=dp),exponent(f)-digits(f))
t = f; f = f + q; q = q - (f - t)
pow = -pow*xx ! computes x^j
! Kahan summation
t = pow*f; err = err + (t - ((sum + t) - sum)); sum = sum + t
t = pow*q; err = err + (t - ((sum + t) - sum)); sum = sum + t
t = sum; sum = sum + err; err = err - (sum - t)
end do
Этот алгоритм приводит к удивительно хорошим результатам.Отличная точность в квадранте 1 и 2 и немного хуже в третьем квадранте.Четвертый квадрант все еще в порядке.Однако он имеет ужасную точность вблизи Pi и 2Pi.