
This document describes the algorithms used in the mathematical constant calculators.

Note: The decision for when to print the next digit is based on numerical observations only, so there is always a chance that a digit will be printed before it has been correctly calculated. The Python versions here all print one million digits correctly, but I cannot prove that they will never fail.


Part 1.  pi = 3.14159...

The pi calculation uses Gosper's series.

Define the polynomials:

p(x) = 5x - 2
q(x) = x(2x - 1)
r(x) = 3(3x + 1)(3x + 2)

Gosper's series is:

pi = p(1) + p(2)q(1) / r(1) + p(3)q(1)q(2) / r(1)r(2) + p(4)q(1)q(2)q(3) / r(1)r(2)r(3) + ...

or, in factored form:

pi = p(1) + q(1) / r(1) * ( p(2) + q(2) / r(2) * ( p(3) + q(3) / r(3) * ( ...

Let the column vector [ u ; v ] represent the rational number u/v

The expression a + b/c * (u/v) can now be rewritten as a matrix-vector product:

[ b , ac  ;  0 , c ] * [ u ; v ]

In this form, the sum of the first n terms of Gosper's series is:

pi(n) = G(1) * G(2) * ... * G(n) * [ 0 ; 1 ]

where G(n) is defined as  [ q(n) , p(n)*r(n)  ;  0 , r(n) ]

The polynomials are calculated from repeated differences.
If the values of a polynomial at equally spaced intervals are listed, the differences between consecutive numbers are the values of a lower degree polynomial.
For example:

1  6  15  28  45    <-- q(x) for x = 1,2,3,4,5
5  9  13  17    <-- 1st differences, linear
4  4  4    <-- 2nd differences, constant

Once you have the first number in each line, the process can be reversed, calculating subsequent numbers with a series of additions.

For the extraction of decimal digits, define the matrices:

M = [ 10 , 0  ;  0 , 1 ]
S(d) = [ 1 , -d  ;  0 , 1 ]

If d is the integer part of u/v,  this can be subtracted with

S(d) * [ u ; v ]

Multiplying the result by 10, ready for the next digit, gives

M * S(d) * [ u ; v ]

Since all of the matrices form one big matrix product, with the digit extraction matrices multiplied on the left, and the series update matrices on the right, the two processes can be interleaved without affecting each other. Each extra term of Gosper's series gives one extra digit, so the code alternates the processes.

Putting it all together gives this Python code, which is the basis for the APGsembly:

#######################################
# product matrix [ a , b  ;  0 , c ], initially G(1)
a, b, c = 1, 180, 60

# polynomials and differences, initially p(1), q(1), r(1)
p, pd = 3, 5
q, qd, q2d = 1, 5, 4
r, rd, r2d = 60, 108, 54

flag = 1

# loop has no upper limit in APGsembly version
for i in range(1000):
   digit = b // c
   b %= c
   print(digit, end = '')

   if flag != 0:
      print('.', end = '')
      flag = 0

   a *= 10
   b *= 10

   p += pd
   q += qd; qd += q2d
   r += rd; rd += r2d

   b += a * p
   b *= r
   a *= q
   c *= r

print()
#######################################


Part 2.  e = 2.71828...

This calculation uses the continued fraction for e, which is faster than the usual infinite series:

e = 2 + 1 / (1 + 1 / (2 + 1 / (1 + 1 / (1 + 1 / (4 + 1 / (1 + 1 / (1 + 1 / (6 + ...

Continued fractions can also be converted to matrix form, with the expression a + 1 / (u/v) becoming

[ a , 1  ;  1 , 0 ] * [ u ; v ] = F(a) * [ u ; v ]

Using this, approximations for e are given by:

e(n) = F(2) * F(1) * F(2) * F(1) * F(1) * F(4) * F(1) * F(1) * ... * F(2n) * [ 1 ; 0 ]

The starting matrix is the product of the first six F matrices, with each update matrix including three more, giving:

e(n) = [ 87 , 19  ;  32 , 7 ] * FF(6) * FF(8) * ... * FF(2n) * [ 1 ; 0 ]

where FF(2n) = F(1) * F(1) * F(2n) = [ 4n + 1 , 2  ;  2n + 1 , 1 ]

Digit extraction works in the same way as for pi.
Each update matrix initially gives three extra digits, but after 12 updates the rate of convergence has improved enough to allow four digits per update.

Since the approximations alternate between too high and too low, correctness could be guaranteed by printing a digit only when successive approximations agree on its value. (The same applies to the next two algorithms.) I'll leave this modification to someone else.

In Python it looks like this:

#######################################
# product matrix [ a , b  ;  c , d ], initially all terms up to F(4)
a, b, c, d = 87, 19, 32, 7

# k is 2n + 1 in the formulas above
k = 7

flag = 1
updates = 0

# loop for 1000 digits; no limit in APGsembly
for i in range(253):
    b += 2 * a
    a = b * k - a
    d += 2 * c
    c = d * k - c
    updates += 1

    for j in range(3 if updates <= 12 else 4):
       digit = a // c
       a -= digit * c
       b -= digit * d

       a *= 10
       b *= 10

       print(digit, end = '')
       if flag == 1:
          print('.', end = '')
          flag = 0

    k += 2

print()
########################################


Part 3.  phi = 1.61803...

This uses another continued fraction, but not the one for phi, which is too slow. Instead it uses

1 + sqrt(5) = 3 + 1 / (4 + 1 / (4 + 1 / (4 + ...

Converting to matrix form, and doubling the bottom row of the starting matrix, gives approximations to phi:

phi(n) = [ 3 , 1  ;  2 , 0 ] * F(4) ^ n * [ 1 ; 0 ]

It makes the APGsembly code simpler if alternate matrix products are stored in mirror image form, that is

[ b , a  ;  d , c ] instead of [ a , b  ;  c , d ]

This makes no difference to the digit extraction, as a/c and b/d have the same integer part.
Each F(4) matrix gives one extra digit.

Python code:

#######################################
a, b, c, d = 3, 1, 2, 0
flag = 1
direction = 0

for i in range(1000):
   if direction == 0:
      b += 4 * a
      d += 4 * c
   else:
      a += 4 * b
      c += 4 * d
   direction = 1 - direction

   digit = a // c
   a -= digit * c
   b -= digit * d

   print(digit, end = '')
   if flag != 0:
      print('.', end = '')
      flag = 0

   a *= 10
   b *= 10

print()
#######################################


Part 4.  sqrt(2) = 1.41421...

Yet another continued fraction, this time:

sqrt(2) = 1 + 1 / (2 + 1 / (2 + 1 / (2 + ...

The first three terms are in the starting matrix, giving:

sqrt2(n) = [ 7 , 3  ;  5 , 2 ] * F(2) ^ n * [ 1 ; 0 ]

As with phi, the product matrix alternates with its mirror image.
Four F(2) matrices give three extra digits.

Python code:

#######################################
a, b, c, d = 7, 3, 5, 2
flag = 1
skip = 2
direction = 0

for i in range(1333):
   if direction == 0:
      b += 2 * a
      d += 2 * c
   else:
      a += 2 * b
      c += 2 * d
   direction = 1 - direction

   if skip == 0:
      skip = 3
   else:
      skip -= 1
      digit = a // c
      a -= digit * c
      b -= digit * d

      print(digit, end = '')
      if flag != 0:
         print('.', end = '')
         flag = 0

      a *= 10
      b *= 10

print()
#######################################
