Modular Testing

November 19, 2013

Today’s exercise is my reminder to myself to do a better job of testing. I was porting some prime-number code to Python, one of several preliminary steps to writing a new essay. I am comfortable enough with Python (mostly), and there is nothing particularly tricky about the code, so I wrote the code, tested it quickly, and went on. Later, a different part of the code was failing, and I couldn’t find the problem. Of course, the problem was in the earlier code that I had quickly tested, and therefore hard to spot, since I was looking in the wrong place. The failing code is shown below.

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

Your task is to write a proper test suite, discover the bug in the code shown above, and fix it. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Follow

Get every new post delivered to your Inbox.

Join 645 other followers