# Concision and Concinnity

## 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 randintfrom operator import mulfrom collections import dequefrom decimal import Decimal, getcontext# Allow enough precision to get an answer getcontext().prec = 1000def 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 += 1def 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 Truedef 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            breakif __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`