Computing irrational numbers

I am genuinely curious, how do people compute decimal digits of irrational numbers in general, and $\pi$ or nth roots of integers in particular? How do they reach arbitrary accuracy?

Solution 1:

I applaud your genuine curiosity. There is no general method. To compute the decimal digits for any particular irrational number you start with a definition of the number.

For $\pi$ there are lots of methods known. @GregoryGrant 's answer points to one of the oldest and most famous. Archimedes started with inscribed and circumscribed regular hexagons, then doubled the number of sides several times to increase the accuracy of his estimate. (See http://en.wikipedia.org/wiki/Pi#Polygon_approximation_era).

Searching for "computing digits of pi" will lead you to information about modern methods.

For roots of polynomials (like nth roots of integers) Newton's method is efficient.

For many irrationals, continued fractions are a good bet, as @Arteom notes.

For some particular irrationals, clever tricks (based on deep theory) give you lots of digits quickly. For example, the fractions 3/2, 7/5, 17/12, 41/29, ... are better and better approximations to $\sqrt{2}$. You should be able to guess the pattern.

Solution 2:

Here's how Archimedes did it 2,500 years ago. $\pi$ is the circumference of a circle with diameter one. What Archimedes did was to inscribe a polygon with $n$ sides inside the circle. Like this for $n=5$.

He then calculated the circumference of the polygon. This gives an approximation to $\pi$ for each $n$ and as $n$ grows it gets closer and closer to $\pi$.

Solution 3:

You generally need to use arbitrary precision arithmetic to compute large numbers of digits of typical irrational numbers. The exception is oddball things like the Champernowne constant 0.12345678910111213141516… :)

There are various arbitrary precision arithmetic packages available. Internally, they use big arrays or lists to hold groups of the digits of the numbers they manipulate. The computer can do simple arithmetic on numbers represented this way using algorithms similar to those taught in grade school which allow you to do arithmetic on multi-digit numbers given the ability to do arithmetic on single digit numbers.

However, for very large numbers the naive multiplication & division algorithms are rather inefficient, so more subtle techniques are used. Multiplication of huge numbers can be performed efficiently using the Fast Fourier transform. Division can be performed by multiplying by the reciprocal of the divisor; reciprocals can be calculated using Newton's method.

Once you have such algorithms at your disposal it's fairly straightforward to calculate $n$th roots using Newton's method, as Lucian mentioned in the comments.

However, it's not necessary to have a full arbitrary precision arithmetic package to calculate large numbers of digits of some of the well-known irrational numbers. It's trivial in most programming languages to add or subtract large numbers that are stored as an array of digit blocks. It's also quite easy to multiply such numbers by a normal-sized number; division by a normal-sized number is similarly easy to implement. And that's all you need to calculate large numbers of terms of many Taylor series to huge precision.

One somewhat famous example of this technique is this gem written in C by the late Dik T. Winter of CWI, which computes $\pi$ to 800 decimal digits.

int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;
     for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,
     f[b]=d%--g,d/=g--,--b;d*=b);}

Unfortunately, that cryptic little program cannot be extended to yield more (correct) digits, due to overflow problems.

But here's a program I wrote a few years ago that can generate large numbers of digits of $e$; it can calculate 10,000 digits in a couple of seconds. The original version of this program was written in C, this version is in Python 2; as you might guess the C version is significantly faster :). Hopefully, the comments explain the algorithm, even if you aren't familiar with Python. The core idea is to treat the Taylor series of $e$ as a number written in factorial base and to simply convert that number into decimal notation.

#! /usr/bin/env python

"""
    ECalc - Big e calculator

    Generate digits of e, using factorial base notation
    Inspired by Dik T. Winter

    Created by PM 2Ring 2003.12.26
    Converted from C to Python 2007.12.10
"""

import sys, math

# Global constants
N = 50                            # Default number of digits
W = 5                             # Digits per block
A = 10 ** W
LR2P = .5 * math.log(2.*math.pi)  # Log Root 2 Pi
L10 = math.log(10.)               # Log 10

def invlfact(n):
    """ Inverse Stirling's approx for log n factorial, using Newton's method """
    x = y = n * L10 - LR2P
    for i in xrange(3):
        x = (y + x) / math.log(x)
    return int(round(x))

def bigE(n):
    """ Generate n digits of e using factorial base notation """
    kmax = invlfact(n) + 1
    print "Calculating e to %d places, using %d cells" % (n, kmax)

    # Table for factorial base digits, initialized with series for e-2 = 1/2! + 1/3! + ...
    # We ignore first two entries
    d = [1] * kmax
    print "e ~= 2."

    j = 1
    while n>0:
        # Get the next W digits by firstly multiplying d by A,
        # propagating carries back down the array,
        # then printing the integer part & removing it.
        kmax = invlfact(n)    # Number of cells needed for this loop
        c = 0                 # Clear carry
        for k in xrange(kmax, 1, -1):
            d[k] = d[k] * A + c
            c = d[k] // k
            d[k] -=  c*k

        # Print block & field separator. May need modifying if you make W large
        jj = (j%10 and 1  or  j%50 and 2  or  j%200 and 3  or  4) - 1
        print "%0*d%s" % (W, c, "\n" * jj),

        j += 1; n -= W

def main():
    # Get number of digits
    n = len(sys.argv) > 1 and int(sys.argv[1]) or N
    bigE(n)

if __name__ == '__main__':
    main()  

According to MathWorld

Around 1966, MIT hacker Eric Jensen wrote a very concise program (requiring less than a page of assembly language) that computed e by converting from factorial base to decimal.

I expect that Eric's algorithm was substantially the same as the one in my code, which I derived from an obfuscated C program written by the late Dik Winter.

Solution 4:

For some irrational numbers, like $\pi$, there are convenient infinite series that converge to them. So for example $$\sum_{n=1}^\infty\frac{1}{n^2}=\frac{\pi^2}{6}$$ By adding up more and more terms of this series you get closer and closer to $\pi^2/6$. You can take an estimate for some value of $n$, multiply by $6$ and then take the square root. That will give you an approximation of $\pi$. This illustrates the idea, but in practice nobody uses this series because there are series that converge to $\pi$ much faster.