How Joseph Method Stole My Day

My friend Method shared this with me via Google Reader this morning, accompanied by the provocative message, "Here you go, Ian." The problem, in short, was to discover the first ten-digit prime in consecutive digits of e. And a day of planned work-catching-up dissolved before it began.

I had a solution in twenty minutes or so, using third-party libraries pycrypto (for isPrime) and mpmath (for e), but that made the whole enterprise trivial and didn't quite seem cricket. Also, it was slightly inefficient, because, being limited by the code of others, I had to calculate e to n digits first, and then check that for primes; I figured if I could calculate digits of e lazily I'd get an extremely minor but nonetheless pride-inducing speed increase.

To Wikipedia! where I learned scads about primality tests (luckily I didn't read the page on the constant e until later, or I would have learned that this was a common puzzle that pretty much everybody tries to solve. Also I would have learned the answer). I remembered how to calculate e using the Taylor series 1 + 1/1! + 1/2! + 1/3!..., so that wasn't a problem—until I ran into precision woes.

That led me down a dark path whereon I attempted to program arbitrary-precision arithmetic from scratch, and learned all about mantissas and bitwise operators over the course of five hours. I got pretty far with that before I realized that I was spending my entire day trying to comprehend what I imagine amounts to at least a month of college-level CS. I went to school for other things, like French poetry.

Also at this point I remembered the existence of Python's own decimal floating-point arithmetic package, decimal, which had apparently dropped completely out of my brain exactly when I needed it. Thankfully, it resurfaced before I had to give up on this psychotic enterprise altogether.

After that it was easy. I post the code below not because I believe it is in any way notable—I'm sure it's old hat to those who went to school for this—but because I just want someone to know that I actually did something today. Thanks, Method.

from random import randint
from operator import mul
from collections import deque
from decimal import Decimal, getcontext

# Allow enough precision to get an answer 
getcontext().prec = 1000

def factorial(x): 
    if x==1: return x
    return reduce(mul, xrange(2, x+1))

def taylor_series():
    """
    e == 1/1! + 1/2! + 1/3! ...
    """
    s = Decimal(1)
    _e = Decimal(1)
    while True:
        _e += (Decimal(1)/factorial(s))
        yield _e 
        s += 1

def digits_of_e():
    """
    Yield digits of e lazily by noticing when
    they stop changing as the Taylor series 
    converges.
    """
    idx = 2 # skip '2.'
    last = None
    for iteration in taylor_series():
        i = str(iteration)
        try:
            cur = i[idx]
        except IndexError:
            continue
        if cur==last:
            idx += 1
            yield cur
        last = i[idx]

def ten_digit_chunks():
    """
    Collect digits of e as they are discovered
    and yield them in ten-digit chunks.
    """
    l = deque()
    for digit in digits_of_e():
        l.append(digit)
        # Keep the stack at 10 by popping
        # from the left
        while len(l)>10:
            l.popleft()
        if len(l)==10: 
            yield ''.join(l)

def is_prime(n, k=20):
    """
    Simple implementation of the Miller-Rabin 
    primality test. Thanks, Wikipedia!

    n is the integer to be tested; k is the 
    precision (number of random possible 
    witnesses to test against)
    """
    # Find s and d such that (2**s)*d==n
    d, s = n-1, 0
    while not d & 1: # d & 1 == d % 2, only faster
        d >>= 1 # d >> 1 == d / 2, only faster
        s += 1
    # Test for k random numbers
    for i in xrange(k):
        _c = False
        a = randint(2, n-2)
        x = pow(a, d, n) # a**d % n
        if x==1 or x==n-1:
            continue
        for r in xrange(1, s):
            x = pow(x, 2, n) # x**2 % n
            if x==1:
                return False
            elif x==n-1:
                _c = True
                break
        if _c: continue
        return False
    return True

def main():
    for chunk in ten_digit_chunks():
        print "Trying %s..." % chunk
        if is_prime(long(chunk)):
            print "="*20
            print
            print "Found! First ten-digit"
            print "prime in consecutive digits"
            print "of e is %s!" % chunk
            print
            print "="*20
            break

if __name__=="__main__":
    main()

UPDATE: In order to get through this ridiculous set of problems to the unsatisfying and insulting end, you can use this implementation of the Brent-Salamin algorithm to find primes in digits of pi, replacing taylor_series in the digits_of_e function in the code above:

def brent_salamin():
    """
    A series whose sum converges on pi
    incredibly quickly.
    """
    a = Decimal("1.")
    b = Decimal("1.")/(Decimal("2.").sqrt())
    t = Decimal("1.")/Decimal("4.")
    p = Decimal("1.")
    while True:
        a1 = Decimal(a + b)/Decimal("2.")
        b1 = Decimal(a*b).sqrt()
        t1 = t - p*((a - a1)**2)
        p1 = 2*p
        yield ((a1+b1)**2)/(t1*Decimal("4."))
        a, b, t, p = a1, b1, t1, p1