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

965 Views Asked by At

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?

4

There are 4 best solutions below

0
On BEST ANSWER

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.

6
On

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.

0
On
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\}$$

0
On

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.