why 'lambda' and 'reduce'?

[Manuel Garcia]

I had a few people email me, asking how many of the 5000 digits of Pi
are accurate, and being surprised that all 5000 are correct.  I can
explain how it works as best I can (I am not a professional
mathematician).

The limits of how many digits it can compute are only the limits of
Python's long integer implementation, and the limits of the computer's
memory.  1,000,000 digits, probably, if you don't mind the wait.
1,000,000,000 digits, don't think so.  It is hardly the most efficient
way to calculate Pi, but it also is not grotesquely inefficient.  In
fact, it runs surprisingly fast in Python 2.3 for 5000 digits.

The code contains absolutely no 'obfuscation', except the use of '+'
for string concatenation.  All the 'lambda's and 'reduce' and whatnot
actually have to be there for the efficiency of the algorithm; to make
the code run as fast as a straightforward implementation.  The whole
of the obfuscation is that Python's built-in functions and long
integers make coding the algorithm unusually 'clean' (if that mess can
be called 'clean').

The code uses a 'quadratic iteration', referring to the speed with
which correct digits are produced.  With each iteration the number of
correct digits is doubled.  Modern computers of massive numbers of
digits of Pi (many billions) use such iterations.  (There are newer
techniques for computing just a specific digit of Pi, without having
to compute the digits before, but if you wanted the unbroken list of
digits, you would use an iteration algorithm.)

A helpful webpage is:
http://numbers.computation.free.fr/Constants/Pi/iterativePi.html
(The iteration I used is #1.2.1, found in 1976 independently by Eugene
Salamin and Richard Brent)

An article in Scientific American (I haven't read it, but a couple of
different websites said it was fun and easy to understand):
Borwein, Jonathan M., and Peter B. Borwein. RAMANUJAN AND PI; February
1988, page 66

If you wish to download a massive number of digits of Pi, to confirm
all this:
http://www.geocities.com/thestarman3/math/pi/picalcs.htm

Geometrical or calculus based algorithms for computing Pi typically
add a fixed number of correct digits with each 'step' (either a single
compass and ruler construction in an infinite construction or a single
term in an infinite sum or product).  These newer techniques are based
on pretty 'heavy-duty' number theory, increase the number of correct
digits times 2, times 3, times 4, times 5, etc.  (I think I saw the
construction of an iteration that increased times 16)  The Japanese
group who hold the record on computing digits of Pi (51 billion
digits) used a 'times 4' (quartic) iteration, optimized for binary
computers.  Their iteration actually computes the reciprocal of Pi,
making it hard to inspect the digits while you are creating them.

I picked the simplest quadratic iteration, because for 5000 digits it
hardly makes any difference.  The '13' in the code is the number of
iterations needed.  For 10000 digits, you would use 14 iterations,
etc.

To maintain accuracy to 5000 digits, I could not use Python floating
point.  I had no choice but to use long integers, much like Tim
Peters' 'FixedPoint'.  The integers are huge (around 10**5000 to
10**10000 in size), but Python has no limits on the size of long
integers.  (Python 2.3 uses the Karatsuba algorithm to multiply long
ints, much faster than 2.2.  I don't think division was speeded up).
To allow use of longs instead of decimals, I had to scale up all
calculations by 10**5010 (I needed 10 extra decimal places to
guarantee the accuracy of the 5000 digits, because I haven't
eliminated rounding errors, but I only need 10 (maybe less, I didn't
really try) because the iteration is pretty numerically stable)

What looked like 'geometrical mean' was actually my having to
implement a square root function with Newton's method, because
Python's built in pow() converts to floating point with fractional
exponents, and 10**5000 is too large for the built in floating point.
The number '15' is the number of iterations Newton's method needs (it
is also quadratic, I threw a couple of extra iterations in there to
keep the Pi iteration happy).

If you go through the expression step by step, you will see no
'lambda' or 'reduce' can be removed, unless you sacrifice speed or
accuracy.  So it is true it is hardly obfuscated, because my hands
were tied almost every step of the way.  'lambda' took the place of
variable assignment, and my bizarre 'reduce' expressions were used for
looping.

Anyway, here is the code.  I measured x*x as running faster than x**2
or pow(x,2) using 2.3's 'timeit' module, but it really doesn't make
any difference, as the running time is swamped by the long integer
division and multiplication.

print (lambda p:p[0]+'.'+p[1:])(
    str((lambda(x,y,t,a):(2L*x*x)//a)(
        (lambda F:(lambda S:reduce(
            lambda(x,y,t,a),_:((x+y)//2L,
                S((x*y)//F),2L*t,
                (a-(t*(((x+y)//2L)**2-
                       (S((x*y)//F))**2))//F)),
            [0]*13,(F,(F*F)//S(2L*F),2L,F//2L)))(
                lambda n:reduce(lambda x,_:(
                    x-x//2L+(n*F)//(2L*x)),
                [0]*15,
                n//2L)))(10L**(5010))))[:5000])

Manuel