Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question

1 Answer 1

up vote 0 down vote accepted

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.