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

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