is there a major inefficiency in this miller-rabin pseudocode in CLRS?

Question

This question may actually have nothing to do with the Miller-Rabin primality testing procedure; it may only easy analysis of some simple pseudocode.

On p 969 of CLRS (Introduction to Algorithms 3ed), an auxiliary function to Miller-Rabin is presented:

WITNESS(a, n)
    let t and u be such that t >= 1, u is odd, and n-1 = 2^t u
    x_0 = MODULAR-EXPONENTIATION(a, u, n)
    for i = 1 to t
        x_i = x_{i-1}^2 mod n
        if x_i == 1 and x_{i-1} != 1 and x_{i-1} != n-1
            return TRUE
    if x_t != 1
        return TRUE
    return FALSE

I copied the above exactly from the textbook.

Now, knowing only that MODULAR-EXPONENTIATION returns a result between 0 and n-1, inclusive, I think the pseudocode above is entirely equivalent to

WITNESS(a, n)
    let t and u be such that t >= 1, u is odd, and n-1 = 2^t u
    x_0 = MODULAR-EXPONENTIATION(a, u, n)
    if x_0 == 1 or x_0 == n-1
        return FALSE
    else
        return TRUE

If so, there's probably something else wrong with the original implementation, since if I'm not mistaken Miller-Rabin witnessing does require some sort of looping. Can someone provide a simple counterexample to show that I am wrong?

Answer 1

The Miller-Rabin primality test is designed to be TRUE for n being a prime, so returning FALSE should only apply to composite numbers. Let's test this with a little Python program.

def wrongwitness(a, n):             #implementation of your shortcut
    u = n - 1
    t = 0
    while u % 2 == 0:               #n - 1 = 2^t * u
        u //= 2
        t += 1

    x_0 = pow(a, u, n)              #x0 = a ^ u (mod n), oops, where is t?

    if x_0 == 1 or x_0 == n - 1:
        return False
    else:
        return True

primes = [5, 7, 11, 13, 17, 19, 23, 29, 31]

for p in primes:         
    for a in range(2, p):           #1 < a < p
        if not wrongwitness(a, p):  #witness returned FALSE, though we have a prime number
            print("Found counter example: a = ", a, "and p = ", p )

This gives us lots of counterexamples for your shortcut implementation as small as a = 2 and p = 5 or a = 3 and p = 7. Actually all (p - 1, p) tuples are counter examples. So no shortcut, you have to test all square roots of a^(n-1) as explained in your text book.

P.S.: But there are ways to reduce the number of calculations, you have to perform. Subsets of witnesses have been identified for n up to 3,317,044,064,679,887,385,961,981. So for n < 1,373,653 it is for instance sufficient just to test a=2 and a=3.