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:

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 just 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
@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
@Chris: Your answer seems plausible to me; why did you delete it? – joriki Jul 27 '12 at 15:23
@joriki To give other people a chance to answer! I thought it was an interesting question and I had fun answering it - I didn't want to deny other people that fun. I'll undelete it in a day or two if no one else responds. – 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.

-