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

-

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.

-

The calculation of r was incorrect; it should be 191.

wait what?, that is wrong the r for 89 is not 191, is 43

first star with the definition of r in the context of the AKS test:

find the smallest r such that such that: ordr(n) > (log2(n))2

where ordr(n) is the multiplicative order of n modulo r, and that is: if gcd(n,r)=1, then is smallest k>=1 such that nk=1 (mod r)

programing that is very easy, here are the first ordr(89) for r>=2:

ord2(89)=1, ord3(89)=2, ord4(89)=1, ord5(89)=2, ord6(89)=2, ord7(89)=6, ord8(89)=1, ord9(89)=2, ord10(89)=2, ord11(89)=1, ord12(89)=2, ord13(89)=12, ord14(89)=6, ord15(89)=2, ord16(89)=2, ord17(89)=4, ord18(89)=2, ord19(89)=18, ord20(89)=2, ord21(89)=6, ord22(89)=1, ord23(89)=22, ord24(89)=2, ord25(89)=10, ord26(89)=12, ord27(89)=6, ord28(89)=6, ord29(89)=28, ord30(89)=2, ord31(89)=10, ord32(89)=4, ord33(89)=2, ord34(89)=4, ord35(89)=6, ord36(89)=2, ord37(89)=36, ord38(89)=18, ord39(89)=12, ord40(89)=2, ord41(89)=40, ord42(89)=6, ord43(89)=42, ord44(89)=1, ord45(89)=2, ord46(89)=22, ord47(89)=23, ord48(89)=2, ord49(89)=42, ord50(89)=10, ord51(89)=4, ord52(89)=12, ord53(89)=13, ord54(89)=6, ord55(89)=2, ord56(89)=6, ord57(89)=18, ord58(89)=28, ord59(89)=58, ord60(89)=2

now (log2(89))2 = ( 6.4757 )2 = 41.9351

then r=43 because ord43(89)=42

Now, the other problem I see is when you explain how you check for n=ab for b>1, you said that you check only primes power, so unless that you do something else behind the scene, that is wrong too.

Take for example this number 8281, there is no prime such that pb=8291 for b>=1, but 8281 is a perfect square, it is 912 (91=7*13)

for a more accurate check you can calculate the b-th root of n, n1/b, round that and ask if that to the power of n is equal to n, that is: round(n1/b)b=n

About the rest I don't know, my knowledge about polynomials operations is not enough because I am also a programmer that want to add the AKS test to my collection...

-