Modular Testing

November 19, 2013

Testing prime-number functions is easy, if you take the time to do it, because the functions can be checked against each other. We begin by testing primes against isPrime:

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
            if isPrime(i):
                print "isPrime not prime", i

Running the test causes noting to be reported, as expected, on the theory that “no news is good news:”

>>> testPrime(200)

For any x, the inverse of x (mod m) exists and is equal to 1 whenever x is not coprime to m. That mathematical fact makes it easy to test the modular inverse function:

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

Again, no news is good news:

>>> testInverse(200)

The third test looks at the jacobi symbol and modular square root functions. If a number a is a quadratic residue to an odd prime modulus p, the jacobi symbol will be 1 and the modular square root x will satisfy the equation (x * x) % p = a:

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

When we run the test with testQuadRes(100), a whole cascade of error messages appears. The errors occur when p is 17, 41, 73, 89 and 97; since all of those numbers are 1 (mod 8), it is obvious that the problem is with the default case of modSqrt. Indeed, the next-to-last line is in error — m should be divided by 2 — and fixing the error makes everything work properly:

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)

In my quick testing at the console, I had not tested any p ≡ 1 (mod 8), and thus missed the bug. It was hard to find but easy to fix, and would never have been an issue if I had tested properly in the first place.

By the way, such testing doesn’t prove that the function is bug-free, though it does make me feel better about my work. Challenge: can anybody spot any remaining bugs?

You can run the program at

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:
                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: Logo

You are commenting using your 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


Get every new post delivered to your Inbox.

Join 728 other followers

%d bloggers like this: