AKS prmality test

[killo]

I have been looking up the AKS primality test and would like to implement it in code. The algorithm is here: http://teal.gmu.edu/courses/ECE746/project/F06_Project_resources/Salembier_Southerington_AKS.pdf
on page 2. The improved algorithm.
I understand most of the steps but some I don't, in particular:
step 4:
for a = 1 to sqrt(r)log(n)
if (( x + a)^n = xn + a (mod f(x), n)), output PRIME
next

I don't know what it means by having two numbers in mod() with a comma between them but I would know what it meant if there was only one number, it would mean that both sides produce the same remainder when divided by that number.
Also it doesn't state a value for x. Later on in the paper it says about how the equation is some kind of polynomial so that a value for x is not needed, but I don't understand how to implement this in code.
I showed the paper to one of my maths teachers but he didn't know what was meant by two numbers in the mod() function and didn't know how to go about proving that one side of the equation in step 4 was equal to the other side.
I'd appreciate any help in this.

[jemidiah]

You're probably running up against abstract algebra you've never heard of. The notation used is saying the left and right sides are equivalent in the quotient ring Z_n[x]/(f(x))....

More practically, though, it's saying this:

1. Compute (x+a)^n using a variant of Exponentiation by Squaring (http://en.wikipedia.org/wiki/Exponentiation_by_squaring).

That is, for a particular number a, figure out the constants in front of each power of x. However, using standard exponentiation by squaring will take too long. At the end, you'll have computed n coefficients, while the whole algorithm is supposed to take polynomial time in the *length* of n (written as, say, a binary number). So, if you use the standard algorithm, your run time will be at least O(n) instead of O((logn)^c) for some constant c, which is much too slow.

The trick, then, is to compute only a few terms. This can be done by performing the arithmetic modulo the polynomial x^r - 1 (this wasn't chosen out of a hat, at all, but the reason it was chosen is technical). That is, at any point in your intermediate squaring steps you get a polynomial at or above degree r, feel free to just subtract off any multiple of (x^r-1) to get rid of the large terms. For example, if your intermediate step gave

3x^9 + 4x^2 + x - 1

and you had previously found r=7, you can subtract off

(x^7-1)*3x^2 = 3x^9-3x^2

to get

x^2 + x - 1

You can continue your exponentiation by squaring from there. The point is to keep your intermediate polynomials short--in this case, of degree <= r.

2. Reduce x^n+a so that it's at most degree r.

Do so by subtracting off multiples of (x^r-1) like you did above until the resulting polynomial has a small enough degree (<= r). Since r is kept small, and there are around n/r repetitions of this step, it's likely this way to compute this step has a large run time order. In fact, from the AKS paper, it looks like r = O(log^3 n), so this step would definitely take too long, being almost linear in time. So, an optimization is needed. Computing a few steps "by hand" gives it to us:

x^n+a - (x^r-1)x^(n-r) = x^n+a - x^n + x^(n-r) = x^(n-r)+a
x^(n-r)+a - (x^r-1)x^(n-2r) = x^(n-2r)+a
...
x^(n-qr)+a

where q is the largest integer such that qr < n. Put another way, q is the largest integer such that q < n/r, also known as the floor of n/r (unless n/r is an integer, in which case it's 1 less than the floor). Since this can be computed very easily, this step can indeed be performed quickly.

3. Check to see if the coefficients you have from (1) and (2) are equal or not, which says if the equation is true or false.


In all of this, every coefficient is an integer mod n. That is, whenever you calculate a coefficient, take that number mod n. This keeps the coefficients small in length, which is another requirement for the overall algorithm being fast.

I wouldn't want to put myself through the pain of implementing that algorithm, but if you're really interested in it, go for it, and good luck!

[killo]

Wow thanks for explaining that to me :) I'll carefully look at this over the next few days and tell you how I get on.
Some preliminary questions:

Where you said:

and you had previously found r=7, you can subtract off
(x^7-1)*3x^2 = 3x^9-3x^2

where did you get the 3x^2 from? did you just choose that to get the polynomial up to the 9th order? How did you go about choosing it?

So now, my general understanding is that you need to keep on subtracting off (x^r-1) from both sides of the equation until they are as small as you can get them, and then once that's done, take both sides mod n and then compare the two sides to see if they're equal?

(Also I have no Idea what a quotient ring is, but it doesn't seem like I need to understand that)

[jemidiah]

Where you said:

and you had previously found r=7, you can subtract off
(x^7-1)*3x^2 = 3x^9-3x^2

where did you get the 3x^2 from? did you just choose that to get the polynomial up to the 9th order? How did you go about choosing it?

Yup, I chose 3x^2 so that the resulting polynomial would have degree <= r. The generalized version of modulus used in this algorithm uses polynomials instead of integers. Under modular arithmetic, two numbers that differ by a multiple of the modulus are considered equal (eg. 4 == 10 mod 3, since they differ by 10-4=6=2*3). In this polynomial version, two polynomials that differ by a multiple of the modulus are considered equal (eg. as above 3x^9 + 4x^2 + x - 1 == x^2 + x - 1 mod x^7-1, since they differ by 3x^9 - 3x^2=3x^2(x^7-1)). Your goal ultimately is to reduce both sides of that equation by subtracting off multiples of x^7-1, just like you'd do with normal modular arithmetic, so you can finally tell if they're equal or not. In this case, you've reduced "enough" when both sides are of degree <= the degree of the modulus.

In general, given a polynomial of degree d greater than r with leading term a, subtract off ax^(d-r)(x^r-1), which will get rid of the highest order term:

ax^d+(lower order terms) - ax^(d-r)(x^r-1)
= ax^d+(lower order terms) - ax^d + ax^(d-r)
= (lower order terms) + ax^(d-r)

You may need to repeat this several times to get the result of degree <= r.

So now, my general understanding is that you need to keep on subtracting off (x^r-1) from both sides of the equation until they are as small as you can get them, and then once that's done, take both sides mod n and then compare the two sides to see if they're equal?

That's pretty close. The repeated squaring I mentioned above is very important in efficiently evaluating the left hand side. Put another way, without repeated squaring you'd be better off dividing n by each integer up to Sqrt(n) to test if it's prime. At each step of repeated squaring you keep the polynomials you have to work with small by taking them mod x^r-1 using reduction steps like the one above.

During the course of evaluating the left hand side, you'll almost certainly get coefficients at or larger than n. The binomial theorem says how big the coefficients might get--(n choose i) for coefficient of x^i--though they'll probably stay smaller since that's not taking into account that we're taking this modulo x^r-1. Those coefficients are calculated in part using n!, which, if n is enormous, would be *gigantic*. Anywho, to keep the numbers small, at each step you compute a new coefficient you have to take it mod n. You could instead do this at the very end, that is, take each coefficient mod n and compare left and right hand side polynomials to see if they're equal--to see if they have the same coefficients--but this is most likely wildly inefficient since in intermediate steps you will have had to do arithmetic on insanely large numbers.


(Also I have no Idea what a quotient ring is, but it doesn't seem like I need to understand that)

I'm trying to hide the details that aren't relevant to actually computing that equivalence. To be honest, I don't know that many of them; I'm not great at abstract algebra, and would need to look up/prove a number of preliminary results to be satisfied I'm ready to understand the reasoning in the AKS algorithm's proof.


Edit: it occurred to me that I meant < r instead of <= r basically every time I've used that phrase in the last two posts.

[jemidiah]

Another thing occurred to me. When I subtracted off ax^(d-r)(x^r-1), effectively I was able to just lower the exponent of the ax^d term by r. I can do this as many times as I want, so long as the exponent stays >= 0. This is basically just a roundabout way of doing modular arithmetic, this time in the exponent, using modulus r. So, you can entirely ignore the whole "subtract off a polynomial" thing in this particular case and just take your exponents mod r.

So, succinctly, to check for equality,
1. Evaluate (x+a)^n using repeating squaring. At each step take coefficients mod n and exponents mod r. The coefficients will be numbers since a will be a number.
2. Reduce x^n+a by taking the exponent mod r and the coefficients (1 and a) mod n (which probably does nothing).
3. Compare polynomials obtained from (1) and (2) by checking if the coefficients are all equal. If so, (x+a)^n = x^n+a [...in Z_n[x]/(x^r-1)...]; otherwise, (x+a)^n != x^n+a.

[killo]

Thanks, I almost understand this fully now, but I have 2 questions:
At each step take coefficients mod n and exponents mod r.
so if I had 23x^25 + 34x^12 + 12x^7
and n = 31, r = 7
Taking coefficients mod n: 23x^25 + 3x^12 + 12x^7
Taking exponents mod r: 23x^4 + 3x^5
(the 12x^7 term dissapears since 7 mod 7 = 0)
Is that right?

Also, the method for working out r in the paper, is:
· q > floor ((log n)^2)
· Compute n^j mod q for j = 1, 2,…,[(logn)^2]
· If residue equals 1 mod q then q++
· Else r = q

(I'm assuming they mean log to base 10)
Here's how I would work out r if n = 31

q = floor(log(31)^2) = floor(2.224...) = 2
31 mod 2 = 1 (so q++)
31 mod 3 = 1 (so q++)
31 mod 4 = 4
31^2 mod 4 = 1 (so q++)
31 mod 5 = 1 (so q++)
31 mod 6 = 1 (so q++)
31 mod 7 = 3
31^2 mod 7 = 2

so r = 7

Is that right?

[jemidiah]

(the 12x^7 term dissapears since 7 mod 7 = 0)

The 12x^7 term becomes 12x^(7 mod 7) = 12x^0 = 12, not 0. Otherwise, you're right. For the rest, I'll have to check it from the paper later, but I'd be a bit surprised if the logarithm is base 10. It's probably base 2 or e, but I'd have to hunt through the paper to see. The original paper might shed some light if this one doesn't.

[killo]

Oh yeah, of course it does, I forgot that x^0 is 1 and not 0.

On wikipedia it says that the logarithm is to base 2 http://en.wikipedia.org/wiki/AKS_primality_test
but that's the initial algorithm they come up with, we're looking at the improved one, so it might be different?

also, if it was log to base 2 then this doesn't seem to fit because:
Lets say I was checking if 31 is prime for example: then in step 2: log(31)^2 = 24.544, so in the best case scenario r = 24
Then in step 4 where it says: "for a = 1 to sqrt(r)log(n)" well sqrt(r)log(n) would equal 24, so that would be: for a = 1 to 24.
24 is alot bigger than sqrt(31) = 5.57, so I might aswell just do a search for factors up to sqrt(31), this is why I think log to base 2 doesn't fit.

It would make a similar problem if it was log to base e (since 2 and 2.7 aren't that far apart)

[jemidiah]

Well, log_2(31) might be vaguely near 31 causing the naive algorithm to outperform AKS, but log_2(120439812357613947561945619), for instance, is quite small compared to that long number. The AKS primality test is asymptotically (i.e. "eventually") much, much faster than the naive algorithm, since the naive algorithm is O(Sqrt(n)) multiplications and AKS is O(log(n)^c) for c <= 12 depending on optimizations. The big-O hides when the naive algorithm starts being slower than AKS, but it does unequivocally say that, at some point, AKS starts being strictly faster than the naive algorithm.

As for the rest of your previous post,
Also, the method for working out r in the paper, is:
· q > floor ((log n)^2)
· Compute n^j mod q for j = 1, 2,…,[(logn)^2]
· If residue equals 1 mod q then q++
· Else r = q

(I'm assuming they mean log to base 10)
Here's how I would work out r if n = 31

q = floor(log(31)^2) = floor(2.224...) = 2
31 mod 2 = 1 (so q++)
31 mod 3 = 1 (so q++)
31 mod 4 = 4
31^2 mod 4 = 1 (so q++)
31 mod 5 = 1 (so q++)
31 mod 6 = 1 (so q++)
31 mod 7 = 3
31^2 mod 7 = 2

so r = 7


The notation q > floor(log^2 n) means you should start at floor(log^2 n) + 1, not floor(log^2 n). The bullet points are, to be honest, not very clear or well done. They require an understanding of the multiplicative order of n, mod r, to interpret, which is circular. What they meant to write but wanted to be more concise than was...


f = floor(log^2 n)
q = f + 1
while True:
for j=1, ..., f
if n^j mod q = 1:
break // i.e. skip the rest of the for loop; this q didn't work
end if
next
if j=f+1: // i.e. the for loop completed without a break
break // i.e. exit the do loop; q worked
q++
do
r = q


In words, this is finding the smallest number r, where r > log^2 n, such that r^i is not 1 for each i between 1 and log^2 n.

As for the base of the logarithm, the original paper uses log base 2. The revised paper continues to use log base 2. The graphs at the end of the implementation paper you linked are in "number of bits" which is basically same as log n, where log is base 2. I looked at the LiDIA documentation which says log is base e, but without more code than that shown in Figure 1 I'd say they simply wrapped that function to make it log base 2.

[killo]

Thanks for your help so far, I've started implementing this, but even the first step runs slower than the naive method:

Naive method (Trial division):
If CheckingNum Mod 2 = 0 Then goto Done 'Composite
sqrtcheckingnum = checkingnum ^ 0.5
i = 3
Do
If CheckingNum Mod i = 0 Then GoTo Done 'Composite
i = i + 2
Loop Until i > SqrtCheckingNum
'Prime
Done:


Step 1 of AKS:
If (n = a^b, a > 1, b > 1), output COMPOSITE.

If CheckingNum Mod 2 = 0 Then GoTo Done2 ' Composite
a = 1
Do
a = a + 2
a1 = a * a
If a1 > CheckingNum Then Exit Do
Do
If a1 = CheckingNum Then GoTo Done2 'Composite
a1 = a1 * a
Loop While a1 <= CheckingNum
Loop
'Number is not a perfect square
Done2:


It looks like the bigger numbers I use, the more the AKS code is going to get slower than the naive method. Maybe there's just something I am not seeing that could allow for a radical speed up, but I've been looking at this all day now.
On a side note, it's interesting to see how the graph of the AKS algorithm goes in 'steps', it must be when there's a group of numbers it comes up to that allows it to exit the routine earlier.
For both methods I took the average of two runs and gave the process 'Realtime' priority, so that it is not affected by other processes. (and ran the compiled version instead of in the IDE :) )

Speed for all numbers up to 1 million:

[jemidiah]

Your method for checking if n=a^b, a,b>1, requires in the worst case adding 1+2+2+2... until the sum squared is greater than n. Assuming no work gets done in the inner loop (which isn't at all true) this check takes O(sqrt(n)) operations, which dwarfs the logarithmic runtime of the other steps. It's no surprise the naive method beats this.

It's not immediately obvious to me how to improve this, but your improvement testing n mod 2 getting rid of half the tests strongly suggests large improvements. The AKS paper (http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf) says this step can be accomplished in O(log^3 n) and references "Joachim von zur Gathen and J¨urgen Gerhard. Modern Computer Algebra. Cambridge University Press, 1999" for the method used. In a very brief search I found a discussion (http://www.physicsforums.com/archive/index.php/t-264152.html) of this same problem. They suggest computing successive ith roots of n (eg. sqrt(n), n^(1/3), n^(1/4), ...). You can stop computing when the ith root becomes less than 2, which happens before you've tested log_2 n. There are several further links I'll let you explore. Good luck :).