Elegant high precision `log1p`?

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

Update:

yode points out in a comment below that there are log1p() and expm1() functions in Mathematica! However they are hidden. They are simply:

Internal`Log1p
Internal`Expm1

They operate only on numeric inputs.

(I don't know at what version these showed up.)

Compare this plot to the one further below using the Log function. The errors really are all zero!

Plot[Internal`Log1p[x]/x - 
  N[N[Log[1 + SetPrecision[x, Infinity]], 34]]/x, {x, 0.01, 1}, 
 PlotPoints -> 100, PlotRange -> All, ImageSize -> Large]

---

I don't think Mathematica has that function. Seems like it should. (Same for expm1().) You should not need to resort to non-machine arithmetic to get the right answer.

Here is something that will do the trick using only machine arithmetic, if the input is a machine number:

log1p[x_] :=
 If[MachineNumberQ[x],
  If[x < 0.5,
   If[# - 1 == 0, x,
      x Divide[Log[#], # - 1]] &[1 + x],
   Log[1 + x]],
  Log[1 + x]]

If it's not a machine number, then it just gives you Log[1+x].

You can compile the machine number portion for much faster execution (by two orders of magnitude!):

log1px = Compile[{{x, _Real}},
  If[x < 0.5,
   If[# - 1 == 0, x,
      x Divide[Log[#], # - 1]] &[1 + x],
   Log[1 + x]], CompilationTarget -> "C"]

The cutoff for x >= 0.5 is where Log[1+x] works just fine, with only the least significant bit varying:

Plot[Log[1 + x]/x - 
  N[N[Log[1 + SetPrecision[x, Infinity]], 34]]/x, {x, 0.01, 1}, 
 PlotPoints -> 100, PlotRange -> All, ImageSize -> Large]

I think that is because 1. and 1.*^-15 are machine-precision numbers, and Mathematica does NOT do precision-tracking on machine-precision calculations. I suggust using the arbitrary-precision numbers instead.

Through[{Precision, Accuracy}@#] & /@ {1., 1.*^-15}

{{MachinePrecision, 15.9546}, {MachinePrecision, 30.9546}}

Now we specify a precision to the numbers to make them arbitrary-precision numbers:

res = Log[1.`20 + 1.`20*^-15]

1.000*10^-15

Through[{Precision, Accuracy}[res]]

{5., 20.}

And it won't be wrong even for low precision cases:

res2 = Log[1.`5 + 1.`5*^-15]

0.*10^-5

Through[{Precision, Accuracy}[res2]]

{0., 5.}