30
$\begingroup$

Sometimes it is hard to understand how numerical expressions are evaluated. I remember reading claims by Wolfram on how smart the Kernel is to evaluate expressions trees numerically by recognizing patterns, yet I don't see how it applies in very simple examples.

This question is about numerics, but to see if symbolics can help in an elegant way.

It is known that when working with finite precision, the function $\log(1+x)$ should have a special implementation for small $x$. That is why functions like log1p exists in many libraries (on top of log). For example:

/*C code*/ log(1. + 1.e-15) == 1.11022e-15

/*C code*/ log1p(1.e-15) == 1.e-15

(The second version is more exact, the first is "wrong")

In Mathematica:

Log[1. + 1.*^-15] == 1.11022*10^-15

(wrong answer)

Mathematica doesn't have such Log1P function. One can say, well, that is because it doesn't need to, because of the symbolic power. In fact one can know the answer.

N[Log[1 + 1/10^15], 100] = 9.999999999999995000000000000003333333333333330833333333333335333333333333331666666666666668095238095*10^-16

But this is not general, if I want to evaluate Log[1 + x] and x has machine precision then I can't force to use something like log1p. Because 1+x will evaluate to a machine precision number.

These are my attempts:

x = 1.*^-15;
Log[1 + x]
Log[1. + x]
N[Log[1. + N[x, 20]], 20]
N[Log[1 + N[x, 20]], 20]

All evaluate to the wrong answer (1.11022*10^-15)

Finally I find this expression,

N[Log[1 + Rationalize[x, 0]], 20]
9.999999999999995000*10^-16

But really? Is it that hard to get $\log(1+x)$ for small $x$ numerically? Do I have to roll my own Log1p?

$\endgroup$
  • 4
    $\begingroup$ "I remember reading claims by Wolfram on how smart the Kernel is" reminded me of blog.stephenwolfram.com/2012/12/… $\endgroup$ – Dr. belisarius Dec 20 '13 at 1:16
  • 3
    $\begingroup$ We also need the inverse of log1p, expm1, for E^x - 1 when x is small. $\endgroup$ – Ray Koopman Dec 20 '13 at 5:34
  • $\begingroup$ @Ray, I suppose the trick here can be used if need be. $\endgroup$ – J. M. will be back soon Jun 7 '15 at 20:55
  • $\begingroup$ @alfC See updated answer. $\endgroup$ – Mark Adler Dec 7 '16 at 15:15
18
$\begingroup$

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]

errors all zero!


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]

things go south below 0.5

$\endgroup$
  • 1
    $\begingroup$ Yes, it helps a lot. Without compiling, running it on a million random numbers in 0..0.1, I got about six seconds with no compilation, 0.4 seconds compiled to the MVM, and 0.07 seconds compiled to C. $\endgroup$ – Mark Adler Dec 20 '13 at 7:30
  • 1
    $\begingroup$ Oh, and using arbitrary-precision numbers takes 87 seconds. $\endgroup$ – Mark Adler Dec 20 '13 at 7:37
  • 6
    $\begingroup$ See Internal`Log1p and Internal`Expm1. :) $\endgroup$ – yode Dec 7 '16 at 14:42
  • 1
    $\begingroup$ @yode When there's a standard name, I try ?*`Log1p or rather ?*`*og1* taking most of the standard name and wildcarding things that might or might not be capitalized and so forth. Sometimes, though, you get hundreds of hits and it's not so helpful. Some people do a similar thing with Names[]. $\endgroup$ – Michael E2 Dec 7 '16 at 15:30
  • 1
    $\begingroup$ Just a side note: Internal`Log1p and Internal`Expm1 are not in v9.0.1. $\endgroup$ – xzczd Nov 3 '17 at 7:48
