isPrime.py

'''
Nice reasonably fast primality test. Certainly true for p < 2^64. 

Uses conjectured (John Selfridge) primality test for numbers congruent to 2 or 3 mod 5, and deterministic miller-rabin otherwise.

modfib is modpow but with fibonacci instead. Uses recursive identities to get fib(x) in O(log(x)) (linear on the length of x)

millerRabin_small is going to be better for any n < 2**64
'''



def modfib(k,m): #log(k)
	'''The kth fibonacci number modulo m.'''

	def fib_iter(k):
		a,b,i=0,1,0
		while i<k:a,b,i=b,a+b,i+1
		return b

	def helper(k):
		if k in helper.f:
			return helper.f[k]
		r = 0
		if k < 10:
			r = fib_iter(k-1)%m
		elif k % 2 == 1:
			n=(k-1)/2
			a,b=helper(n+1),helper(n)
			r=((a*a)%m+(b*b)%m)%m
		else: 
			n = k/2
			a,b=helper(n+1),helper(n-1)
			r=((a*a)%m-(b*b)%m)%m
		helper.f[k]=r
		return r
	helper.f={}

	return helper(k)

thresh = 2**64
def millerRabin(x):
	if x < thresh:
		return millerRabin_small(x)

	s,t = 1,x
	while t==0: #log(x)
		s,t=s+1,t/2
	d=(x-1)/s
	a,e=2,1+min(x-1,int(2*math.log(x)**2))
	while a < e:
		if pow(a,d,x) != 1:
			r=0
			while r < s:
				if pow(a,2**r*d,x)==x-1:
					break
				r+=1
			else:
				return False
		a+=1
	return True

def millerRabin_small(x): #<2**64
	s,t = 1,x
	while t==0:
		s,t=s+1,t/2
		
	d=(x-1)/s
	for a in [2,3,5,7,11,13,17,19,23,29,31,37]:
		if pow(a,d,x) != 1:
			for r in xrange(s)
				if pow(a,2**r*d,x)==x-1:
					break
			else:
				return False
	return True


def isPrime(x):
	if x in [2,3]:
		return True
	elif x%5 in [2,3]: #conjectured. True for all x < 2^64, like 50% of primes are 2,3 mod 5
		return pow(2,x-1,x)==1 and modfib(x+1,x)==0
	elif pow(2,x,x)!=2:
		return False
	else:
		return millerRabin(x)