Projecteuler.net Problem 381

p381.py

#!/usr/bin/env python 
"""
  Projecteuler 
  Problem 381 
  
  Wilson's theorem
"""

import fastlib 
import math
import multiprocessing
import cPickle as pickle
n = 10**8
#n = 100

def f1(x):
    return x if fastlib.is_prime(x) else 0
try:
    print 'Loading prime numbers...'
    ps = pickle.load(open('./prime.dat'))
except:
    print 'Generating prime numbers...'
    p = multiprocessing.Pool(4)
    ps = filter(lambda x: x!=0, p.map(f1, xrange(5,n+1)))
    print 'Len(ps) =', len(ps)
    if raw_input('Are you going to write to prime.dat (Y/N)?') == 'Y':
        pickle.dump(ps, open('./prime.dat','w'))
print 'Calculating S(n)...'


def invmod(a,m):
    sm = m
    p = [0, 1]
    q = []
    i = 0
    while True:
        r = m % a
        if i >= 2:
            p.append((p[-2] - p[-1]*q[-2]) % sm) 
        q.append(m // a)
        m=a
        a=r
        if r == 0:
            if len(q) < 2:
                p.append(1)
            else:
                p.append((p[-2] - p[-1]*q[-2]) % sm) 
            return p[-1]
        i+=1
def r1(r):
    return r-1
def r2(r):
    return 1
def r3(r):
    return (r-1)/2
def r4(r):
    return (invmod(r-3,r) * r3(r)) % r
def r5(r):
    return (invmod(r-2,r) * invmod(r-3,r) * invmod(r-4,r)) % r

def S(p):
    r = [r1(p), r2(p), r3(p), r4(p), r5(p)]
    return sum(r) % p

p = multiprocessing.Pool(4)
s = p.map(S, ps)
print 'Reducing...'
print sum(s)