Help:=proc(): print(`ModPol(p,x,h), ord1(a,r)`): 
print(`Check31(M), AKS(n)`):
end:


#ModPol(p,x,h): inputs a polynomial p in the variable
#x and another polynomial h, in x, finds
#p mod h
ModPol:=proc(p,x,h) 
rem(p,h,x):

end:

with(numtheory):

#ord1(a,r): the smallest pos. integer k such that
#a^k mod r =1 (gcd(a,r)=1)
ord1:=proc(a,r) local S,i:
if gcd(a,r)<>1 then
  RETURN(1):
fi:

S:=sort(divisors(phi(r))):


for i from 1 to nops(S) do

 if a^S[i] mod r =1 then
    RETURN(S[i]):
 fi:
od:
ERROR(`Nothing found`):


end:

#Check31(M) checks Lemma 3.1. in the AKS paper
#for m between 7 and M 
Check31:=proc(M) local m:
evalb({seq(evalb(lcm(seq(i,i=1..m))>=2^m),m=7..M)}={true}):
end:

#AKS(n): the Agrawal-Kayal-Saxena famous
#primaity-testing algoritm:

AKS:=proc(n) local r,a,asya,X:

if not type(n,integer) or n<=1 then
 ERROR(`bad input`):
fi:

#step 1
for r from 2 to trunc(evalf(log(n)/log(2))) do
 if type(root(n,r), integer) then 
RETURN(false)
   
 fi:
od:

#step 2
for r from 2 while ord1(n,r) <= evalf(4*log[2](n)^2) do od:

#step 3
for a from 2 to r do
 if 1<igcd(a,n) and igcd(a,n)<n then
    RETURN(false):
 fi:
od:

#step 4

if n<=r then
  RETURN(true):
fi:

#step 5

for a from 1 to trunc(2*sqrt(phi(r))*log[2](n)) do
  
asya:=(X+a)^n-X^n-a :

asya:=rem(asya,X^r-1,X) mod n:

if asya<>0 then
  RETURN(false):
fi:

od:

true:

end: