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

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.