20
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Very elegant solution. I hope you don't mind that I added a plot to illustrate the result. (Probably the same as most of other answers). $\endgroup$ – alfC Dec 20 '13 at 7:34
  • $\begingroup$ Same answer as mine, except it does needless operations for large x. $\endgroup$ – Mark Adler Dec 20 '13 at 7:36
  • $\begingroup$ @Mark Yes, I omitted a check for large x. I also used 1. instead of the 1 that would handle cases in which x has more than machine precision. As for being the same, note that w === 1. is not the same as w - 1 == 0, and Divide[a,b] is not the same as a/b for machine precision numbers. $\endgroup$ – Ray Koopman Dec 20 '13 at 8:20
  • 3
    $\begingroup$ @MarkAdler according to Daniel Lichtblau recently in comments on another answer, a/b converts into a*1/b before evaluation, whereas Divide does the division directly. This could be important when we work with machine-precision numbers, but is not documented anywhere AFAIK. $\endgroup$ – Oleksandr R. Dec 20 '13 at 8:57
  • 1
    $\begingroup$ FYI, I get very slightly better results with mine using Divide instead of /. The root-mean-square error goes down from 1.19*(10^-16) to 1.13*(10^-16). I have updated mine to use Divide. $\endgroup$ – Mark Adler Dec 21 '13 at 5:04
9
$\begingroup$

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

$\endgroup$
8
$\begingroup$

Problem

The problem with Log[1. + 1.*^-15] not yielding 1.*^-15 is not due to Log, but to MachinePrecision inputs, which I think the OP implied in the question statement:

1 + 1.*^-15
% - 1
(*
  1.
  1.11022*10^-15
*)

So Log[1 + 1.*^-15] does return the right answer, 1.11022*10^-15, for the actual input.

Solution

Here is a simple way to get log1p-type evaluation:

log1p[x_] := x Hypergeometric2F1[1, 1, 2, -x];

log1p[1.*^-15]
(*  1.*10^-15  *)

One needs to be careful, because log1p[x] evaluates to Log[1+x] when x is symbolic and you lose the precision:

log1p[x] /. x -> 1.*^-15
(*  1.11022*10^-15  *)

To prevent this, one can use the Precision tricks in the other answers, or use ?NumericQ and clear the previous definition. For instance

log1p[x] /. x -> SetPrecision[1.*^-15, $MachinePrecision] // N
(*  1.*10^-15  *)

or

ClearAll[log1p];
log1p[x_?NumericQ] := x Hypergeometric2F1[1, 1, 2, -x];

log1p[x] /. x -> 1.*^-15
(*  1.*10^-15  *)

Of course, in this second method, you lose the symbolic equivalence to Log[1+x]. But all the other current solutions suffer the same drawbacks of one or the other definitions of log1p that are given here.

Addendum: expm1

The function

expm1[x_] := x Hypergeometric1F1[1, 2, x];

can be used like log1p above.

$\endgroup$
  • 3
    $\begingroup$ This is pretty clever! I know for a fact that the hypergeometric functions in Mathematica were carefully written to be able to deal with arguments at or near $1$, and this carefulness is exploited here to compute $\log(1+x)$. For expm1, x MittagLefflerE[1, 2, x] works as well. $\endgroup$ – J. M. will be back soon Jun 8 '15 at 5:49
  • $\begingroup$ @J.M. Thanks. I wasn't familiar with MittagLefflerE. It seems slower, though. $\endgroup$ – Michael E2 Dec 11 '16 at 18:20
