primes.py

import numpy as np
import pylab


def is_prime(n):
    status = True
    if n < 2:
        status = False
    else:
        for i in range(2,n):
            if n % i == 0:
                status = False
    return status


def get_first_n_primes(n):
    count = 0
    i = 1
    while True:
        i += 1
        if is_prime(i):
            count += 1
            yield i
            if count >= n:
                return


get_odds = lambda ps: (np.arange(ps.shape[0]) % 2) * -2 + 1


primes = np.array(list(get_first_n_primes(40)))
# Set 2 to 1
primes[0] = 1

sequence = 4. / (get_odds(primes) * primes)

pylab.close()
pylab.title("RE: Pi hiding in prime regularities")
pylab.plot(sequence.cumsum())
pylab.plot(np.pi * np.ones_like(primes))
pylab.legend(['Cumlitive sum of 4/x', 'PI'])
pylab.xticks(np.arange(primes.size), ["%+d" % prime for prime in get_odds(primes) * primes], rotation=-90)
pylab.xlabel("Primes")
pylab.show()