#!/usr/bin/env sage -python
#
# SolovayStrassen: Performs the Solovay-Strassen primality test
#
# Stephen Forrest, http://wandership.ca/

import sys 
from sage.all import * 

def ilog(x, base=math.e):
    i = 0
    if x < 1 and base < 1:
        raise ValueError, "logarithm cannot compute"
    while x < 1:
        i -= 1
        x *= base
    while x >= base:
        i += 1
        x /= base
    return i

def SolovayStrassen(n,err):
    b = -ilog(err,2)
    R = IntegerModRing(n)
    for i in range( b ):
        m = randint(2,n-1)
        if gcd(m,n) > 1 or R( kronecker_symbol(m,n) ) <> R(m) ** ((n-1)/2):
            return(false)
    return(true)

# Examples
print SolovayStrassen(1999, 0.0001)
print SolovayStrassen(2001, 0.0001)