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.
Best Answer
Ok, I've got it. The calculation of r was incorrect; it should be 191. Working code, including the full AKS prover, is at http://programmingpraxis.codepad.org/6ZHrsEmx. I'll have a full write-up at my blog later this week: Part 1 and Part 2.