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.