Directly calculating Graham’s Number is impossible, to the point that even writing out approximations of what the calculation should look like exceed the capacity of the universe in all but the simplest forms. But we can have some fun in calculating what some of the least significant digits look like by observing the convergent properties of “Power Towers”.

Using a simple formula to calculate the rightmost digits of Graham’s Number is straight forward:

def calc_graham(d):
        x = 3
        for n in range(0,d):
                x = (3 ** x) % (10 ** d)
        return x

This is, of course, horrifically slow. Depending upon the number of digits d we want to find, we’ll be calculating 3 to the power of a number of length d digits d times. Trying to calculate anything more than 6 or 7 digits will bring most computers to a crawl, much less trying to compute 100 or more digits.

Modular Exponentiation to the rescue. Modular exponentiation is a technique used throughout computer science (especially cryptography) to simplify the computation of equations of the type:

C = b^e mod m

Our calculation of the rightmost digits in the Graham Number fit this exactly, all we need to do is whip up a modpow method for Python integers:

def modpow(base, exponent, modulus):
        result = 1
        while exponent > 0:
                if exponent & 1 == 1:
                        result = (result * base) % modulus
                exponent = exponent >> 1
                base = (base * base) % modulus

This method is based upon an algorithm in Bruce Schneier’s Applied Cryptography. With this algorithm in place, the entire program to calculate the rightmost digits of Graham’s Number becomes:

#!/usr/bin/python
import sys
import time

def calc_graham(d):
        x = 3
        td = 10**d
        print "Iterating..."
        for n in range(0,d):
                print " step %d" % (n + 1)
                x = modpow(3, x, td)
        return x

def modpow(base, exponent, modulus):
        result = 1
        while exponent > 0:
                if exponent & 1 == 1:
                        result = (result * base) % modulus
                exponent = exponent >> 1
                base = (base * base) % modulus
        return result

if __name__ == '__main__':
        digits = int(sys.argv[1])
        partitionSize = 50
        if len(sys.argv) > 2:
                partitionSize = int(sys.argv[2])

        st = time.time()
        graham = calc_graham(digits)
        ed = time.time()

        print "calculation took %f seconds" % (ed - st)
        print "right most %d digits of G are:\n" % (digits)

        grahamString = "%d" % (graham)
        while len(grahamString) > 0:
                print grahamString[:partitionSize]
                grahamString = grahamString[partitionSize:]

Calculating the least significant digits of Graham’s Number in this way is still limited in efficiency by the number of digits being computed, but is orders of magnitude faster than the straight forward approach. And the last 500 digits of G?

24259506950647383956574791365193517983345353625214300354012602677162267216041981
06522631693551887803881448314065252616878509555264605107117200099709291249544378
88749606288291172506300130362293491608025459461494578871427832350829242102091825
89675356043086993801689249889268099510169055919951195027887178308370183402364745
48882222161573228010132974509273445945043433009010969280253527518332898844615089
40424826501819385156253579639961899396790549663800322234872396701848518643905910
4575627262464195387