Summation on HP 42S

[Albert Chan]

Hi, Frido Bohn

For reference, 1000! ~ 4.023872600770937735437024339230039857194E+2567

To avoid overflow, sum log10 of the integers.
A straight sum will loses a few significant digits.

>>> from math import *
>>> n = 1001
>>> frac = lambda x: x - int(x)
>>> sum(log10(i) for i in xrange(2, n))
2567.6046442221304
>>> 10 ** frac(_)
4.0238726007486756

If error terms is collected, result is 1000X more accurate

>>> sum = error = 0
>>> for i in xrange(2, n):
...         v = log10(i)
...         old = sum
...         sum += v
...         error += v - (sum - old) ### collect what is not added to sum
...
>>> sum, error
(2567.6046442221304, 2.4005242238445135e-12)
>>> 10 ** (frac(sum) + error)
4.0238726007709165

Update: we could collect powers of 10, without log10()

lua> p, e = 1, 0
lua> for x = 1, 500, 2 do
        p = x * (x+1) * (1000-x) * (1001-x) * p
        if p > 1e200 then e=e+300; p=p*1e-300 end
        end
lua>
lua> p, '* 10^' .. e
4.0238726007709384e+167 * 10^2400