## 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         else:             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 http://programmingpraxis.codepad.org/977M1PA4.

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
```