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.
Programming Praxis 
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.
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