#########################################################
# AKS, A Deterministic Primality Test
#########################################################
# Purpose: Decide in polynomial time whether
# a given positive integer n is prime or not.
# Method: "PRIMES is in P" (M. Agrawal, N. Kayal,
# N. Saxena). Note that we also use improvements due
# to D. Bernstein and F. Morain.
# Author: Peter Luschny, 2002/09/27, Ver. 3.1
# Implementation: Maple V, Rel. 5, Ver. 5.0
# Remerciements: For discussion, hints and bug hunting.
# Gottfried Barthel, Alfred Heiligenbrunner and Carl Devore.
# Updates: www.luschny.de/math/primes/aks.txt
#########################################################
# An Outline of the AKS-Algorithm.
#########################################################
# Pseudocode:
#
# aks_isprime := proc(n)
# if(n<2) then return('ERROR') fi;
# if(n<4) then return('PRIME') fi;
# if((n mod 2) = 0) then return('COMPOSITE') fi;
# if(isSquare(n)) then return('COMPOSITE') fi;
# r := 2;
# loop
# r := nextprime(r);
# if(r divides n) then return('COMPOSITE') fi;
# if(r >= sqrt(n)) then return('PRIME') fi;
# q := the largest prime factor of (r-1) ;
# maxa := 2*sqrt(r)*log[2](n);
# if((q >= 2*maxa) and (n^((r-1) div q) <> 1 (mod r)))
# then exit_loop fi;
# end_loop;
# for a from 1 to maxa
# if( (x-a)^n <> (x^n-a) (mod x^r-1,n) )
# then return('COMPOSITE') fi;
# end_for;
# for i from 3 to log[2](n)
# if(i_throot(n)^i = n) then return('COMPOSITE') fi;
# end_for;
# return('PRIME');
# end_proc.
#
#=========================================================
#
# Four variants of the AKS-algoritm are implemented
# within a common framework:
# 'isprime_AKS' as described by Agrawal/Kayal/Saxena,
# 'isprime_AKS_B' the variant given by D. Bernstein
# 'isprime_AKS_SG' assuming conjecture 5 in the AKS paper
# 'isprime_AKS_BP' assuming conjecture 6 in the AKS paper
# You will see, that the Bernstein variant is more efficient
# than the original one. And that, if the conjectures are true,
# they would lead to much faster solutions which might even gain
# relevance in practice.
# We outline the Bernstein variant of the algorithm.
# The algorithm has three parts, which correspond here to the
# three functions "powerTest", "polyTest" and "isPerfectPower".
# Part 1:
# "powerTest(n)" returns .
# if d = -1, then n is prime,
# because no k <= sqrt(n) divides n.
# id d > 0, then d is a divisor of n.
# if d = 0, we do not know yet, whether n
# is prime or not - but two positive
# integers r and s were calculated, such that
#
# (1) binomial(q+s-1,s) > n^(2floor(sqrt(r)))
# with q the greatest prime divisor of r-1,
# (2) n has no prime factors <= s and
# (3) n^((r-1)/q) mod r is not in {0,1}.
# Part 2:
# "polyTest(n,r,s)" checks the following relation:
# (x-a)^n = x^n-a (mod x^r - 1,n) for all 1<=a<=s.
# If not all of these relations are true, n is not prime,
# otherwise n is the perfect power of some prime n=p^k.
# Note that this function is the function which consumes
# most of the computig time.
# Part 3:
# If the causa is not yet finita, n = p^k, p prime, and
# we have to check k=1 or k>1. This is done with the
# function "isPerfectPower(n)", which returns k.
# Caveat: This primality test is slow!
# For example: To prove n=10000019 prime (variant
# Bernstein), it takes 6 and 1/2 hours on my machine.
# The time is consumed in the function 'polyTest',
# where the for-loop is performed 1785 times and
# each call to Powmod consumes about 13 seconds.
#
#########################################################
# USAGE
#########################################################
# Check n=10007 and use the original AKS method.
# isprime_AKS(10007);
# Check n=10007 and use the AKS_Bernstein method.
# isprime_AKS_B(10007);
# Similarly isprime_AKS_SG(10007) and isprime_AKS_BP(10007);
#
# (verbose := false and verboseOnPrime := false) <=> quiet
# (verbose := false and verboseOnPrime := true) <=> some infos
# (verbose := true and verboseOnPrime := true) <=> all infos
# timing := true: <=> shows time used in subtests and total time
#########################################################
# Implementation
#########################################################
readlib(powmod): readlib(iroot): readlib(issqr):
floorRoot := proc(n) local r;
#----------------------------
r := isqrt(n);
if (r*r > n) then r := r-1 fi;
r end:
#maxPrimeFactor := proc(n) local L;
#---------------------------------
#L := op(2,numtheory[ifactors](n)):
#op(1,op(nops(L),L)) end:
# Carl Devore suggested the following
# improved coding:
maxPrimeFactor := n -> numtheory[ifactors](n)[2,-1,1]:
#---------------------------------
# Called "for-loop" in the paper
polyTest := proc(n,r,s) local a,rhs,lhs;
#---------------------------------------
#rhs := Powmod(x^n-a,1,x^r-1,x) mod n;
rhs := (Power(x,n mod r)-1) mod n;
for a from 1 to s do
lhs := Powmod(x-a,n,x^r-1,x) mod n;
if(lhs <> rhs) then RETURN(false) fi;
rhs := rhs - 1;
od;
RETURN(true) end:
# AKS original, the "while-loop"
powerTestA := proc(n)
local r,s,q,d,log2n,sqrtn,sqrtr,loop;
#------------------------------------
log2n := evalf(log[2](n));
sqrtn := floorRoot(n);
d := 0; s := 0; r := 2;
loop := true;
while(loop) do
r := nextprime(r);
if(irem(n,r) = 0) then d := r; break; fi;
if(r >= sqrtn) then d := -1; break fi;
q := maxPrimeFactor(r-1);
sqrtr := floorRoot(r);
s := floor(evalf(2*sqrtr*log2n));
if(q >= 2*s) then
if ((Power(n,iquo(r-1,q)) mod r) <> 1)
then loop := false;
fi
fi
od;
if(verbose) then
if(d = 0) then print(`n=`,n,``,r,q,s); fi;
fi;
(r,s,d) end:
# AKS-Bernstein, the "while-loop"
powerTestB := proc(n)
local r,s,q,d,m,c,cMin,cMax,logn,sqrtn,sqrtr,loop;
#-------------------------------------------------
# Note that I calculate cMin and cMax and start
# searching s only if cMin <= cMax, which means
# only if the existence of c(s) is assured.
logn := evalf(ln(n));
sqrtn := floorRoot(n);
d := 0; s := 0; r := 2;
loop := true;
while(loop) do
r := nextprime(r);
if(irem(n,r) = 0) then d := r; break fi;
if(r >= sqrtn) then d := -1; break fi;
q := maxPrimeFactor(r-1);
if((Power(n,iquo(r-1,q)) mod r) > 1) then
sqrtr := floorRoot(r);
cMin := evalf(2*sqrtr*logn);
cMax := evalf(1/8*(-ln(Pi)+ln(2)*(4*q-1)
-ln(2*q-1))*(8*q-5)/(2*q-1));
if(cMin <= cMax) then
c := cMax;
for s from q-1 by -1 to 1 do
c := evalf(c - ln(q+s-1) + ln(s));
if(is(cMin > c))
then loop := false; s := s+1; break fi;
od;
fi;
fi;
od;
if(verbose) then
if(s > 0) then
c := evalf(lnGAMMA(q+s-1)-lnGAMMA(s)-lnGAMMA(q-1));
print(`n=`,n,``,r,q,s);
print(``,cMin,c,cMax);
fi;
fi;
(r,s,d) end:
# Some values for checking Bernstein's powerTestB
# n
# 10000019 -> 2207, 1103, 1045; 1482.8, 1489.6, 1524.1
# 100000007 -> 2879, 1439, 1387; 1952.5, 1960.0, 1989.8
# 1000000007 -> 3623, 1811, 1785; 2486.7, 2494.2, 2505.3
# 10000000019 -> 4547, 2273, 2188; 3085.4, 3093.2, 3145.7
# 100000000003 -> 5387, 2693, 2651; 3697.9, 3706.0, 3727.9
# 1000000000039 -> 6599, 3299, 3169; 4476.2, 4484.0, 4567.8
# 10000000000037 -> 7523, 3761, 3676; 5148.5, 5156.6, 5208.2
# 100000000000031 -> 8699, 4349, 4310; 5995.9, 6004.1, 6023.3
# AKS with conjecture 5, Sophie Germain, the "while-loop"
powerTestSG := proc(n)
local r,s,q,d,log2n,log2nSqr64,sqrtn,sqrtr,loop;
#----------------------------------------------
log2n := evalf(log[2](n));
log2nSqr64 := floor(log2n*log2n*64.0);
sqrtn := floorRoot(n);
d := 0; s := 1; r := 2;
loop := true;
while(loop) do
r := nextprime(r);
if(irem(n,r) = 0) then d := r; break; fi;
if(r >= sqrtn) then d := -1; break fi;
# Just immediately reject any r such that r-1 does
# not have a large enough maximal prime factor.
# Just check isprime((r-1)/2).
# It hasn't been proven that this will work -- that there
# are sufficiently many co-Sophie Germain primes.
# [Carl Devore on comp.soft-sys.math.maple]
#
# Note, that in this case the condition
# q = (r-1)/2 >= 4sqrt(r)log[2](n)
# implies r >= 64(log[2](n))^2.
if(r >= log2nSqr64) then
q := (r-1)/2;
if(isprime(q)) then
if ((Power(n,iquo(r-1,q)) mod r) <> 1)
then loop := false;
fi
fi
fi
od;
if(verbose) then
if(d = 0) then print(`n=`,n,``,r,q,s); fi;
fi;
(r,s,d) end:
# AKS with conjecture 6, the "while-loop"
powerTestBP := proc(n) local r,s,d,n2,sqrtn;
#-------------------------------------------
n2 := n*n-1; sqrtn := floorRoot(n);
d := 0; s := 1; r := 2;
do r := nextprime(r);
if(irem(n,r) = 0) then d := r; break fi;
if(r >= sqrtn) then d := -1; break fi;
if(irem(n2,r) <> 0) then break fi;
od;
if(verbose) then
if(d = 0) then print(`n=`,n,`r=`,r); fi;
fi;
(r,s,d) end:
isPerfectPower := proc(n) local i,log2n;
# returns i iff n = m^i, i>1, 1 otherwise
#----------------------------------------
log2n := floor(evalf(log[2](n)));
for i from 2 to log2n do
if(iroot(n,i)^i = n) then RETURN(i) fi;
od;
RETURN(1) end:
AKS := proc(n,algo) local r,s,d,tm,tp,p;
# option algo = {ORIG, BERN, CONSG, CONBP }
#-------------------------------------------
# First some trivial cases out of the way.
if(n<2) then
if(verbose) then print(`Input not valid`,n) fi;
RETURN(false)
fi;
if(n<4) then
if(verboseOnPrime) then print(n,`is prime`) fi;
RETURN(true)
fi;
if(type(n,even)) then
if(verbose) then print(`2 | `,n) fi;
RETURN(false)
fi;
if(issqr(n)) then
if(verbose) then print(isqrt(n),`| `,n) fi;
RETURN(false)
fi;
# The first test, the powerTest, will find most
# small factors by trial division or prove n prime
# by finding no factors <= sqrt(n).
# Otherwise it returns values r and s for the next test.
if(verbose) then lprint(`Using algorithm `, algo) fi;
tm := time():
if(algo=ORIG ) then (r,s,d) := powerTestA(n);
else if(algo=BERN ) then (r,s,d) := powerTestB(n);
else if(algo=CONSG) then (r,s,d) := powerTestSG(n);
else if(algo=CONBP) then (r,s,d) := powerTestBP(n);
fi fi fi fi;
if(timing) then
print(`Time used in powerTest`,time()-tm,`sec.`);
fi;
if(d = -1) then
if(verboseOnPrime) then
print(n,`is prime (trial division)`)
fi;
RETURN(true);
fi;
if(d > 0) then
if(verbose) then print(d,` | `,n) fi;
RETURN(false);
fi;
# We had no success and have to perform a second test,
# the polyTest. Things are not easy any more.
# We have to work with polynomials in Z/nZ[x]/(x^r-1)
# and perform s tests.
# While waiting for the result, better remember that
# your life time is bounded by polynomial time too ;-)
# Gottfried Barthel calls it the "black hole" of the
# algorithm :-) He claims that n=3640471 is the first
# number, which passes this line (Bernstein variant).
tp := time():
p := polyTest(n,r,s);
if(timing) then
lprint(`Time used in polyTest`,time()-tp,`sec.`);
lprint(`Total time used `,tm + tp,`sec.`)
fi;
if(p) then
if(algo <> CONBP) then
# Now n = p^r, p prime, r >= 1
r := isPerfectPower(n);
if(r > 1) then
if(verbose) then
print(n,`is perfect power, exponent =`,r)
fi;
RETURN(false);
fi;
fi;
if(verboseOnPrime) then print(n,`is prime (poly test)`) fi;
RETURN(true);
else
if(verbose) then
print(n,`is composite (poly test)`);
fi;
RETURN(false)
fi end:
##############################################################
# The public interface
##############################################################
isprime_AKS := proc(n) AKS(n,ORIG) end: # Original
isprime_AKS_B := proc(n) AKS(n,BERN) end: # Bernstein
isprime_AKS_SG := proc(n) AKS(n,CONSG) end: # Conjecture 5
isprime_AKS_BP := proc(n) AKS(n,CONBP) end: # Conjecture 6
verbose := false: # gives some additional informations (or not)
verboseOnPrime := true: # verbose mode only if n is prime
timing := false: # shows time used in subtests and total time
##############################################################
# Some Test Functions
##############################################################
testSmall := proc(n, algo) local i;
for i from 1 by 2 to n do
if(is(algo(i) <> isprime(i)))
then print(i,algo,"ERRORRR!!!");
fi;
od: ` ` end:
testSmall(100,isprime_AKS);
testSmall(100,isprime_AKS_B);
testSmall(100,isprime_AKS_SG);
testSmall(100,isprime_AKS_BP);
testBig := proc(n, algo) local i;
for i from 3636601 by 2 to (3636600+n) do
algo(i);
od: ` ` end:
testBig(100,isprime_AKS_B);
testBig(100,isprime_AKS_SG);
testBig(100,isprime_AKS_BP);
testPower := proc(n,algo) local i,p,t;
for i from 1 to n do
t := time(): p := algo(nextprime(10^i)); t := time()-t;
print(`----> Total time`,t,`sec.`);
od: end:
testPower(6,isprime_AKS_B);
testPower(6,isprime_AKS_SG);
testPower(6,isprime_AKS_BP);
#########################################################
# A AKS-FAQ for the beginner
#########################################################
#
# The conjecture that there is a polynomial time primality
# test has been confirmed by Agrawal, Kayal and Saxena
# (Indian Institute of Technology in Kampur). The paper
# has a date of August 6, 2002. But was distributed to
# several number theorists earlier. Hendrik Lenstra just
# wrote to the Number Theory mailing list saying that the
# proof is correct, clever and elegant; and it is elementary
# except for one result from analytic number theory
# needed to establish the running time.
#
# Edwin Clark on sci.math, 2002-08-06 18:35:32 PST
#========================================================
# Does this imply P=NP ?
# In other words: is it proven that primality test
# is NP-hard?
# ---
# No. Before now, PRIMES was known to be in (NP intersect
# coNP) which is a superset of P. PRIMES has not been proven
# complete for the class (NP intersect coNP), though.
#
# P=NP would clearly imply that (NP intersect coNP)=P, but it
# is also possible that P != NP but (NP intersect coNP)=P.
#
# Stephen Forrest on sci.math
#
# Note that if P = (NP intersect coNP) then
# factoring can be done in polynomial time.
#
#========================================================
#
# Lance Fortnow, Friday, September 20, 2002, writes:
# http://fortnow.com/lance/complog/
#
# I saw Carl Pomerance yesterday give a wonderful
# presentation on the AKS primality algorithm. He
# made an interesting point about the algorithm.
# The algorithm runs in time O(n^12) where n is the
# number of bits of the input to be tested. The big O
# notation hides a constant, i.e., the algorithm uses
# c n^12 for some constant c. That c is unknown!
#
# The AKS algorithm uses a result by Fouvry on the
# distribution of certain kinds of primes. Fouvry's
# result uses another result that is proven as such:
# First it is proved assuming the Extended Riemann
# Hypothesis is true. If the ERH fails, then the place
# where it fails can be used to prove the result.
# The constant c will depend on where the ERH fails.
# To determine c would require settling the Extended
# Riemann Hypothesis!
#========================================================
#
# How did the primality tests work until now?
#
# The best primality tests today use randomization and
# therefore could in principle give the wrong answer,
# or else could only be proven to run in polynomial-time
# if we assumed the Extended Riemann Hypothesis.
# In contrast, the AKS test for primality, does not involve
# randomization, and can be proven to run in polynomial time,
# without having to assume anything like the ERH.
#========================================================
#
# The paper "PRIMES is in P" does contain a negated
# congruence (x-a)^n =/= x^n - a (mod x^r-1, n).
# What this means is that (x-a)^n and x^n - a
# (which are elements of the polynomial ring Z[x])
# are not congruent modulo the ideal of Z[x] generated
# by x^r - 1 and n.
#
# In general two elements r and s of a ring R are congruent
# modulo the ideal I of R iff r - s is an element of R.
# The AKS statement says that (x-a)^n - (x^n - a) is not
# an element of the ideal (x^r - 1, n) of Z[x].
# This ideal is the set of polynomials (x^r-1)f(x) + n g(x)
# where f and g lie in Z[x].
#
# So the AKS statement says that (x-a)^n - (x^n - a)
# cannot be expressed as (x^r-1)f(x) + n g(x) for
# polynomials f(x) and g(x) with integer coefficients
#
# Robin Chapman in sci.math
#
# Remark: With Maple you might play with the following
# function as a warming up.
#
# checkCongruence := proc(n,r) local lh,rh,mlh,mrh;
# rh := powmod(x^n-1,1,x^r-1,x): mrh := rh mod n;
# lh := powmod(x-1,n,x^r-1,x): mlh := lh mod n;
# lprint(n, r, mlh - mrh);
# mlh = mrh end:
#========================================================
#
# Given a poly P(X) with integer coef's,
# one may compute P(X) mod (n, X^r - 1) simply
# by reducing each of monomial of P as follows
#
# e (e mod r)
# c X -> (c mod n) X
#
# In other words perform polynomial arithmetic
# modulo the equations n = 0 and X^r = 1,
# i.e. compute modulo the ideal (n, X^r - 1).
#
# Bill Dubuque on sci.math
#
# Remark: This can be used to simplify the 'polytest'.
# rhs = Powmod(x^n-1,1,x^r-1,x) mod n;
# = (Power(x,n mod r)-1) mod n;
#========================================================
#
# I just scanned the paper and it says its time complexity
# is O((log n)^12), but there is a bar over the O which
# I'm not familliar with.
# ---
# The bar above the O can be taken as a +o(1) in the
# exponent, as log(D) to any constant power will be
# dominated by D.
#
# Phil Carmody on sci.math
#
# Remark: In PostScript parlance 'the O with the bar
# above' is called 'Otilde'.
#========================================================
#
# For those curious about the practicality of the
# new algorithm:
#
# Using n = 628363443011, a very reasonably sized prime,
# unless there's a mistake in my computation, the
# smallest r that works is 97463. That's with rounding
# down on the logs and sqrts. Rounding up, r = 103043.
# Computing the log and the sqrt in floating point,
# multiplying, and rounding normally, r = 98387.
#
# So even finding the r is quite time consuming.
# But after finding the r, you still have to compute
# 2*sqrt(r)*lg(n) = 24471 polynomials, each of which has
# roughly 97463 terms.
#
# To be fair, the algorithm only guarantees a relatively
# small r when n is very large. But I doubt if the
# polynomial computations can ever be practical.
#
# Carl Devore on sci.math.symbolic (edited by P.L.)
# (Note: Carl speaks about the *original* AKS algorithm.)
#========================================================
#
# My question is: Which is the smallest composite
# number, which survives the first test (i.e. the
# "while-loop", the "powerTest")
#
# (a) in the AKS-case
# (b) in the Bernstein-case?
# -------
#
# In the AKS-case:
#
# p1 p2 n=p1*p2 r
# -------------------------------
# 65551 65557 4297326907 65543
#
# Aldo D. in sci.crypt.
#
# A *preliminary* search showed for the Bernstein
# case a much smaller n
#
# n := 3694003 = 1913*1931;
#
# r := 1907; irem(n,r) = 144;
# ifactor(r-1) = { 2, 953 };
# q := 953; s := 931;
# n^((r-1)/q) mod r = 1666;
# log[2](binomial(q+s-1,s)) = 1877.066;
# 2*floor(sqrt(r))*log[2](n) = 1876.241;
#
# Peter L. in sci.crypt.
#
# Remark:
# The significance of these constants is, that they allow
# an 'early out' in the algorithm for smaller values of n.
# (Not implemented.)
#========================================================
#
# What is 'conjecture 5' (SG)?
# ---
#
# * The number of co-Sophie Germain primes is asymptotic
# * to Cx/(log(x)^2), where C is the twin prime constant.
# * C is estimated to be approximately 0.6601618..
#
# This conjecture has been verified for r <= 10^10.
#
# A prime is called co-Sophie-Germain prime
# if both r and (r-1)/2 are primes.
#
# Remark:
# In Maple these primes are called 'safe' primes.
# To see the first few, execute
# [seq(numtheory[safeprime](i),i=1..2998)]:
# sort(convert(convert(%,set),list));
#========================================================
#
# What is 'conjecture 6' (BP)?
# (Bhattacharjee and Pandey, 2001)
# ---
#
# * (1) if r does not divide n, and
# * (2) if r does not divide n^2-1, and
# * (3) if (x-1)^n = x^n-1 (mod x^r - 1,n) holds,
# * then n is prime!
#
# This conjecture has been verified for n <= 10^10 and r <= 100.
#========================================================
#
# Chris Caldwell on "The Prime Pages"
# http://www.utm.edu/research/primes/prove/prove4_3.html
#
# The key to AKS' result is another simple version
# of Fermat's Little Theorem:
#
# Theorem: Suppose that a and p are relatively prime
# integers with p > 1.
#
# p is prime if and only if (x-a)^p = (x^p-a) (mod p)
#
# Proof. If p is prime, then p divides the binomial
# coefficients binomial(p,r) for r = 1, 2, ... p-1.
# This shows that (x-1)^p = (x^p-a^p) (mod p), and the
# equation above follows via Fermat's Little Theorem.
# On the other hand, if p > 1 is composite, then it
# has a prime divisor q. Let q^k be the greatest power
# of q that divides p. Then q^k does not divide
# binomial(p,q) and is relatively prime to a^(p-q), so the
# coefficient of the term x^q on the left of the equation
# in the theorem is not zero, but it is on the right.
#
#########################################################
# Addendum: The PB-algorithm, a forerunner of AKS
# (BP -> Rajat Bhattacharjee and Prashant Pandey)
#########################################################
#
# Only for reasons of readebility
alias(composite=false): alias(invalid=false):
alias(prime=true):
# The BP primality test
isprime_BP := proc(n) local r,rmax,rhs,lhs;
# -----------------------------------------
if(n<2) then RETURN(invalid) fi;
if(n=2) then RETURN(prime) fi;
rmax := ithprime(floor(log(n)));
# The case r=2 trivally holds.
# So we might start with r=3.
r := 2;
do r := nextprime(r);
if(r > rmax) then RETURN(prime) fi;
if(irem(n,r) = 0) then RETURN(composite) fi;
rhs := Powmod(x^n+1,1,x^r+1,x) mod n;
lhs := Powmod(x+1,n,x^r+1,x) mod n;
if(lhs <> rhs) then RETURN(composite) fi;
od;
end:
# Der Beweis der Korrektheit des Algorithmus beruht auf
# folgender *Vermutung* von Pandey und Bhattacharjee
#
# (1) r>1 prim, n>1
# (2) r teilt nicht n
# (3) (x+1)^n = x^n+1 (mod x^r + 1, n)
# => n^2 = 1 (mod r) oder n ist prim
#
# Beweis der Korrektheit des Algorithmus (mit Vermutung):
#
# Pruefe fuer alle Primzahlen 2<=p<=(log(n)-ten Primzahl)
# (2) und (3). Ist dabei (2) oder (3) irgendwann einmal
# verletzt, dann ist n zusammengesetzt. Und wir sind fertig.
#
# Haben wir dagegen alle durchgeprueft und sind nicht ueber
# (2) oder (3) gestolpert, dann gilt nach dem Chinesischen
# Restsatz n^2 = 1 mod s, wobei s das Produkt der ersten
# log(n)+1 Primzahlen ist. Also ist n=1 oder n=s-1.
# Nach Definition von s ist n < (s-1). Und n > 1.
# Mais c'est absurde. Also ist n prim.
#
# Und der Algorithmus laeuft in polynomialer Zeit, weil die
# Kongruenzrelation (3) in polynomialer Zeit berechnet werden
# kann und die "Schleife" hoechstens log(n) mal durchlaufen wird.
#
#############################################################
# Post Scriptum
#############################################################
#
# Hopefully this Maple implementation makes the
# AKS-algorithm more accessible to the non-expert.
#
# It cannot, however, be considered as a full implementation,
# because of the use of the functions "nextprime", "ifactors"
# and "Powmod". Moreover, "nextprime" is build on the function
# "isprime", and the Maple docs say: "The function isprime is
# a probabilistic primality testing routine." :-(
#
# However, it would not be difficult to reimplement
# these functions in a completly deterministic manner without
# affecting the polyomial time bounds of the AKS-algorithm.
#
# Any comments, suggestions and bug reports are welcome.
# Please send mail to: peter AT luschny DOT de
#
#############################################################
# References
#############################################################
#
# M. Agrawal, N.Kayal, N. Saxena: PRIMES is in P. Preprint.
# http://www.cse.iitk.ac.in/news/primality.pdf
#
# D. Bernstein: An exposition of the Agrawal-Kayal-Saxena
# primality-proving theorem. Preprint.
# http://cr.yp.to/papers/aks.ps
#
# F. Morain: Primalité théorique et primalité pratique
# ou AKS vs. ECPP.
# http://www.lix.polytechnique.fr/Labo/Francois.Morain/aks-f.pdf
#
# P.Pandey, R.Bhattacharjee: Primality Testing.
# http://www.cse.iitk.ac.in/research/btp2001/primality.ps.gz
#
# J.F.Voloch: Improvemets to AKS.
# ftp://ftp.ma.utexas.edu/pub/papers/voloch/aks.pdf
#
########################## FIN ###############################