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

0 comments

Post a Comment