codepad

Language:

# modular testing

def gcd(m,n):
    while n <> 0:
        m, n = n, m % n
    return m

def primes(n):
    ps, sieve = [], [True] * (n + 1)
    for p in range(2, n + 1):
        if sieve[p]:
           ps.append(p)
           for i in range(p * p, n + 1, p):
               sieve[i] = False
    return ps

def isPrime(n):
    if n % 2 == 0:
        return n == 2
    d = 3
    while d * d <= n:
        if n % d == 0:
            return False
        d = d + 2
    return True

def inverse(x, m):
    a, b, u = 0, m, 1
    while x > 0:
        q = b // x
        x, a, b, u = b % x, u, x, a - q * u
    if b == 1: return a % m
    raise ArithmeticError("must be coprime")

def jacobi(a, p):
    a = a % p; t = 1
    while a != 0:
        while a % 2 == 0:
            a = a / 2
            if p % 8 in [3, 5]:
                t = -t
        a, p = p, a
        if a % 4 == 3 and p % 4 == 3:
            t = -t
        a = a % p
    return t if p == 1 else 0

def modSqrt(a, p):
    a = a % p
    if p % 4 == 3:
        x = pow(a, (p+1)/4, p)
        return x, p-x
    if p % 8 == 5:
        x = pow(a, (p+3)/8, p)
        c = pow(x, 2, p)
        if a != c:
            x = (x * pow(2, (p-1)/4, p)) % p
        return x, p-x
    d = 2
    while jacobi(d, p) != -1:
        d += 1
    t, s = p-1, 0
    while t % 2 == 0:
        t, s = t / 2, s + 1
    at, dt, m = pow(a, t, p), pow(d, t, p), 0
    for i in range(0, s):
        if pow(at * pow(dt, m), pow(2, s-i-1), p) == p - 1:
            m = m + pow(2, i)
    x = (pow(dt, m) * pow(a, (t+1)/2, p)) % p
    return x, p-x

def testPrime(n):
    ps = primes(n)
    for i in range(2,n):
        if i in ps:
            if not isPrime(i):
                print "prime not isPrime", i
        else:
            if isPrime(i):
                print "isPrime not prime", i

testPrime(200)

def testInverse(n):
    for m in range(2, n):
        for i in range(0, m):
            if gcd(i,m) == 1:
                j = inverse(i,m)
                if (i*j)%m <> 1:
                    print "inverse fails", i, m

testInverse(100)

def testQuadRes(n):
    for p in primes(n)[1:]:
        for a in range(1,p):
            if jacobi(a,p) == 1:
                x, y = modSqrt(a,p)
                if (x*x)%p <> a:
                    print "failure", a, p, x
                if (y*y)%p <> a:
                    print "failure", a, p, y

testQuadRes(100)

Private [?]

Run code