Solovay-Strassen primality testing with Python.

solovay_strassen.py

def solovay_strassen(n, k=10):
	if n == 2:
		return True
	if not n & 1:
		return False

	def legendre(a, p):
		if p < 2:
			raise ValueError('p must not be < 2')
		if (a == 0) or (a == 1):
			return a
		if a % 2 == 0:
			r = legendre(a / 2, p)
			if p * p - 1 & 8 != 0:
				r *= -1
		else:
			r = legendre(p % a, a)
			if (a - 1) * (p - 1) & 4 != 0:
				r *= -1
		return r

	for i in xrange(k):
		a = randrange(2, n - 1)
		x = legendre(a, n)
		y = pow(a, (n - 1) / 2, n)
		if (x == 0) or (y != x % n):
			return False

	return True

# benchmark of 10000 iterations of solovay_strassen(100**10-1); Which is not prime.

#10000 calls, 2440 per second.
#571496 function calls (74873 primitive calls) in 4.100 seconds

can you please explain line 14,18 and 26

Check on num 10493450440980479677. It is Prime, but test says that it is composite

Check on num 10493450440980479677. It is Prime, but test says that it is composite

For large numbers, use the decimal library:

def solovay_strassen(n, k=10): from decimal import Decimal n, k = Decimal(n), Decimal(k)