Update!

Since my last update in 2010, I've been working tirelessly day and night — hardly sleeping — to achieve my lifetime goal of calculating 1,000,000,000,000 decimal digits of pi. My latest strategy involves waiting for desktop scientific calculators to mature just a bit more, then I'll type 'pi' into the calculator and export the results.

It's gonna be epic.

100,000

I've passed the 100,000 digit milestone.  My program appears below.  I know that there are some serious optimizations that should take place.  The slowest part of the program seems to be the Euclid GCD algorithm, though other factors will become limiting as the terms grow longer.  Next stop is 1,000,000 digits.


import sys


#
# ASSORTED CONSTANTS AND FUNCTIONS
#


EXP_BASE = (-640320)**3


saved_k = 0
saved_n1 = 1 # (6k)!
saved_n2 = 13591409 # 13591409 + 545140134*k
saved_d = 1 # (3k)! * (k!)**3 * EXP_BASE**k


def factorial(n, final = 0):
  f = 1
  while n > final:
    f *= n
    n -= 1
  return f


def add(a, b):
  'a and b are fraction tuplets'
  n = a[0]*b[1] + b[0]*a[1]
  d = a[1]*b[1]
  return reduce(n, d)


def intsqrt(x, bot = 0, top = None, print_updates = False):
  if x == 1: return 1 # special case
  if not top: top = x
  i = 1
  expected = 3*len(str(x)) if print_updates else None
  while True: # return when done
    if print_updates:
      if not i % 1000:
        sys.stderr.write('Performing iteration %s of about %s.\n' % (i, expected))
      i += 1
    a = (bot + top)/2
    b = a*a
    if b > x:
      top = a
      continue
    if b < x:
      if a == bot: return a
      bot = a
      continue
    return a


def gcd(x, y): # "The laws of nature are but the mathematical thoughts of God."
  while x:
    y, x = x, y % x
  return y


def reduce(n, d):
  g = gcd(n, d)
  return (n/g, d/g)


def calculate_term(k):
  global saved_k, saved_n1, saved_n2, saved_d


  # calculate the saved parts
  if k == saved_k + 1:
    # we can reuse the saved parts
    saved_n1 *= factorial(6*k, 6*(k - 1))
    saved_n2 += 545140134
    saved_d *= factorial(3*k, 3*(k - 1)) * k**3 * EXP_BASE
  else:
    # calculate from scratch
    saved_n1 = factorial(6*k, 3*k) # we kind of cheat here...
    saved_n2 = 13591409 + 545140134*k
    saved_d = factorial(k)**3 * EXP_BASE**k
  saved_k = k


  # reduce and return
  # (note that we could reduce returned value further, but we
  # can't reduce saved_n2, so we let add do any reduction)
  (saved_n1, saved_d) = reduce(saved_n1, saved_d)
  return (saved_n1*saved_n2, saved_d)


def sum_terms(start_term, end_term, print_updates = False):
  k = start_term
  s = (0, 1)
  while k <= end_term:
    if print_updates:
      sys.stderr.write('Performing iteration %s of %s.\n' % (k - start_term + 1, end_term - start_term + 1))
    s = add(s, calculate_term(k))
    k += 1
  return (s[0], s[1])


def C(digits, print_updates = False):
  'Digits is number of digits needed.  Divide by 10**(digits - 8).'
  if digits < 8: raise Exception('digits must be at least eight')
  n = digits + 7 # to get accurate digits in the least significant positions
  return 426880*intsqrt(10005 * 100**(n - 8), print_updates = print_updates)/10000000


def pi(s, c):
  return c*s[1]/s[0]


def usage():
  print 'Usage: pi.py [-t ] [-c ] [-s ]'
  print '             [--pi ] [--test]'
  print '             [-v ]'
  print ' -t calculates the sum of terms from to '
  print ' -c calculates the constant 426880 * intsqrt(10005) * 10**'
  print ' -s calculates the sum of two terms (loads from files)'
  print ' --pi calculates pi (loads from files)'
  print ' --test run program self-test'
  print ' -v verifies pi'
  print
  exit()


#
# KNOWN PI STRING FOR TESTING ONLY
#


known_pi_string = '314' # shortened for blog


#
# MAIN PROGRAM
#


if len(sys.argv) < 2: usage()


if sys.argv[1] == '-t':
  if len(sys.argv) != 4: usage()
  print sum_terms(int(sys.argv[2]), int(sys.argv[3]), print_updates = True)


elif sys.argv[1] == '-c':
  if len(sys.argv) != 3: usage()
  print C(int(sys.argv[2]), print_updates = True)


elif sys.argv[1] == '-s':
  if len(sys.argv) != 4: usage()
  first = open(sys.argv[2]).read()
  second = open(sys.argv[3]).read()
  #TODO security check inputs here
  print add(eval(first), eval(second))


elif sys.argv[1] == '--pi':
  if len(sys.argv) != 4: usage()
  s = open(sys.argv[2]).read()
  c = open(sys.argv[3]).read()
  #TODO security check s here
  print pi(eval(s), int(c))


elif sys.argv[1] == '--test':
  pi_string = str(pi(sum_terms(0, 100), C(1000)))
  i = 0
  m = min(len(known_pi_string), len(pi_string))
  while i < m and pi_string[i] == known_pi_string[i]:
    i += 1
  print '%s out of %s digits passed' % (i, m)


elif sys.argv[1] == '-v':
  if len(sys.argv) != 3: usage()
  pi = open(sys.argv[2]).read()
  if len(pi) > len(known_pi_string):
    sys.stderr.write('Problem: only %s known digits\n' % len(known_pi_string))
    exit()
  i = 0
  while i < len(pi) and pi[i] == known_pi_string[i]: i += 1
  print '%s out of %s digits passed' % (i, len(pi))


else:
  usage()


sys.stderr.write('Completed %s\n' % ' '.join(sys.argv[1:]))

One Trillion Digits

Pi was first calculated to over one million decimal digits in 1973, the year before I was born.  The current world record is almost three trillion digits, and I have little doubt that the quadrillion digit mark will be achieved this decade.  Why, then, would one set a personal goal of calculating one trillion digits of pi?  Because it's there, smirking down at us, daring that we master it.  I'm not the first climber to attempt to scale this mountain, and I certainly will not be the last.