# Numerically approximating limit with large numbers, bumping up against machine precision

Consider $\lim_{x\rightarrow\infty}(1+\frac1x)^{x}=e$. Using standard calculus techniques, this limit can be evaluated, however, approximating it directly with numerical code can be difficult depending on machine precision. For example, here is a graph from octave:

Near the right edge of the graph, before machine precision kills the computation completely, we get $\left(1+\frac{1}{{9(10)^{15}}}\right)^{9(10)^{15}}=7.37725371726828$ which is clearly incorrect/inexact. This stems from the approximation $1+\frac{1}{{9(10)^{15}}}=\mathtt{3ff0000000000001}_{16}=1+2^{-52}$ (one plus machine epsilon), again which is obviously inexact.

Is there a way to get an approximate value of $e$ out of this formula numerically, say to more than 8 or so decimal places (in double precision float)? The accuracy peaks around $x=10^8$ with $(1+1/x)^x-e\approx -3(10)^{-8}$. Are there general numerical techniques for evaluating this or similar limits where large and small numbers are combined?

• Interesting that the accuracy maximizes very close to $\sqrt{\epsilon}$, with an argument at $~\sqrt{\epsilon}$. I'll be you can prove that this formula does indeed have a lowest relative error of $\sqrt{\epsilon}$. – user14717 May 8 '18 at 5:26
• Actually, I have a pretty convincing argument that $\sqrt{\epsilon}$ is the best you can do with this formula. Let $x = 2^{k}$ for $k \le 52$ and then do exponentiation by squaring and you should be able to see it. – user14717 May 8 '18 at 5:32
• @user14717 interesting observation. I bet you are correct, and that makes sense. – jdods May 8 '18 at 5:38

Yes. Richardson's techniques can be applied and we can almost recover the douple precision representation of $e$.

Let $$\phi(n) = (1 + 1/n)^n$$ and compute $$a_k = \phi(2^k), \quad k=1,2,3,\dotsc.$$ Here are the first 30 values of $a_k$ and some auxiliary numbers.

 k |          approximation |     fraction |   error estimate
1 |  2.250000000000000e+00 |   0.00000000 |   0.00000000e+00
2 |  2.441406250000000e+00 |   0.00000000 |   1.91406250e-01
3 |  2.565784513950348e+00 |   1.53890434 |   1.24378264e-01
4 |  2.637928497366600e+00 |   1.72402823 |   7.21439834e-02
5 |  2.676990129378183e+00 |   1.84692702 |   3.90616320e-02
6 |  2.697344952565099e+00 |   1.91903568 |   2.03548232e-02
7 |  2.707739019688020e+00 |   1.95831169 |   1.03940671e-02
8 |  2.712991624253434e+00 |   1.97884059 |   5.25260457e-03
9 |  2.715632000168991e+00 |   1.98933967 |   2.64037592e-03
10 |  2.716955729466436e+00 |   1.99464945 |   1.32372930e-03
11 |  2.717618482336880e+00 |   1.99731960 |   6.62752870e-04
12 |  2.717950081189666e+00 |   1.99865851 |   3.31598853e-04
13 |  2.718115936265797e+00 |   1.99932894 |   1.65855076e-04
14 |  2.718198877721971e+00 |   1.99966439 |   8.29414562e-05
15 |  2.718240351930294e+00 |   1.99983217 |   4.14742083e-05
16 |  2.718261089904603e+00 |   1.99991608 |   2.07379743e-05
17 |  2.718271459109306e+00 |   1.99995804 |   1.03692047e-05
18 |  2.718276643766046e+00 |   1.99997902 |   5.18465674e-06
19 |  2.718279236108013e+00 |   1.99998951 |   2.59234197e-06
20 |  2.718280532282396e+00 |   1.99999475 |   1.29617438e-06
21 |  2.718281180370437e+00 |   1.99999738 |   6.48088041e-07
22 |  2.718281504414670e+00 |   1.99999869 |   3.24044233e-07
23 |  2.718281666436840e+00 |   1.99999934 |   1.62022170e-07
24 |  2.718281747447938e+00 |   1.99999967 |   8.10110983e-08
25 |  2.718281787953491e+00 |   1.99999985 |   4.05055522e-08
26 |  2.718281808206267e+00 |   1.99999993 |   2.02527768e-08
27 |  2.718281818332656e+00 |   1.99999987 |   1.01263891e-08
28 |  2.718281823395851e+00 |   2.00000009 |   5.06319431e-09
29 |  2.718281825927448e+00 |   2.00000000 |   2.53159715e-09
30 |  2.718281827193247e+00 |   2.00000000 |   1.26579858e-09


