1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/examples/prime.py	Tue Feb 22 18:58:21 2011 +0100
     1.3 @@ -0,0 +1,74 @@
     1.4 +isProbPrime = True
     1.5 +
     1.6 +import sys
     1.7 +from math import log
     1.8 +from random import randint
     1.9 +
    1.10 +def main():
    1.11 +	args = parse_args()
    1.12 +	n = args.n
    1.13 +	if args.s:
    1.14 +		if safe_is_prime(n):
    1.15 +			print "%i is most probably a prime number." % n
    1.16 +		else:
    1.17 +			print "%i is not a prime number." % n
    1.18 +		return 
    1.19 +	elif args.a:
    1.20 +		a = args.a
    1.21 +	else:
    1.22 +		a = randint(2, n-1)
    1.23 +	if is_prime(n, a):
    1.24 +		print "%i is probably a prime." % n	
    1.25 +	else:
    1.26 +		print "%i is not a prime." % n
    1.27 +
    1.28 +def safe_is_prime(n):
    1.29 +	aset = set()
    1.30 +	while len(aset) < int(log(n)):
    1.31 +		aset.add(randint(2, n-1))
    1.32 +	for a in aset:
    1.33 +		if not is_prime(n, a):
    1.34 +			return False
    1.35 +	return True
    1.36 +	
    1.37 +def is_prime(n, a=None):
    1.38 +	if not a:
    1.39 +		a = randint(2, n-1)
    1.40 +	global isProbPrime
    1.41 +	isProbPrime = True
    1.42 +	#print "Prime Test for %i with a = %i:" % (n, a)
    1.43 +	r = power(a, n-1, n)
    1.44 +	if r != 1 or not isProbPrime:
    1.45 +		#print "=> %i is not a prime." % n
    1.46 +		return False
    1.47 +	else:
    1.48 +		#print "=> %i is probably a prime." % n
    1.49 +		return True
    1.50 +
    1.51 +def power(a, p, n):
    1.52 +	global isProbPrime	
    1.53 +	if p == 0:
    1.54 +		return 1
    1.55 +	x = power(a, p/2, n)
    1.56 +	r = (x*x)%n
    1.57 +	#print "power(%i, %i, %i)" % (a, p, n)
    1.58 +	#print "result = %i" % r
    1.59 +	if r == 1 and x != 1 and x != n-1:
    1.60 +		isProbPrime = False
    1.61 +		#print "isProbPrime = False"
    1.62 +	if p%2 == 1:
    1.63 +		r = (a*r)%n
    1.64 +		#print "updated result = %i" % r
    1.65 +	return r
    1.66 +
    1.67 +from argparse import ArgumentParser
    1.68 +
    1.69 +def parse_args():
    1.70 +    parser = ArgumentParser(description='Applies the randomised prime number test.')
    1.71 +    parser.add_argument('n', type=int, help='number to be tested')
    1.72 +    parser.add_argument('-a', type=int, help='provide a yourself')
    1.73 +    parser.add_argument('-s', action='store_true', help='multiple tests to minimise probability of failure')
    1.74 +    return parser.parse_args()
    1.75 +
    1.76 +if __name__ == "__main__":
    1.77 +    main()