Implementing AKS Primality Prover

I am a computer programmer interested in prime numbers. I have implementations of several algorithms related to prime numbers at my blog. I want to add an implementation of the AKS primality prover to my collection, but I am having trouble, and my knowledge of math is insufficient to make sense of some of the things I read, so I ask for help here. I have two questions: how to compute the exponentiation of a polynomial mod another polynomial and an integer, and how to compute the r in the AKS algorithm. I begin with polynomial exponentiation, using an example.

Consider the problem of squaring polynomial x^3 + 4 x^2 + 12 x + 3 modulo (x^5 - 1, 17). Polynomial multiplication is exactly the same as grade-school multiplication, except there are no carries, so the process looks like this:

                 1    4   12    3
1    4   12    3
---  ---  ---  ---
3   12   36    9
12   48  144   36
4   16   48   12
1    4   12    3
---  ---  ---  ---  ---  ---  ---
1    8   40  102  168   72    9


Thus, 1 x^3 + 4 x^2 + 12 x + 3 squared is 1 x^6 + 8 x^5 + 40 x^4 + 102 x^3 + 168 x^2 + 72 x + 9. To reduce the result modulo 1 x^5 - 1 we divide by the grade-school long-division algorithm and take the remainder, which gives 1 x + 8 with a remainder of 40 x^4 + 102 x^3 + 168 x^2 + 73 x + 17:

                                          1    8
+ ---  ---  ---  ---  ---  ---  ---
1 0 0 0 0 -1 |   1    8   40  102  168   72    9
1    0    0    0    0   -1
---  ---  ---  ---  ---  ---
8   40  102  168   73    9
8    0    0    0    0   -8
---  ---  ---  ---  ---  ---
40  102  168   73   17


We can confirm the calculation by multiplying and adding the remainder:

       1    0    0    0    0   -1
1    8
---  ---  ---  ---  ---  ---
8    0    0    0    0   -8
1    0    0    0    0   -1
---  ---  ---  ---  ---  ---  ---
1    8    0    0    0   -1   -8
40  102  168   73   17
---  ---  ---  ---  ---  ---  ---
1    8   40  102  168   72    9


Then we simply reduce each of the coefficients of the remainder modulo 17, giving the result 6 x^4 + 0 x^3 + 15 x^2 + 5 x + 0. The whole computation can be organized as shown below. Note how the division and reduction modulo x^5 - 1 is accomplished, eliminating the high-order coefficients and adding them back to the low-order coefficients; we are relying here on the simple form of the polynomial modulus, and would have to do more work if it was more complicated:

                 1    4   12    3   multiplicand
1    4   12    3   multiplier
---  ---  ---  ---
3   12   36    9    3 * 1 4 12 3 * x^0
12   48  144   36        12 * 1 4 12 3 * x^1
4   16   48   12              4 * 1 4 12 3 * x^2
1    4   12    3                   1 * 1 4 12 3 * x^3
---  ---  ---  ---  ---  ---  ---
1    8   40  102  168   72    9   1 4 12 3 * 1 4 12 3
-1   -8                   1    8   reduce modulo x^5 - 1
---  ---  ---  ---  ---  ---  ---
40  102  168   73   17   reduce modulo 17
6    0   15    5    0   final result


Now that we can perform modular multiplication of a polynomial, we return to the modular exponentiation of a polynomial, which is done using the same square-and-multiply algorithm as modular exponentiation of integers, except that we use polynomial modular multiplication instead of integer modular multiplication. Here is the algorithm on integers; to adapt it to polynomials, just replace the integer modular multiplications with polynomial modular multiplications.

func powermod(b, e, m) # b^e (mod m)
r := 1
while e > 0
if e is odd
r := r * b (mod m)
e := floor(e/2)
b := b * b (mod m)
return r


So, the first question: Have I correctly stated the algorithm for modular exponentiation of polynomials?

For the second question, I start with the statement of the AKS algorithm, as given at the Prime Pages:

Input: Integer n > 1
Output: either PRIME or COMPOSITE

if (n has the form a^b with b > 1)
then output COMPOSITE

r := 2

while (r < n) {

if (gcd(n,r) is not 1)
then output COMPOSITE

if (r is prime greater than 2)
then {

let q be the largest factor of r-1

if (q > 4 * sqrt(r) * log2 n)
and ( n^{(r-1)/q} is not 1 (mod r) )
then break

}

r := r+1

}

for a = 1 to 2 * sqrt(r) * log2 n {

if ( (x-a)^n is not (x^n-a) (mod x^r-1,n) )
then output COMPOSITE

}

output PRIME;


Again I will work with a specific example, trying to prove that n = 89 is prime. We first consider if n is a perfect power of the form a^b with b > 1. When a = 2, the powers of 2 are 4, 8, 16, 32, 64 and 128, so 2 fails. When a = 3, the powers of 3 are 9, 27, 81 and 243, so 3 fails. When a = 5, the powers of 5 are 25 and 125, so 5 fails. When a = 7, the powers of 7 are 49 and 343, so 7 fails. When a = 11, or any higher prime, a^2 is greater than 89, so we are finished with the perfect power tests.

Second we compute the value of r. Since r and q = (r-1)/2 must both be prime, and q must be greater than 4 * log n = 26, the only possibility is r = 59, but 4 * sqrt(59) * log2(89) = 199, so there is no early break from the loop, and r must be 89.

So, the second question: Have I correctly computed r = 89?

I think that must be incorrect. With r = 89, a will range from 1 to 122. Let's take a = 17 as an example. Modular exponentiation of the polynomial (x - 17) ^ 89 (mod x^89 - 1, 89) is 73; that is, the number 73, with all coefficients of powers of x equal to 0. But that doesn't equal x^89 - 17, suggesting that 89 is composite. Of course, 89 is prime, so something is wrong.

Code is available at http://codepad.org/4jwgScdX. Please let me know what I have done wrong. Thank you for reading this far.

-