AKSPrimeQ

Roger Germundsson, Wolfram Research, 8/11/2002

This notebook contains several implementations of a recently proposed deterministic primality test by Agrawal, Kayal and Saxena [web]. The main claim in this article is that with this test they can prove polynomial complexity of primality testing in terms of the number of digits in the input.

Basic Algorithm

This is the the pseudo code version of their algorithm as given in the paper. Using Mathematica pseudo code.

AKSPrimeQ[ n_Integer?Positive ] :=
Catch[
 Module[ {r, q, a},
	If[ "n is of the form (a^b) and b>0", Throw[ False ] ];
	
	r=2;
	While[ r < n,
		If[ GCD[n,r] ≠ 1, Throw[ False ] ];
		If[ PrimeQ[r],
			q = "largest prime factor of r-1";
			If[ q ≥ 4 (Sqrt[ r]) Log[n] ∧ Mod[ (n^(((r - 1)/q))), r ] ≠ 1, Break[] ];
		];
		r++;
	];
	
	For[ a=1, a ≤ 2 (Sqrt[ r]) Log[n], a++,
		If[ (((x - a))^n) ≠  
        (x^n) - a mod (x^r) - 1,  n, Throw[ False ]
	];
	True
 ]
]

As it stands many things need to be cast into an efficient constructive form. But the bulk of the time is spent in the last loop. Imagine for example a typical n of the order (10^30) you realize that you will have compute the (10^30)th power of a polynomial.

This may or may not be a practical algorithm. Most likely it is not, however it can probably be the basis of some efficient algorithms.

Utilities

There are some things we need and which are not performance critical by a long shot. We collect these here.

Implementation

Returns true iff n = (a^b) for b > 1. Note we really don't need the PrimeQ test in the loop, but clearly that is correct and we are essentially betting that it is faster to the PrimeQ check then the arithmetic for the typical ranges of j that will come up.

exactPowerQ1[ 
          n_Integer ]  :=  
Catch[ 
Module[ {nn = 
                N[n, Log[10. , n] + 10]}, 
If[ 
              Round[ (nn^(1/2)) ]^2  ==  n,  
              Throw[ True ] ]; 
Do[ 
              If[ PrimeQ[j]  && Round[(nn^(1/j))]^j ==  n,  
                Throw[ True ] ],  {j,  3,  Floor[Log[2, n]],  
                2}]; 
False
]
];

This is an alternative way of implementing the same idea.

exactPowerQ2[ 
          n_Integer ]  := 
Catch[ 
Module[ {nn = N[n, Log[10. , n] + 10],  
              logn  =  Floor[Log[2, n]],  i = 1,  
              p = Prime[1]}, 
While[ 
              p  ≤  logn, 
If[ 
                Round[ (nn^(1/p))]^p == n,  Throw[ True ], 
                p = Prime[ ++i ] ]
]; 
False
]
];

This implements another subalgorithm that finds the bound r used in the algorithm.

findR[ 
          n_ ]  :=  
Module[ {r,  
            s}, 
Do[
r = 
              s; 
If[ GCD[n, s] ≠ 1,  
              Throw[ False ] ]; 
If[ 
              PrimeQ[ s ], 
q = 
                Max[ FactorInteger[s - 1][[All, 
                    1]] ]; 
If[ 
                q  ≥  4 (Sqrt[s]) Log[n]  ∧  
                  PowerMod[ n, ((s - 1)/q),  s ]  ≠ 1,  
                Break[] ];
], 
{s,  
              2,  n - 
                1}
]; 
r
];

Examples

They generate the same thing.

Table[ exactPowerQ2[n],  {n,  100,  200}] ===  
      Table[ exactPowerQ1[n],  {n,  100,  200}]

They appear to be about the same performance. In all cases the speed of this function appears very small compared to the primality test as a whole.

Timing[ Do[ 
          exactPowerQ2[ n ],  {n,  (10^3),  2 (10^3)}]; ]

Timing[ Do[ 
          exactPowerQ1[ n ],  {n,  (10^3),  2 (10^3)}]; ]

This shows that this performance is more then adequate for our purposes.

Timing[ exactPowerQ1[ 1234747474747474 ] ]

Timing[ exactPowerQ2[ 1234747474747474 ] ]

This is the routine that finds a bound for r used in the algorithm. For the most part not performance critical.

findR[ Prime[10]] // Timing

findR[ Prime[100] ] // Timing

findR[ Prime[1000]] // Timing

Method – Reference Implementation

The basic idea here is to get a simple reference implementation that we can compare other implementations to.

Implementation

This is a slow reference implementation for the test (((x - a))^n) == (x^n) - a mod (x^r) - 1,  
      n.

modularEqualQ0[ a_,  r_,  
          n_ ]  := 
Module[ {x,  pr,  
            pl}, 
pl  = 
            PadRight[
              CoefficientList[ 
                PolynomialRemainder[ (((x - a))^n), (x^r) - 1, x],  x],  
              r]; 
pr  =  
            PadRight[ 
              CoefficientList[ 
                PolynomialRemainder[ (x^n) - a    ,  (x^r) - 1, x],  
                x],  r]; 
Mod[pl, n]  ===  
            Mod[pr, n]
];

AKSPrimeQ0[ 
          n_Integer?
            Positive ]  := 
Catch[ 
Module[ {r,  a}, 
If[ 
              exactPowerQ1[n],  Throw[False] ]; 
r  =  
              findR[ n ]; 

Do[ 
If[ ! modularEqualQ0[ a,  r,  n ],  
                Throw[ False ] ], 
{a,  1,  
                Floor[2 (Sqrt[r]) 
                    Log[n]]}
]; 
True
]
];

Examples

Some simple examples.

AKSPrimeQ0[ Prime[ 10 ] ]

AKSPrimeQ0[ Prime[ 10 ]  +  1 ]

Table[ {i,  AKSPrimeQ0[ i ],  PrimeQ[i]},  {i,  10,  
        20} ]

Performance is of course not that good for this implementation.

AKSPrimeQ0[ Prime[ 30 ] ] // Timing

Method – Closed Form Polynomial Powers

The basic idea here is to get a closed form formulae for the polynomial power used in the test section. Which is to test:

(((x - a))^n)  == (x^n) - a mod (x^r) - 1,  
    n

Another way to see that is to say that we are computing a canonical form in the quotient ring

(Subscript[ℤ, n])[
        x]/⟨(x^r) - 
          1⟩

In this ring we have:

(x^k) mod (x^r) - 
      1 ⟹ (x^Mod)[k, r]

And hence by using the binomial theorem and the periodicity condition above:

(((x - a))^n) mod (x^r) - 1
==(∑_(i = 0)^n) Binomial[n, 
            i] (x^i) (((-a))^(n - i))
\(\(\(==\)\(\ \)\(∑\_\(i = 0\)\%\(r - 1\)\((∑\_\(j = 0\)\%\(q - 1\)Binomial[n, j\ r\  + \ i]\ \((\(-a\))\)\^\(n - j\ r - i\))\) 
            x\^i\  + ∑\_\(i = 0\)\%\(n - q\ r + 1\)Binomial[n, 
              q\ r\  + \ 
                i] \((\(-a\))\)\^\(n - q\ r - i\)\)\)\)

(x^n) - a mod (x^r) - 1 == (x^Mod)[n, r] - a

This illustrates these formulae.

Block[ {n = 10,  q = 3,  
        r = 3}, 
Collect[ (∑_(i = 0)^(r - 1))((∑_(j = 0)^(q - 1))Binomial[n, 
                      j r  +  i] (((-a))^(n - j r - i))) 
              (x^i)  +  (∑_(i = 0)^(n - q r + 1))Binomial[n, 
                q r  +  i] (((-a))^(n - q r - i)) (x^i),  
        x ]
]

PolynomialRemainder[ (((x - a))^10),  (x^3) - 1,  x ]

This idea is implemented below.

Implementation 1 (List Implementation)

In this case we generate coefficient lists of length r.

Checking the polynomial identity based on closed form formulae.

modularEqualQ1[ a_,  r_,  
          n_ ]  := 
Module[ {q = Quotient[n, r],  pl,  
            pr}, 
pl  =  
            Table[ 
Mod[ (∑_(j = 0)^(q - 1))Mod[ Binomial[n, j r  +  i] PowerMod[-a,  n - j r  - i,  n],  n],  
                  n ]  + 
If[ 
                  q r  +  i ≤ n,  
                  Mod[ Binomial[n, q r  +  i ] PowerMod[-a,  
                        n - q r  - i,  n],  n],  
                  0], 
{i,  0,  
                r - 1} 
]; 
pr  =  
            Table[ KroneckerDelta[ Mod[n, r] - i] - 
                a KroneckerDelta[ i ],  {i,  0,  
                r - 1} ]; 
Mod[pl,  n]  ===  
            Mod[pr,  n]
];

AKSPrimeQ1[ 
          n_Integer?
            Positive ]  := 
Catch[ 
Module[ {r,  a}, 
If[ 
              exactPowerQ1[n],  Throw[False] ]; 
r  =  
              findR[ n ]; 

Do[ 
If[ ! modularEqualQ1[ a,  r,  n ],  
                Throw[ False ] ], 
{a,  1,  
                Floor[2 (Sqrt[r]) 
                    Log[n]]}
]; 
True
]
];

Examples 1 (List Implementation)

Some simple examples.

AKSPrimeQ1[ Prime[ 10 ] ]

AKSPrimeQ1[ Prime[ 10 ]  +  1 ]

Table[ {i,  AKSPrimeQ1[ i ],  PrimeQ[i]},  {i,  10,  
        20} ]

This is actually slower then doing the polynomial implementation.

AKSPrimeQ1[ Prime[ 30 ] ] // Timing

Implementation 2 (No List Implementation)

We can improve on the previous idea, by realizing that we don't need to maintain the full list of coefficients. But rather just compare coefficient by coefficient and exit the first time this is no longer the case.

This implements the idea of comparing coefficient by coefficient and exit as soon as there is a discrepancy.

modularEqualQ2[ a_,  r_,  
          n_ ]  := 
Catch[
Module[{q = 
                Quotient[n, r],  
              t}, 
Do[ 
(
t  =  (Mod[ 
Sum[ 
                            Mod[ Binomial[n, j r + i] PowerMod[-a,  
                                  n - j r - i,  n],  n ],  {j,  0,  
                              q - 1}] + 
                          If[q r + i ≤ n,  
                            Mod[Binomial[n,  q r  +  i] 
                                PowerMod[-a,  n - q r - i,  n],  
                              n ],  0 ],  
                        n ]  ===  
Which[ 
                        i == 0,  
                        Mod[ -a,  n ], 
i ==  
                          Mod[n, r],  1,  
True,  
                        0 ]); 
If[ t  ===  False,  
                  Throw[ 
                    False ] ]
),  
{i,  0,  
                r - 1}
]; 
t
]
];

AKSPrimeQ2[ 
          n_Integer?
            Positive ]  := 
Catch[ 
Module[ {r,  a}, 
If[ 
              exactPowerQ1[n],  Throw[False] ]; 
r  =  
              findR[ n ]; 

Do[ 
If[ ! modularEqualQ2[ a,  r,  n ],  
                Throw[ False ] ], 
{a,  1,  
                Floor[2 (Sqrt[r]) 
                    Log[n]]}
]; 
True
]
];

Examples 2 (No List Implementation)

Some simple examples.

AKSPrimeQ2[ Prime[ 10 ] ]

AKSPrimeQ2[ Prime[ 10 ]  +  1 ]

Table[ {i,  AKSPrimeQ2[ i ],  PrimeQ[i]},  {i,  10,  
        20} ]

This represents only marginal time saving in this case.

AKSPrimeQ2[ Prime[ 30 ] ] // Timing

Method – Fast Polynomial Powers

The basic idea here is to use polynomial arithmetic, but to make use of the special multiplications we are doing and the ring that we are doing them in.

The first idea is to use that the fact that multiplication by a fixed polynomial x - a in the ring (Subscript[ℤ, n])[
          x]/⟨(x^r) - 1⟩ is really a linear operation, which in our case there are simple list operations in Mathematica for. For suppose we have:

(((x - a))^k) == (Subscript[α, 0])[k] + (Subscript[α, 1])[k] 
          x + ⋯ + (Subscript[α, r - 1])[k] 
          (x^(r - 1))

then we have:

(((x - a))^(k + 1)) == (x - a) (((x - a))^k) ==  
      x (((x - a))^k) - a (((x - a))^k)

Hence it is enough know how x (((x - a))^k) and -a (((x - a))^k) changes the coefficient vector α[k]. Here we have:

x (((x - a))^(k - 1)) ⟺ β[
            k]== (({{0, 0, ⋯, 1}, {1, ⋱, 0, 0}, {⋮, ⋱,  , 0}, {0, 0, 1, 0}}))α[k - 1]==RotateRight[α[k - 1]]

\(\(-a\)\ \((x - a)\)\^\(k - 1\)⟺γ[k]\)== (({{-a, 0, ⋯, 0}, {0, -a,  , 0}, {⋮, ⋱, ⋱, 0}, {0, 0, 0, -a}}))α[k - 1]== -a α[k - 1]

And hence we have:

This idea is used in implemenation 1 below.

Another related idea is we could implement a fast power method, i.e. a power tree method by using repeated squaring. This means that we typically only need to perform Log[2, n] multiplications, but we will need a general multiplication and not only the special multiplication above. In fact we can do this very neatly by using ListConvolve with appropriate cyclic padding behavior.

Implementation 1 – Special Polynomial Multiplication

This is a reference implementation using polynomial operations.

iteratedProduct0[ a_,  r_,  n_ ]  :=  
Mod[ 
          CoefficientList[ 
            PolynomialRemainder[ (x - a)^n,  x^r - 1,  x ],  x ],  
          n ];

This implements the linear system above using list rotation and list multiplication.

iteratedProduct1[ a_,  r_,  
          n_ ]  := 
Module[{α0}, 
α0  = 
            Table[ Which[ i == 0,  Mod[-a,  n],  
                i == 1,  1,  True,  0],  {i,  0,  
                r - 1}];  
Nest[ 
            Mod[ RotateRight[#] -  a #,  n ]  &,  α0,  
            n - 1 ]
];

This in turn gives another test for the primality testing algorithm above ... implemented below.

modularEqualQ3[ a_,  r_,  
          n_ ]  := 
Module[{pl,  
            pr}, 
pl  =  
            iteratedProduct1[ a,  r,  n ]; 
pr  =  
            Table[ KroneckerDelta[ Mod[n, r] - i] + 
                Mod[-a,  n] KroneckerDelta[ i ],  {i,  0,  
                r - 1} ]; 
pl  ===  
            pr
];

AKSPrimeQ3[ 
          n_Integer?
            Positive ]  := 
Catch[ 
Module[ {r,  a}, 
If[ 
              exactPowerQ1[n],  Throw[False] ]; 
r  =  
              findR[ n ]; 

Do[ 
If[ ! modularEqualQ3[ a,  r,  n ],  
                Throw[ False ] ], 
{a,  1,  
                Floor[2 (Sqrt[r]) 
                    Log[n]]}
]; 
True
]
];

Examples 1 - Special Polynomial Multiplication

Some simple examples.

AKSPrimeQ3[ Prime[ 10 ] ]

AKSPrimeQ3[ Prime[ 10 ]  +  1 ]

Table[ {i,  AKSPrimeQ3[ i ],  PrimeQ[i]},  {i,  10,  
        20} ]

This represents more time savings ...

AKSPrimeQ3[ Prime[ 30 ] ] // Timing

Implementation 2 – Fast and Special Polynomial Multiplication

Polynomial multiplication modulo (x^r) - 1and modulo n. The list lenghts are assumed to be r. This correspond to cyclic convolution.

PolynomialTimes[ al_List,  bl_List,  
          n_Integer ]  :=  
Mod[ 
          ListConvolve[al,  bl,  {1, 1}],  n ];

This implements multiplication by x - a modulo (x^r) - 1.

PolynomialTimes1[ al_List,  a_Integer,  
          n_Integer ]  :=  
Mod[ 
          RotateRight[al] -  a al,  n ];

This implements the power tree method for this polynomial to the (k^th) power. It basically uses the Horner form of the binary expansion of the exponent. It computes the powers from the inside out.

FastPolynomialPower[ a_,  r_,  
          n_ ]  := 
Module[ {nbits,  
            pow}, 
nbits  =  
            Rest[ IntegerDigits[ n,  2 ] ]; 
pow  =  
            Table[ Which[ i == 0,  Mod[-a, n],  i == 1,  
                1,  True,  0],  {i,  0,  
                r - 1}]; 
While[ 
            Length[ nbits ]  >  0, 
pow  =  
              PolynomialTimes[ pow,  pow,  n ]; 
If[ 
              First[nbits]  ===  1,  
              pow  =  
                PolynomialTimes1[ pow,  a,  
                  n ] ]; 
nbits  =  
              Rest[ nbits ];
]; 
pow
];

This gives yet another version of the equation check.

modularEqualQ4[ a_,  r_,  
          n_ ]  := 
Module[{pl,  
            pr}, 
pl  =  
            FastPolynomialPower[ a,  r,  n ]; 
pr  =  
            Table[ KroneckerDelta[ Mod[n, r] - i] + 
                Mod[-a,  n] KroneckerDelta[ i ],  {i,  0,  
                r - 1} ]; 
pl  ===  
            pr
];

AKSPrimeQ4[ 
          n_Integer?
            Positive ]  := 
Catch[ 
Module[ {r,  a}, 
If[ 
              exactPowerQ1[n],  Throw[False] ]; 
r  =  
              findR[ n ]; 

Do[ 
If[ ! modularEqualQ4[ a,  r,  n ],  
                Throw[ False ] ], 
{a,  1,  
                Floor[2 (Sqrt[r]) 
                    Log[n]]}
]; 
True
]
];

Examples 2 – Fast and Special Polynomial Multiplication

Some simple examples.

AKSPrimeQ4[ Prime[ 10 ] ]

AKSPrimeQ4[ Prime[ 10 ]  +  1 ]

Table[ {i,  AKSPrimeQ4[ i ],  PrimeQ[i]},  {i,  10,  
        20} ]

This represents only marginal time saving in this case.

AKSPrimeQ4[ Prime[ 30 ] ] // Timing

The difference is more noticeable for large examples.

AKSPrimeQ4[ Prime[ 50 ] ] // Timing

AKSPrimeQ3[ Prime[ 50 ] ] // Timing

AKSPrimeQ2[ Prime[ 50 ] ] // Timing

AKSPrimeQ1[ Prime[ 50 ] ] // Timing

Slightly larger ...

AKSPrimeQ4[ Prime[ 100 ] ] // Timing

AKSPrimeQ3[ Prime[ 100 ] ] // Timing

AKSPrimeQ2[ Prime[ 100 ] ] // Timing

AKSPrimeQ1[ Prime[ 100 ] ] // Timing

Method – DanL and DTerr

The implementation below uses finite field polynomial operations available only in Mathematica 5.0. The implementation was done by David Terr and Dan Lichtblau.

Only included for reference and reference timings ...

Implementation

Note: This requires Mathematica 5.0 ...

$Version

isPowerQ[n_] := 
      Module[{logn = Floor[Log[2, n]], b, nn = N[n, Log[10. , n] + 10]}, 
        If[Round[(Sqrt[nn])]^2 == n, Return[True]]; 
        For[j = 3, j < logn, 
          j += 2, If[PrimeQ[j], b = Round[(nn^(1/j))]; 
              If[(b^j) == n, Return[True]]];]; False]

\(MyPrimeQ[
        n_] := 
Module[{r, q, a, done = False, b, logn, f, 
          g1, g}, 
If[isPowerQ[n], 
          Return[False]]; 
logn = 
          4. *Log[2. , n]; 
r = 
          Ceiling[4*logn^2]; 
If[\(! OddQ[
              r]\), \(r++\)]; 
Print["init r ", 
          r]; 
While[\(! done\), 
If[
            GCD[n, r] ≠ 1, Print[r]; 
            Return[False]]; 
If[
            PrimeQ[r], 
q = 
              Max[\(FactorInteger[\((r - 1)\)/2]\)[\([All, 
                    1]\)]]; 
If[
              q > Sqrt[N[r]]*logn && 
                PowerMod[n, \((r - 1)\)/q, r] ≠ 
                  1, 
\(done = 
                  True;\)
];
]; 
r\  += 2;
]; 
b = 
          Sqrt[N[r]]*logn/2. ; 
Print["r ", 
          r]; 
Print["b ", b]; 
h = 
          Join[{n - 1}, Table[0, {r - 1}], {1}]; 
g1 = 
          Algebra`PolynomialPowerModList[{0, 1}, n, h, 
            n]; 
For[a = 1, 
          a < b, \(a++\), 
f = 
            Algebra`PolynomialPowerModList[{n - a, 1}, n, h, 
              n]; 
g = 
            Algebra`PolynomialPlusModList[g1, {n - a}, 
              n]; 
If[f =!= g, 
            Return[False]];
]; 
Return[
          True];
]\)

Examples

$Version

Small example.

MyPrimeQ[ Prime[10] ] // Timing

AKSPrimeQ4[ Prime[ 10 ] ] // Timing

AKSPrimeQ3[ Prime[ 10 ] ] // Timing

AKSPrimeQ2[ Prime[ 10 ] ] // Timing

AKSPrimeQ1[ Prime[ 10 ] ] // Timing

Slightly larger example.

MyPrimeQ[ Prime[50] ] // Timing

AKSPrimeQ4[ Prime[ 50 ] ] // Timing

AKSPrimeQ3[ Prime[ 50 ] ] // Timing

AKSPrimeQ2[ Prime[ 50 ] ] // Timing

AKSPrimeQ1[ Prime[ 50 ] ] // Timing