# Computing a large exp(x) in a numerically robust way.

I'm trying to compute $\lfloor e^x \rfloor$, where x is a 64-bit integer. The problem is that the result of the computation may be close to 2^64. In this range, 64-bit floating point numbers will be sparser than 64-bit integers, so it would be a bad idea to use something like the exp library function in C, which returns a double. Instead I'd like to use a method which computes the 64-bit integer result directly.

Is there a formula or well-conditioned algorithm for computing this floor value as an integer, without losing precision by going through floating point?

-
If $x$ is a 64-bit integer, then $\lfloor e^x \rfloor$ is as large as $e^{2^{64}}$, which has approximately $8 \times 10^{18}$ digits. That's about 2.8 exabytes of information. –  heropup Feb 12 '14 at 1:11
@heropup x is not that large. It's as much as $ln(2^{64}-1)$. –  Brian Gordon Feb 12 '14 at 1:12
@heropup the result, $e^x$, may be as large as $2^{64}$. Not $x$ itself –  Omnomnomnom Feb 12 '14 at 1:12
@BrianGordon Given that $x$ has to be less than roughly 40 (plus/minus), why not simply build a table with e.g. Wolfram Alpha and store that in your app? –  Steven Stadnicki Feb 12 '14 at 1:13
@StevenStadnicki That was precisely what I was about to suggest, given the condition that $e^x < 2^{64}$. I don't know why $x$ needs to be specified as a 64-bit integer. It is the result that needs that specification. –  heropup Feb 12 '14 at 1:14

By far the easiest way to do this is just to set your favorite CA system (say, Wolfram Alpha) to computing the precise value of $e^x$ for the $\approx 40$ values of $x$ s.t. $e^x\lt 2^{64}$, and store the results in a table that's hard-coded into your app.

-

First thing that comes to mind:

Use the series representation $e^x = \sum_{i=1}^\infty \frac{x^n}{n!}$. Because $21!>2^{64}$, calculating

After some more thorough consideration it seems that in order for $$\left\lfloor e^x \right\rfloor = \left\lfloor \sum_{n=0}^{n} \frac{x^n}{n!} \right\rfloor$$ to produce an exact value, we need $n = O(x)$. That is, the higher $x$ is, the more terms we need, so that $n \geq 106$ would be required for precision up to $x = 40$. This is not as useful as I had guessed.

-
This doesn't allow for the possibility that later terms might 'conspire' to push the floor value up a little bit (from ...3516.999999... to ...3517.000001, for instance), though obviously you can bound the contribution from the later terms and confirm that this is impossible in specific cases. And this will require bignum intermediate storage. –  Steven Stadnicki Feb 12 '14 at 1:15
The given formula doesn't even work. First, the index of summation doesn't match. Second, the floor of the sum is not equal to the sum of the floors. You can't stop at $21!$. If you evaluate the LHS for $x \ge 9$, the result is not equal to the RHS: $$\lfloor e^9 \rfloor - \left\lfloor \sum_{n=0}^{20} \frac{9^n}{n!} \right\rfloor = 8103 - 8099 \ne 0.$$ The discrepancy gets rapidly worse as $x$ increases. –  heropup Feb 12 '14 at 1:27
@heropup thank you for noting the typo. Note also that I said to take the floor of the sum, not the sum of the floors, which I agree would yield the incorrect result. –  Omnomnomnom Feb 12 '14 at 1:32
Yes, I did notice it. Your answer is still mathematically incorrect. Try computing the RHS and LHS of your equation for $x = 30$. They aren't even close to equal. The reason why I mentioned the sum of the floors is because you are assuming that the remaining terms beyond $n > 20$ are too small to contribute to the value of the sum; i.e., you are implicitly assuming the sum of the floors equals the floor of the sum, thereby underestimating the true value. –  heropup Feb 12 '14 at 1:33
@heropup well, that's embarrassing. Thanks for pointing that out. –  Omnomnomnom Feb 12 '14 at 1:45
Table[Floor[Exp[k]], {k, 0, 44}]


gives

$$\{1,2,7,20,54,148,403,1096,2980,8103,22026,59874,162754,442413,1202604,3269017,8886110,24154952,65659969,178482300,485165195,1318815734,3584912846,9744803446,26489122129, 72004899337,195729609428,532048240601,1446257064291,3931334297144,10686474581524,29048849665247,78962960182680,214643579785916,583461742527454,1586013452313430,4311231 547115195,11719142372802611,31855931757113756,86593400423993746,235385266837019985,639843493530054949,1739274941520501047,4727839468229346561,12851600114359308275\}$$

-

Perhaps some sort of CORDIC algorithm is the best option. Computing huge factorials and fractions to the required precision will be painful.

If the performance is not so critical, and programmer time is costly, perhaps using a high-precision floating point library, like MPFR, is the best bet.

-