Натуральный логарифм с использованием ряда в Haskell дает неточные результаты
Поэтому я написал две функции для вычисления натурального логарифма переменной x, после увеличения верхней границы добавочной суммы до 33000 функции все еще возвращают неточные результаты, тестируемые в ghci, по сравнению с функцией журнала по умолчанию, импортированной из Prelude, вот определение кода:
lnOfx :: Float -> Float
lnOfx x = netSum f 1 33000
where f i = 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm
where oddTerm = 2*i-1
netSum f minB maxB = sum [f i | i <- [minB .. maxB]]
lnOfx2 :: Float -> Float
lnOfx2 x = netSum f 1 33000
where f i = (1/i)*((x-1)/x)**i
netSum f minB maxB = sum [f i | i <- [minB .. maxB]]
И результаты теста:
log 3
1.0986122886681098
lnOfx 3
1.0986125
lnOfx2 3
1.0986122
log 2
0.6931471805599453
lnOfx 2
0.6931472
lnOfx2 2
0.6931473
Так почему же результаты отличаются и как правильно (как и функция log из Prelude) правильно рассчитать натуральный логарифм?
1 ответ
Математика с плавающей точкой сложна. Одна из вещей, которая может вызвать потерю точности, - это сложение чисел с очень разными величинами. Например, в вашем алгоритме, начиная примерно i=25
условия в вашей сумме достаточно малы, чтобы они перестали иметь значение:
-- 25t term:
let f x i = let oddTerm = 2*i-1 in 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm
let y = f 3.0 25
-- summation up to 24 item
let s = 1.098612288668109
-- this will return True, surprisingly!
s + y == s
Одна вещь, которую вы можете сделать, чтобы смягчить это, это добавить числа в обратном порядке, чтобы маленькие числа складывались вместе, прежде чем они будут добавлены к большим числам.
lnOfx :: Float -> Float
lnOfx x = netSum f 1 33000
where f i = 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm
where oddTerm = 2*i-1
netSum f minB maxB = sum $ reverse [f i | i <- [minB .. maxB]]
В моих тестах этого было достаточно, чтобы print (lnOfx 3.0)
а также print (log 3.0)
показал все те же цифры.
Но в целом я бы порекомендовал прочитать книгу по численному анализу, чтобы узнать больше об этой проблеме.