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.

About these ads

Pages: 1 2

2 Responses to “Modular Testing”

  1. phillip said

    def modSqrt(a, p):
    a = a %p
    if a ==0:
    return 0,p
    i = 1
    while i <p/2:
    x =i*i%p
    if x == a:
    return x,p-x
    i += 1

    not sure if this is the error, but probably equally as fast as your code.

  2. Paul said

    Here some tests for all functions. The modSqrt function fails. All other seem OK. A lot of random testing is used and for the inverse functions (inverse and modSqrt) a simple check is done, if the function does, what it promises. The example cases for the jacobi symbol were gathered on various places on the web.

    import prprpr as PR
    from random import randrange, choice
    from math import sqrt
    
    primes = PR.primes(1000)
    
    def test_primes():
        act = PR.primes(100)
        exp = map(int, """2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 
            73 79 83 89 97""".split())
        assert act == exp
    
    def test_isPrime():
        N = 100000
        MAX = 1000000
        def is_prime(n):
            mx_prime = sqrt(n)
            for p in primes:
                if p > mx_prime:
                    break
                if n % p == 0:
                    return False
            return True
        for i in xrange(N):
            n = randrange(2, MAX - 1)
            act = PR.isPrime(n)
            exp = is_prime(n)
            assert act == exp, "i={:d} n={:d} act={:d} exp={:d}".format(i, n, 
                                                                act, exp) 
    
    def test_inverse():
        N = 1000
        for i in xrange(N):
            p = choice(primes[1:])
            a = randrange(1, p - 1)
            act = PR.inverse(a, p)
            assert (act * a) % p == 1
    
    def test_jacobi():
        cases = [(1001, 9907, -1), (222222, 304679, -1), (222222, 324899, 1),
                 (222222, 333333, 0), (1411, 317, -1), (1735, 507, 1)]
        for a, p, exp in cases:
            act = PR.jacobi(a, p)  
            assert exp == act
        
        
    def test_modSqrt():  # --> FAILS
        N = 1000
        for i in xrange(N):
            p = choice(primes)
            a = randrange(1, p - 1)
            x, y = PR.modSqrt(a, p)
            assert (x * x) % p == a, "i={:d} a={:d} p={:d} x={:d} y={:d}".format(i,
                                         a, p, x, y)
            assert (y * y) % p == a
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 629 other followers

%d bloggers like this: