Return to Answer

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}]

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}]

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.,3.25,.25]*2^-52 against the value returned by setting the precision to 25, 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}]

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}]

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., x, Divide[x*Log@w, w-1.]]]

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

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.,3.25,.25]*2^-52 against the value returned by setting the precision to 25, 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}]