print(`For old help type: Help():`):
print(`For new help type: HelpNew():`):

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

HelpNew:=proc(): print(`Mult1(p,q,x,n,r)`): 
print(`Mult11(p,a,x,r), Power1(p,x,m,r,n), NewAKS(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:
#Dramatically improved by Josh Loftus (the one  who (so far)
#didn't murder anyone)
AKS:=proc(n) local r,a,asya,X, J,Asya:

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 trunc(evalf(4*log[2](n)^2)) while 
        ord1(n,r) <= trunc(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

Asya:=expand((X+J)^n-X^n-J) mod n:
Asya:=rem(Asya,X^r-1,X) mod n:


for a from 1 to trunc(2*sqrt(phi(r))*log[2](n)) do
  
asya:=subs(J=a,Asya) mod n:
 
if asya<>0 then
  RETURN(false):
fi:

od:

true:

end:



#NewAKS(n): the Agrawal-Kayal-Saxena famous
#primaity-testing algoritm:
#Dramatically improved by Josh Loftus (the one  who (so far)
#didn't murder anyone)
NewAKS:=proc(n) local r,a,asya,X, J,Asya:

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 trunc(evalf(4*log[2](n)^2)) while 
        ord1(n,r) <= trunc(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

#Asya:=expand((X+J)^n-X^n-J) mod n:
#Asya:=rem(Asya,X^r-1,X) mod n:


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

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

od:

true:

end:






















#Mult11(p,a,x,r): p*x^a mod x^r-1 mod n
Mult11:=proc(p,a,x,r) local i:

add(coeff(p,x,i)*x^((i+a) mod r),i=0..degree(p,x)):


end:

#Mult1(p,q,x,n,r): Inputs two polynomials p and q
#of degree <=r-1 and coefficients mod n
#and outputs their product mod x^r-1
Mult1:=proc(p,q,x,n,r) local i:

add(coeff(p,x,i)*Mult11(q,i,x,r) mod n
,i=0..degree(p,x)) mod n:
end:

#Power1(p,x,m,r,n): p(x)^m mod (n, x^r-1)
Power1:=proc(p,x,m,r,n) local P:

if m=0 then 
 RETURN(1):
elif m mod 2=1 then
Mult1( Power1(p,x,m-1,r,n),p,x,n,r):

else
Mult1( Power1(p,x,m/2,r,n), Power1(p,x,m/2,r,n)  ,x,n,r):
fi:


end: