Cancellation and loss of precision

Let $f(x)=1/(e^x-1)$ and suppose we want to compute $f(x)$ for small, postive values of $x$. Since $e^x\approx 1$, we expect that $e^x$ has the form $1.00\cdots00abc\cdots$ for some digits. Since we have a limited precision, we expect to lose some significance, which will be magnified when we take the reciprocal. It's easy to see why this should happen. Suppose that we compute $x_1-x_2$ with $x_1 = 1.23456789$ and $x_2=1.23456787$. We assume that $x_1$ and $x_2$ represent real quantities; thus, these are only decimal approximations but all digits are correct. We might say that the numbers have $9$ digit precision. If we compute their difference, however, we obtain:

\begin{align} &1.23456789 \\ -&1.23456787 \\ &\overline{\color{red}{0.0000000}2} \end{align}

The result only has one digit precision. The problem truly becomes evident, if we take the reciprocal to get $50000000.0$. We'll show that this could be very wrong using the mpmath library in a bit. First, let's explore another example.

We can experiment with this effect using Python's decimal library. Let's load the library and set the working precision fairly low.

In [1]:
from decimal import Decimal, getcontext
getcontext().prec = 7

Let's define a small small number and evaluate our function at that point.

In [2]:
x = Decimal('0.001')
try1 = 1/(x.exp()-1)
try1
Out[2]:
Decimal('999.0010')

Another way to estimate the value of the expression is to use a Taylor approximation: $$e^x-1 = x+\frac{1}{2}x^2 + \frac{1}{6}x^3.$$ On the face of it, it doesn't seem like this should be as good. But it does avoid loss of precision. Let's try it.

In [3]:
try2 = 1/(x + x**2/2 + x**3/6)
try2
Out[3]:
Decimal('999.5002')

Noticeably different:

In [4]:
try1-try2
Out[4]:
Decimal('-0.4992')

But which is better? One approach to figure this out is to check against numpy.

In [5]:
from numpy import exp
1/(exp(0.01)-1)
Out[5]:
99.500833331945515

So, the series approximation actually looks better! When precision is really important, there's a high precision library called mpmath that let's us lift the precision way up.

In [6]:
from sympy.mpmath import mp
mp.dps = 100
x = mp.mpf('0.001')
1/(mp.exp(x)-1)
Out[6]:
mpf('999.5000833333319444444775132266865079573846929673856064295049000610275066616287755611553317493600010126')

We can use mpmath to explore our introductory example. Suppose that, to more precision, \begin{align} x_1 &= 1.23456789101112131415 \ x_2 &= 1.23456787654321012345 \end{align} Let's compute $1/(x_1-x_2)$ to high precision.

In [7]:
mp.dps = 30
x1 = mp.mpf('1.23456789101112131415')
x2 = mp.mpf('1.23456787654321012345')
1/(x1-x2)
Out[7]:
mpf('69118477.9073569268615934310542364')

Oops - very different from $1/0.00000002 = 50000000.0$!