# 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)

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/2) * pow(a, (t+1)/2, p)) % p
    return x, p-x

testQuadRes(100)

failure 2 17 1
failure 2 17 16
failure 4 17 1
failure 4 17 16
failure 8 17 1
failure 8 17 16
failure 9 17 1
failure 9 17 16
failure 13 17 1
failure 13 17 16
failure 15 17 1
failure 15 17 16
failure 16 17 1
failure 16 17 16
failure 2 41 31
failure 2 41 10
failure 4 41 18
failure 4 41 23
failure 5 41 23
failure 5 41 18
failure 8 41 25
failure 8 41 16
failure 9 41 40
failure 9 41 1
failure 20 41 4
failure 20 41 37
failure 21 41 4
failure 21 41 37
failure 23 41 10
failure 23 41 31
failure 25 41 37
failure 25 41 4
failure 31 41 16
failure 31 41 25
failure 32 41 40
failure 32 41 1
failure 33 41 25
failure 33 41 16
failure 36 41 23
failure 36 41 18
failure 39 41 31
failure 39 41 10
failure 40 41 1
failure 40 41 40
failure 3 73 64
failure 3 73 9
failure 6 73 4
failure 6 73 69
failure 9 73 8
failure 9 73 65
failure 12 73 55
failure 12 73 18
failure 18 73 37
failure 18 73 36
failure 19 73 32
failure 19 73 41
failure 23 73 16
failure 23 73 57
failure 24 73 8
failure 24 73 65
failure 25 73 37
failure 25 73 36
failure 27 73 1
failure 27 73 72
failure 35 73 2
failure 35 73 71
failure 36 73 16
failure 36 73 57
failure 38 73 2
failure 38 73 71
failure 41 73 55
failure 41 73 18
failure 46 73 1
failure 46 73 72
failure 48 73 37
failure 48 73 36
failure 49 73 8
failure 49 73 65
failure 50 73 16
failure 50 73 57
failure 54 73 32
failure 54 73 41
failure 57 73 4
failure 57 73 69
failure 61 73 55
failure 61 73 18
failure 65 73 64
failure 65 73 9
failure 67 73 4
failure 67 73 69
failure 69 73 2
failure 69 73 71
failure 70 73 64
failure 70 73 9
failure 71 73 32
failure 71 73 41
failure 72 73 1
failure 72 73 72
failure 5 89 9
failure 5 89 80
failure 9 89 53
failure 9 89 36
failure 10 89 42
failure 10 89 47
failure 11 89 73
failure 11 89 16
failure 17 89 69
failure 17 89 20
failure 18 89 10
failure 18 89 79
failure 20 89 18
failure 20 89 71
failure 21 89 40
failure 21 89 49
failure 22 89 44
failure 22 89 45
failure 25 89 81
failure 25 89 8
failure 34 89 55
failure 34 89 34
failure 36 89 17
failure 36 89 72
failure 40 89 84
failure 40 89 5
failure 42 89 68
failure 42 89 21
failure 44 89 57
failure 44 89 32
failure 47 89 21
failure 47 89 68
failure 49 89 5
failure 49 89 84
failure 50 89 22
failure 50 89 67
failure 53 89 72
failure 53 89 17
failure 55 89 34
failure 55 89 55
failure 57 89 11
failure 57 89 78
failure 68 89 49
failure 68 89 40
failure 69 89 71
failure 69 89 18
failure 71 89 79
failure 71 89 10
failure 72 89 20
failure 72 89 69
failure 73 89 85
failure 73 89 4
failure 79 89 47
failure 79 89 42
failure 80 89 36
failure 80 89 53
failure 81 89 50
failure 81 89 39
failure 84 89 80
failure 84 89 9
failure 85 89 87
failure 85 89 2
failure 87 89 25
failure 87 89 64
failure 88 89 88
failure 88 89 1
failure 2 97 49
failure 2 97 48
failure 3 97 65
failure 3 97 32
failure 4 97 73
failure 4 97 24
failure 6 97 81
failure 6 97 16
failure 8 97 85
failure 8 97 12
failure 9 97 54
failure 9 97 43
failure 11 97 53
failure 11 97 44
failure 12 97 89
failure 12 97 8
failure 16 97 91
failure 16 97 6
failure 18 97 27
failure 18 97 70
failure 22 97 75
failure 22 97 22
failure 24 97 93
failure 24 97 4
failure 25 97 66
failure 25 97 31
failure 27 97 18
failure 27 97 79
failure 31 97 72
failure 31 97 25
failure 32 97 94
failure 32 97 3
failure 33 97 50
failure 33 97 47
failure 36 97 62
failure 36 97 35
failure 43 97 88
failure 43 97 9
failure 44 97 86
failure 44 97 11
failure 47 97 64
failure 47 97 33
failure 48 97 95
failure 48 97 2
failure 49 97 2
failure 49 97 95
failure 50 97 33
failure 50 97 64
failure 53 97 11
failure 53 97 86
failure 54 97 9
failure 54 97 88
failure 62 97 36
failure 62 97 61
failure 64 97 47
failure 64 97 50
failure 65 97 3
failure 65 97 94
failure 66 97 25
failure 66 97 72
failure 70 97 79
failure 70 97 18
failure 72 97 31
failure 72 97 66
failure 73 97 4
failure 73 97 93
failure 75 97 22
failure 75 97 75
failure 79 97 70
failure 79 97 27
failure 81 97 6
failure 81 97 91
failure 85 97 8
failure 85 97 89
failure 86 97 44
failure 86 97 53
failure 88 97 43
failure 88 97 54
failure 89 97 12
failure 89 97 85
failure 91 97 16
failure 91 97 81
failure 93 97 24
failure 93 97 73
failure 94 97 32
failure 94 97 65
failure 95 97 48
failure 95 97 49
failure 96 97 96
failure 96 97 1