1

I'm trying to compute

$$e^{-a} \times \sum_{n=0}^N \frac{a^n}{n!}$$

where $e^{1}=2.7183$, N is an integer > 0, and a is a real value > 0. For example, N=a=1500.

But for large n, the n! goes to infinity, and for large a, the numerator a^n goes to infinity. Can any tricks be applied to compute this? What about in the case of either a>n or a< n?

UPDATE

Based on @DavidK answer, could I get extra precision by taking the log of the expression? For example, the natural log of the above expression is:

$$-a + log(\sum_{n=0}^N \frac{a^n}{n!})$$

then when I get a result from this expression, take its inverse natural log? If so, the problem boils down how to compute

$$log(\sum_{n=0}^N \frac{a^n}{n!})$$

without overflow (any ideas?).

user46688
  • 189
  • 5
    Don't compute $a^n$ and $n!$ separately. Compute the terms as in $({a \over n} )( {a \over n-1}) \cdots ({a \over 1})$. – copper.hat Oct 23 '15 at 20:54

3 Answers3

1

Note that $\sum_{n=0}^N {a^n \over n!} = 1+ {a \over 1} (1+{a\over 2}(1+\cdots{a \over N-1}(1+ {a \over N}))\cdots))$.

copper.hat
  • 172,524
  • Thanks, your comment above (in the original posting) seems more intuitive for me to program in a loop for each n, then sum over n=0 to N. – user46688 Oct 23 '15 at 21:07
  • Thanks @copper.hat, I'm running into the problem where the summation exceeds the max double precision value. Since a is > 1, I tried moving the exp(-a) into the summation to reduce the terms that get summed, which helped a little, but still I'm exceeding the limit. Could I do a log somewhere to complete the math then do an anti-log at the end? – user46688 Oct 23 '15 at 22:02
  • @user46688 Unless the double-precision representation of exp(-a) is zero or denormal, exp(a) should be representable in double precision, and so should all the terms in your sum, which is just a partial sum of a series representation of $e^a$. Note that multiplication by a denormal results in loss of precision, so you really should not be doing this calculation if all the terms can't be represented in your chosen floating-point precision. Perhaps you should consider quadruple precision. – David K Oct 23 '15 at 22:17
  • @DavidK, yes exp(-a) is zero when the summation goes to infinity. I'm coding in Javascript. I haven't heard the term denormal before... are you referring to the log discussion? As an example: a=1512, N=30 yields 0; but a=1512, N=300 yields NaN. – user46688 Oct 23 '15 at 22:25
  • You will need more than double precision, even if denormalised, to represent $e^{1512}$. – copper.hat Oct 23 '15 at 22:28
  • @copper.hat, that should be exp(-1512) which is 0. What goes to infinity is the multiplication of terms shown in your comment above in the original post. – user46688 Oct 23 '15 at 22:30
  • Well, the summation converges to $e^a$, so either you are computing (something close to) $e^a$ or $e^{-a}$, so you will run into precision issues. The correct answer (in infinite precision) is 1 in the limit, of course. – copper.hat Oct 23 '15 at 22:31
  • Is this for a numerical methods class? Are they trying to illustrate an issue or have you solve it? – copper.hat Oct 23 '15 at 22:32
  • Interesting problem. – copper.hat Oct 23 '15 at 22:51
  • 1
    I had more to say than would fit in a comment, and it seemed actually to be a possibly useful answer, so I wrote it as one. You might want to ask the real-life problem that inspired this question as another question. (If you do, link to the new question in a comment on this question so we'll be alerted.) – David K Oct 23 '15 at 23:15
  • Did I say something wrong? – copper.hat Oct 24 '15 at 21:28
  • Nope, I'm just re-evaluating based on the 2 new answers, as I haven't solved this yet, although your answer was the best as of yesterday. Wish I could give out +1 for your comments individually but I have to wait for more reputation points apparently. – user46688 Oct 24 '15 at 21:41
  • No problem, just curious :-). – copper.hat Oct 24 '15 at 21:42
  • You've been awesome, thanks. – user46688 Oct 24 '15 at 21:42
  • You might add to the question the size of $a$ you are interested in (as in $a=1512$). – copper.hat Oct 24 '15 at 21:43
  • Done, thanks... – user46688 Oct 24 '15 at 21:47
  • It occurred to me that since I ended up using the formula from this answer as the starting point for the second half of my answer, this answer deserves some credit too, so +1. – David K Oct 26 '15 at 19:52
1

For any real number $a$ and for any $N \geq 0$ the product $$e^{-a} \times \sum_{n=0}^N \frac{a^n}{n!},$$ that is, the ideal mathematical object that your computation is apparently supposed to represent, is a positive real number. For $N=0$ the value of the product is simply $e^{-a}$, but as you add terms (that is, for larger and larger $N$) the product approaches $1$, but is always at least a tiny bit smaller.

If you have used the expression exp(-a) in an attempt to represent $e^{-a}$, and the result of that expression is zero, it is not possible to get the correct result by computing the sum. That is, there is no number you could set the sum equal to, which, when multiplied by zero, would give you a result that is a positive number less than or equal to $1$.

So the difficulty occurs not because your computation of the sum overflows. The difficulty occurs because your computation of $e^{-a}$ underflows.

