numbertheory.py

import time, math, sys, random
	
def gcd(x0, y0):
	if x0 == y0: 
		return -1

	while y0:
		x0, y0 = y0, x0 % y0
	return x0
	
def gcde(r0p, r1p):	
	if r0p < 1 or r1p < 1 :
		#return [math.max(r0p, r1p), 1, 1]
		return -1
	# Stop param mutation
	r0, r1 = r0p, r1p
	x0, y0, x1, y1 = 1, 0, 0, 1
	
	while 1:
		q = r0 / r1
		r = r0 % r1
		
		if r < 1:
			break
		
		x0, x1 = x1, x0 - q * x1
		y0, y1 = y1, y0 - q * y1
		r0, r1 = r1, r
	
	return [r1, x1, y1]
	
def inv(m, a):
    d, x, y = gcde(m, a)
    if d == 1:
    	return y % m
    	
    raise Exception, 'Error: gcd('+`m`+', '+`a`+') does not equal 1!'
    
def div(m, a, b):
	try:
		iab = inv(m, (a * b) % m) 
		return iab
	except Exception:
		#raise
		print 'Error: gcd('+`m`+', '+`a * b`+') does not equal 1!'
		
def expm(mp, ap, kp):
	m, a, k, r = mp, ap, kp, 1
	# Number of bits in k
	i = int(math.log(k, 2)) + 1
	while i >= 0:
		r = (r * r) % m
		# Is odd?
		if (k >> i) & 1:
			r = (r * a) % m
		i = i - 1
		
	return r
	
def factors(np):
	n = np
	
	fs, i = [], 2

	while i < n / i:
		while n % i == 0:
			fs.append(i)
			n = n / i
		i = i + 1
	if n > 1: 
		fs.append(n)
		
	return fs
	
def phi(np):
	n = np
	fs = factors(n)
	fsl, phi = len(fs), 1
	
	if not fsl:
		return 0
		
	for i in fs:
		phi = (i - 1) * phi 
		
	return phi
	
def fermat(n, t = 50):
	# Must handle error inputs and first four
	# positive integers.
	if n <= 0:
		return False # should be exception
	elif n <= 3:
		return True
	
	# No point wasting our time with even numbers
	if n % 2 == 0:
		return False

	while t:
		a = random.randint(2, n - 2)
		r = expm(n, a, n - 1)

		if r != 1:
			return False
		
		t = t - 1
		
	return True
	
def mr(n, t = 50):
	# No point in checking even numbers
	if n % 2 == 0:
		return 0
	# Rest is stolen ...		
	k, z = 0, n - 1

	# compute m,k such that (2**k)*m = n-1
	while not z % 2:
		k += 1
		z //= 2
	m = z

	# try tests with max random integers between 2,n-1
	ok = 1
	trials = 0
	p = 1
	while trials < t and ok:
		a = random.randint(2,n-1)
		trials += 1
		test = pow(a,m,n)
		if (not test == 1) and not (test == n-1):
			# if 1st test fails, fall through
			ok = 0
			for r in range(1,k):
				test = pow(a, (2**r)*m, n)
				if test == (n-1):
					ok = 1 # 2nd test ok
					break
		else: ok = 1  # 1st test ok
		if ok==1:  p *= 0.25
            
	if ok:  return 1 - p
	else:   return 0
	
		
def prime(d):
	while 1:
		r = random.getrandbits(d)
		
		#if fermat(r, 20):
		if mr(r, 20) > .99999999:
			return r

def safeprime(d):
	while 1:
		p = prime(d)
		p2p1 = (2 * p) + 1
		# I remeber seeing specific criteria to test 
		# 2p + 1 primes which may be more efficient
		# than another mr() 
		if mr(p2p1, 20) > .99999999:
			return p2p1