## 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 "Found! First ten-digit"

print "prime in consecutive digits"

print "of e is %s!" % chunk

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