8
$\begingroup$
LogLogPlot[{Internal`Log1p[x], Log[1 + x]}, {x, 1*^-17, 1*^-14}]

enter image description here

ps:Of course,maybe you need Internal`Expm1,too.

$\endgroup$
  • $\begingroup$ Well-found! +1 :-) $\endgroup$ – dr.blochwave Dec 7 '16 at 14:43
  • 4
    $\begingroup$ Note Internal`Log1p and Internal`Expm1 are new in V11. $\endgroup$ – Michael E2 Dec 7 '16 at 14:54
  • $\begingroup$ @MichaelE2 I don't know it is added in V11. :) $\endgroup$ – yode Dec 7 '16 at 14:56
  • 2
    $\begingroup$ I just checked V10.4.1. (Indeed, I'm pretty sure I checked last year, when I wrote my answer.) Perhaps this Q&A made WRI think about adding them. (But IMO, they are the wrong functions to add. What one really needs is Log1p[x]/x etc., like Sinc[x].) $\endgroup$ – Michael E2 Dec 7 '16 at 14:59
  • $\begingroup$ @MichaelE2 I don't think WRI will response SE so positive.I posted so many bug post but those problem still. :) $\endgroup$ – yode Dec 7 '16 at 15:05
6
$\begingroup$

Perhaps

SetAttributes[log1p, NumericFunction];

N[log1p[x_?MachineNumberQ], _] := N@log1p[SetPrecision[x, $MachinePrecision]];

N[log1p[x_?NumericQ], {MachinePrecision, MachinePrecision}] := 
  N@N[log1p[x], $MachinePrecision];

N[log1p[x_?NumericQ], a_] := N[Log[1 + x], a]

EDIT This probably makes more sense

ClearAll[log1p]

log1p[x_?MachineNumberQ] := N@log1p[SetPrecision[x, $MachinePrecision]];
log1p[x_?InexactNumberQ] := Log[1 + x];
$\endgroup$
  • $\begingroup$ Nice defination. But the underlying Precision and Accuracy can go wrong for some cases? $\endgroup$ – Silvia Dec 20 '13 at 6:35
  • $\begingroup$ How is this function used? I get back the unevaluated log1p[1.*^-15. $\endgroup$ – alfC Dec 20 '13 at 6:51
  • $\begingroup$ @alfC, only through N, but you are right, I should have made it evaluate automatically. I am not sure why I went this road. Probably sleep depravation $\endgroup$ – Rojo Dec 20 '13 at 7:34
1
$\begingroup$
Clear[log1p]

log1p[x_, 
  n : _Integer?Positive : 
   2] := (Series[Log[1 + y], {y, 0, n}] // Normal) /. y -> x

log1p[1.0*^-15]

9.999999999999995*^-16

$\endgroup$
  • 2
    $\begingroup$ One might consider using PadeApproximant[] instead… $\endgroup$ – J. M. will be back soon Jun 8 '15 at 0:08
  • $\begingroup$ @J.M. - Thanks. $\endgroup$ – Bob Hanlon Jun 8 '15 at 0:14
1
$\begingroup$
z = 10.^-15;
N[z/(z/2 + 1), 100]

log1pcontracted[z_,n:_Integer?Positive:2]:=(2z)/(2+z+((-10-5 z) z^2)/(60+z (60+11 z)+Fold[#2[[1]]/(#2[[2]]+#1)&,0,Transpose[{Table[-(324+m (1224+m (1812+m (1312+m (464+64 m)))))z^4/(5+4m),{m,n,0,-1}],Table[(1260+z (1260+154 z)+m (2288+z (2288+284 z)+m (1344+z (1344+168 z)+m (256+z (256+32 z)))))/(5+4m),{m,n,0,-1}]}]]));
log1p[z_,n:_Integer?Positive:2]:=2z/(2+z+Fold[#2[[1]]/(#2[[2]]+#1)&,0,Transpose[{Table[-(m-1)^2z^2,{m,n,2,-1}],Table[(2m-1)(2+z),{m,n,2,-1}]}]]);
log1pcontractedcompiled=Compile[{n,z},(2z)/(2+z+((-10-5 z) z^2)/(60+z (60+11 z)+Fold[#2[[1]]/(#2[[2]]+#1)&,0,Transpose[{Table[-(324+m (1224+m (1812+m (1312+m (464+64 m)))))z^4/(5+4m),{m,n,0,-1}],Table[(1260+z (1260+154 z)+m (2288+z (2288+284 z)+m (1344+z (1344+168 z)+m (256+z (256+32 z)))))/(5+4m),{m,n,0,-1}]}]]))];
log1pcompiled=Compile[{n,z},2z/(2+z+Fold[#2[[1]]/(#2[[2]]+#1)&,0,Transpose[{Table[-(m-1)^2z^2,{m,n,2,-1}],Table[(2m-1)(2+z),{m,n,2,-1}]}]])];

N[log1p[z], 100]

N[log1pcontracted[z], 100]

N[log1pcompiled[2, z], 100]

N[log1pcontractedcompiled[2, z], 100]
$\endgroup$
  • $\begingroup$ Is this the continued fraction? If so, consider using a forward evaluation method so as not to waste effort nor need to have to specify a truncation. $\endgroup$ – J. M. will be back soon Jul 18 '15 at 17:32
  • $\begingroup$ yes, continued fraction $\endgroup$ – Andreas Lauschke Jul 18 '15 at 21:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.