def logarithm_series

stelliferous@0.0.2.0 · math_fp32.bend:146

import Base

146# For s near r/(2+r), d = r-(2+r)*s gives the exact real identity
147# log(1+r) = 2*atanh(s) + log(1+d/(1+s)). The odd series ends at s^11.
148# Split r and s into 12-bit heads so their head product is exact; retain
149# the low products in d. (1-s)*(1+s^2+s^4) approximates 1/(1+s).
150# FastTwoSum retains the low bits when adding the exponent's ln(2) part.
151def logarithm_series(+r: F32,+n: F32) -> F32:
152  +s = (r / (2.0 + r : F32) : F32)
153  +square = (s * s : F32)
154  +r_high = from_bits(U32.and(F32.bits(r),4294963200))
155  +r_low = (r - r_high : F32)
156  +s_high = from_bits(U32.and(F32.bits(s),4294963200))
157  +s_low = (s - s_high : F32)
158  +product_low = ((r_high * s_low + r_low * s_high : F32) + r_low * s_low : F32)
159  +correction = (((r - 2.0 * s : F32) - r_high * s_high : F32) - product_low : F32)
160  +p = ((((0.09090909090909091 * square + 0.1111111111111111 : F32) * square + 0.14285714285714285 : F32) * square + 0.2 : F32) * square + 0.3333333333333333 : F32)
161  +main = (2.0 * s : F32)
162  +rest = ((main * square : F32) * p + correction * ((1.0 - s : F32) * (1.0 + square * (1.0 + square : F32) : F32) : F32) : F32)
163  +high = (n * 0.693359375 : F32)
164  +sum = (high + main : F32)
165  +tail = (main - (sum - high : F32) : F32)
166  (sum + (rest + (n * F32.neg(0.000212194440) + tail : F32) : F32) : F32)