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

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.

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.