If you use a floating-point representation capable of representing $e^{-a}$ as a normalized non-zero value (that is, a value that may be very small, but is sufficiently different from zero that the computer can store it without losing too many digits of precision), then the sum should not overflow when you try to compute it.

Edit: copper.hat's technique is a relatively well-known way to add up terms when each term is a known multiple of the previous one. You may be able to use logarithms with an adaptation of that technique to make a reasonably accurate calculation of your desired product where $e^{-a}$ would underflow.

Using the method without logarithms, you might compute a sequence of values like this:

\begin{align} x_0 &= \frac{a}{N} \\ x_1 &= \frac{a}{N - 1}(1 + x_0) \\ x_2 &= \frac{a}{N - 2}(1 + x_1) \\ & \quad \vdots \\ x_i &= \frac{a}{N - i}(1 + x_{i-1}) \\ & \quad \vdots \\ x_{N-1} &= \frac{a}{1}(1 + x_{N-2}) \\ x_N &= 1 + x_{N-1} \\ \end{align}

The problem is that for very large $a$, the sequence grows very quickly near the end, and may overflow your chosen numeric type if you do this in a typical programming language without special "large number" support. So instead you can work with the sequence $y_0, \ldots, y_N$ defined by $y_i = \log(x_i)$. Assuming $\log$ is natural log, you can use the fact that $1 + x = x\left(1 + \frac 1x\right)$ to rewrite

\begin{align} \log(x_i) &= \log\left(\frac{a}{N - i}(1 + x_{i-1})\right) \\ &= \log(a) - \log(N - i) + \log(x_{i-1}) + \log\left(1 + \frac 1{x_{i-1}}\right) \end{align} as $$ y_i = \log(a) - \log(N - i) + y_{i-1} + \log(1 + \exp(-y_{i-1})). $$

Let $y_0 = \log(a/N)$, use the formula above to compute $y_i$ for $i = 1, 2, \ldots, N - 1$, and finally let $y_N = y_{N-1} + \log(1 + \exp(-y_{N-1})).$

At some point during this computatation, if $a$ and $N$ are both large enough, $\exp(-y_{i-1})$ will underflow, with the result that $\log(1 + \exp(-y_{i-1}))$ will evaluate to $0$ rather than to its correct (but very small) value. This is a bias toward underestimating the sum, but since $\log(1 + \exp(-y_{i-1})) < \exp(-y_{i-1})$, and since $\exp(-y_{i-1})$ will be tiny and (probably) very rapidly decreasing at that point, the bias may not be too much.

At the end, you can evaluate your product, using the fact that $$y_N = \log\left(\sum_{n=0}^N \frac{a^n}{n!}\right).$$

David K
  • 98,388
  • your analysis is correct. Sometimes, however, the summation overflows before $e^{−a}$ underflows. Thus, couldn't I do something like take the log of the expression, that is $log(e^{-a})=−a$ and add it to the log of the summation, then take the inverse log of the result? The problem then is how to take the log of the summation. I updated the original question to raise this possibility. – user46688 Oct 24 '15 at 21:57
  • The log of the summation is indeed where the difficulty would occur. I have added some thoughts to the end of the answer about how to deal with that calculation. – David K Oct 25 '15 at 00:45
  • By the way, I was counting denormalized values as a kind of "underflow"; that is, with IEEE-754 the exponent of $10$ in a positive number can be much more negative than it can be positive, but the most negative powers are achieved by discarding bits of precision, so really by the time the sum will have overflowed, $e^{-a}$ will have too little precision to be useful as a factor even if it has not underflowed all the way to zero. But it all leads to needing some other way to deal with very large numbers (e.g. logarithms) anyway. – David K Oct 25 '15 at 00:53
  • Thanks @DavidK, is $y_n$ as you show above (e.g. the log of the summed terms) equal to $y_N$? In other words, is that a typo (should $y_n$ really be written $y_N$)? If not, how to relate the two? – user46688 Oct 26 '15 at 18:55
  • Indeed, that was a typo, as you suspected. Oops! I have corrected it. – David K Oct 26 '15 at 19:00
  • Great, this really REALLY helped @DavidK, thank you so much. – user46688 Oct 26 '15 at 19:33
0

From its definition, you are guaranteed that, if $a > 0$, $e^{-a} \sum_{n=0}^N \frac{a^n}{n!} \le 1 $. (If $a < 0$ and $|a|$ is large, the sum will have catastrophic cancellation, so there is a problem from the start.)

Therefore, have your current expression be of the form $t_k =e^{-m}s_k $, where $s_k = \sum_{n=0}^k \frac{a^n}{n!} $.

Initially, $k=0$, $m = 0$, and $s_k = 1 $.

Each time $k$ is incremented, and the new $s_k$ and $t_k$ are computed, check if $t_k > 1 $. If so, increment $m$ until $t_k =e^{-m}s_k \le 1 $.

Once $k = N$, compare $m$ with $a$. If $m \ne a$, multiply by $e^{m-a}$. This will guarantee that you will never get overflows (though $\frac{a^n}{n!}$ might get small for $n >> a$, but that is another problem).

A problem with this is that the many multiplies by $\frac1{e}$ to get $e^{-m}$ might cause roundoff error to accumulate. This can be alleviated by techniques such as precomputed tables for $e^{-j}$.

marty cohen
  • 107,799