# Explain this code to compute $\log(1+x)$

It's well known that you need to take care when writing a function to compute $\log(1+x)$ when $x$ is small. Because of floating point roundoff, $1+x$ may have less precision than $x$, which can translate to large relative error in computing $\log(1+x)$. In numerical libraries you often see a function log1p(), which computes the log of one plus its argument, for this reason.

However, I was recently puzzled by this implementation of log1p (commented source):

def log1p(x):
y = 1 + x
z = y - 1
if z == 0:
return x
else:
return x * (log(y) / z)


Why is x * (log(y) / z) a better approximation to $\log(1+x)$ than simply log(y)?

• Where did you find this implementation, was it from a library that's built into a language or from another source? – Ryan Stull Jul 27 '12 at 15:11
• they are same because from first line $x=y-1$ and this multiplied by $(log(1+x))$/(y-1) is just log(1+x),i think division helps us to minimize round-off errors or something like this – dato datuashvili Jul 27 '12 at 15:15
• @data With floating point arithmetic, after doing y=1+x and z=y-1 you do not necessarily have z == x, because of round-off error. – Chris Taylor Jul 27 '12 at 15:21
• @Ratz It's on line 224 here. I translated from C++ to Python. – Chris Taylor Jul 27 '12 at 15:26
• This is a nice example where an aggressive optimizer shouldn't be allowed to replace x * (log(y) / z) with log(1+x)! – user2468 Jul 27 '12 at 15:33

The answer involves a subtlety of floating point numbers. I'll use decimal floating point in this answer, but it applies equally well to binary floating point.

Consider decimal floating point arithmetic with four places of precision, and let

$$x = 0.0001234$$

Then under floating point addition $\oplus$, we have

$$y = 1 \oplus x = 1.0001$$ and $$z = y \ominus 1 = 0.0001000$$

If we now denote the lost precision in $x$ by $s = 0.0000234$, and the remaining part by $\bar{x}$, then we can write

$$x = \bar{x} + s$$ $$y = 1 + \bar{x}$$ $$z = \bar{x}$$

Now, the exact value of $\log(1+x)$ is

$$\log(1+x) = \log(1+\bar{x}+s) = \bar{x}+s + O(\bar{x}^2) = 0.0001234$$

If we compute $\log(y)$ then we get

$$\log(1+\bar{x}) = \bar{x} + O(\bar{x}^2) = .0001000$$

On the other hand, if we compute $x \times (\log(y)/z)$ then we get

$$(\bar{x}+s) \otimes (\log(1+ \bar{x}) \div \bar{x}) = (\bar{x}+s)(\bar{x}\div \bar{x}) = \bar{x}+s = 0.0001234$$

so we keep the digits of precision that would have been lost without this correction.

• For a brief moment I was contemplating the wrong notion that simply returning $x$ in the else block would suffice, until I reminded myself that the else block has to work not only for $x$ that is small. – Frenzy Li Jun 3 '18 at 1:36