Skip to main content
6 of 6
added 49 characters in body

David Goldberg ("What every computer scientist should know about floating-point arithmetic", ACM Computing Surveys, Vol 23, No 1, March 1991, p 12, Th 4) gives pseudocode that is equivalent to

log1p[x_Real] := With[{w = 1 + x}, If[w - 1 == 0, x, x * Log @ w/(w - 1)]]

EDIT - Following Mark Adler's comments, I checked the binary representation of the results (using RealDigits[#,2,53]) for x in Range[1.,5.,.25]*2^-52 against the value returned by setting the precision to 35, and he is right on both counts: the comparison should be w-1 == 0, not w === 1., and the division should use the a/b form, not the Divide[a,b] form. I have changed the code accordingly.


LogLogPlot[{log1p[x], Log[1 + x]}, {x, 1*^-17, 1*^-14}]

log1p

Ray Koopman
  • 3.3k
  • 15
  • 13