[Edu-sig] More fun with Euler's Number

[Kirby Urner]

>
>The approximation re terms needed, given request for d digits, could
>probably use some fine tuning no doubt.

Here's another stab at generating e to more-than-we-usually-get
decimals.

There's this algebraic transformation, which I haven't quite
figured out yet, that lets you rewrite 1 + 1/1! + 1/2! + 1/3!
+ 1/4! + 1/5! as 1+(1+(1+(1+(1+1/5)/4)/3)/2)/1.  This translates
into a very simple algorithm for e (see below), used by
Jean-Michel Ferrard in a program for the TI (calculator).
He manages to get 600 significant digits out of the thing
(a TI-92 or 89).

But this algorithm requires knowing in advance how many terms
you need in order to produce n significant digits.  Jean-Michel
knows 294 is the magic number for 600 digits, but the best
approximation I was able to find, for the number of terms,
was here: http://numbers.computation.free.fr/Constants/E/e.html
where is written:  "To have d decimal digits of e, the first
K ~ d log(10)/log(d) terms of the series are sufficient"
(Section 7).

But that's approximate, and it doesn't work for a lot of the
values of d that I try.  My earlier approach was to boost it
higher -- but how much higher would I need to go?  I went
too far.

If we're really not sure in advance how many terms to
compute, how do we stop from running through too many?
I thought a good way to apply the brakes would be to
check the least significant bits of the long integer
being completed, and to say "if these bits don't change
through 10 iterations, then we're done".  I also go
to a few digits beyond what the user requests, for
good measure, but only display the number of decimals
(= terms) actually asked for.

Here's a program implementing this strategy, along with a
different form of the algorithm derived in Section 4 of the
above web resource.

def gete(terms):
     e = 10 ** (terms + 5)
     sig = u = n = 1
     cnt = 0
     while cnt < 10:
        e = (e * (u+1))//u
        u = (n+1)*(u+1)
        n += 1
        leastsig = e & 1111

        if leastsig == sig:  cnt += 1
        else:
           sig = leastsig
           cnt = 0

     se = str(e)
     return se[0]+"."+se[1:terms+1]

I asked for the first 6000 digits of e, and was able to compare them
with an online resource to ascertain my result indeed ended with the
right digits:  ...6110250517100   It goes through about 2090
iterations to reach the terminal condition, whereas the log formula
would suggest 1588.

Takes about 50 seconds to complete on my 1.2 Ghz AMD w/ 640 MB RAM
(ymmv).

However, now that I know I should use about 2090 iterations, if
I use the algorithm in Jean-Michel Ferrard's paper (it was
emailed to me -- not on the web AFAIK), then I get 6000 digits
in a blindingly short time, under 1 second!  Sheesh, that's a
*huge* difference.  Here's that algorithm:

   def gete2():
       "Return the value of e with 6000 digits"
       u = x = 10**6000
       for k in range(2090,0,-1):
           u = x + u//k
       return str(u)

Clearly, knowing the number of terms in advance, and using
the above (which implements the algebraic transformation I
mention up top), is the superior method, by a long shot.
I need to get a better grip on this "knowing how many terms
in advance" issue.

Of course none of these methods are what the professionals
are using to compute e to millions or billions of digits.
That's a whole other ballgame, using 'binary splitting' and
the like (involves Fourier Transforms -- 'Who is Fourier?').
Still, this game is fun enough, if not ground-breaking
(and it's educational).

Kirby