The second column contains $a_k$, i.e. the approximation of $e$. The third column contains Richardson's fraction, i.e., $$f_k = \frac{a_{k-1} - a_{k-2}}{a_{k}-a_{k-1}}$$ In exact arithmetic it would converge to $2$ from below. The convergence would be monotonic for sufficient large $k$ with $2-f_k = O(2^{-k})$. This pattern is observed until $k=26$ and it is clearly broken for $k=27$. The fourth column contains Richardson's error estimate, i.e, $$e_k = a_k - a_{k-1}$$ We can improve the accuracy of our approximation $a_k$ by adding the error estimate $e_k$, i.e., $$b_k = a_{k+1} + e_{k+1}, \quad k=1,2,3,\dotsc$$ This gives us the numbers

 k |          approximation |     fraction |   error estimate
1 |  2.632812500000000e+00 |   0.00000000 |   0.00000000e+00
2 |  2.690162777900696e+00 |   0.00000000 |   1.91167593e-02
3 |  2.710072480782852e+00 |   2.88051902 |   6.63656763e-03
4 |  2.716051761389766e+00 |   3.32978233 |   1.99309354e-03
5 |  2.717699775752015e+00 |   3.62817263 |   5.49338121e-04
6 |  2.718133086810941e+00 |   3.80330557 |   1.44437020e-04
7 |  2.718244228818849e+00 |   3.89871541 |   3.70473360e-05
8 |  2.718272376084548e+00 |   3.94858986 |   9.38242190e-06
9 |  2.718279458763880e+00 |   3.97409856 |   2.36089311e-06
10 |  2.718281235207324e+00 |   3.98699962 |   5.92147814e-07
11 |  2.718281680042452e+00 |   3.99348732 |   1.48278376e-07
12 |  2.718281791341928e+00 |   3.99674053 |   3.70998254e-08
13 |  2.718281819178145e+00 |   3.99836941 |   9.27873881e-09
14 |  2.718281826138617e+00 |   3.99918481 |   2.32015755e-09
15 |  2.718281827878913e+00 |   3.99959248 |   5.80098488e-10
16 |  2.718281828314009e+00 |   3.99979280 |   1.45032134e-10
17 |  2.718281828422786e+00 |   3.99989385 |   3.62589958e-11
18 |  2.718281828449980e+00 |   4.00006532 |   9.06460092e-12
19 |  2.718281828456779e+00 |   3.99967342 |   2.26633527e-12
20 |  2.718281828458479e+00 |   3.99947753 |   5.66657832e-13
21 |  2.718281828458903e+00 |   4.00837696 |   1.41368398e-13
22 |  2.718281828459010e+00 |   3.96265560 |   3.56751665e-14
23 |  2.718281828459037e+00 |   4.01666667 |   8.88178420e-15
24 |  2.718281828459043e+00 |   4.28571429 |   2.07241631e-15
25 |  2.718281828459044e+00 |   4.66666667 |   4.44089210e-16
26 |  2.718281828459046e+00 |   1.00000000 |   4.44089210e-16
27 |  2.718281828459045e+00 |  -3.00000000 |  -1.48029737e-16
28 |  2.718281828459045e+00 |         -Inf |   0.00000000e+00
29 |  2.718281828459045e+00 |          NaN |   0.00000000e+00


In exact arithmetic, Richardson's fraction would converge toward $4$ from below and monotonically so for sufficient large $k$. This pattern is observed for the computed values until $k=17$. It is broken at $k=18$. This does not imply that the error estimates, i.e., $$e_k' = \frac{b_k - b_{k-1}}{3}$$ are unreliable for $k>17$, but the accuracy of the error estimates suffers as $k$ is increased.

What can we do? We are lucky, and $$b_{17} + e_{17}' \approx 2.718281828459045$$ differs from the floating point representation of $$e \approx 2.718281828459046$$ by one unit in the last place. It is worth noting that these numbers required $a_k = \phi(2^k)$ for $k \leq 18$.

Underlying the use of Richardson's techniques is the existence of an asymptotic error expansion of the form $$e - \phi(n) = c_1 n^{-1} + c_2 n^{-2} + \dotsc.$$ The properties of Richardson's fractions and estimates flows directly from this expansion. When adding Richardson's error estimate to the current approximation you typically reduce the error, but you lose control, unless you compute a new error estimate. You gradually lose the ability to do this as the order of the approxmation increases and you run out of bits.

• This is great. I don't think I've ever heard of Richardson's fractions. It seems just sticking with powers of two corrects most machine precision errors at least up until machine epsilon. – jdods May 8 '18 at 13:50
• No worries, I am not aware of single textbook which investigates when Richardson's error estimates are reliable. The fractions have a provable behavior in exact arithmetic. Deviations from the correct behavior indicates that rounding errors have become significant and that the error estimates are loosing quality. This is particularly important when computing integrals or solving ODEs as there is little point in reducing the stepsize when rounding errors start to dominate. – Carl Christian May 8 '18 at 13:56