##Written by Kellen Myers, Rutgers University, as homework for Dr. Z.'s Experimental Mathematics class
# New AKS implementation myAKS(n)
# Old AKS implementation AKS(n)

with(numtheory):
with(numtheory,order):

myord1:=proc(n,r)
 if gcd(n,r)<>1 then return(1): fi:
 return(order(r,n)):
end:

myAKS:=proc(n) local r,a,g,A,p,p2,i:
 #PART 0
 if not type(n,integer) or n<=1 then error(`bad input`): fi:

 #PART 1
 for r from 2 to trunc(log[2](n)) do
  if type(root(n,r), integer) then return(false): fi:
 od:
 
 #PART 2 & 3 & 4 -- all these steps can be combined
 r:=2:
 while r <= trunc(evalf(4*log[2](n)^2)) do
   r:=r+1: g:=igcd(r,n):                                # STEP 2
   if evalb(g <> 1 and g <> n) then return(false): fi:  # STEP 3
 od:
 while myord1(n,r)<=evalf(4*log[2](n)^2) do
   if r>=n then return(true): fi:                       # STEP 4
   r:=r+1: g:=igcd(r,n):                                # STEP 2
   if evalb(g <> 1 and g <> n) then return(false): fi:  # STEP 3
 od:

 #PART 5 - includes repeated squaring, etc.
 B2=convert(n,base,2);
 p:=0:
 p2:=x+A:
 for i from 1 to nops(B2) do
 p2:=modp(rem(modp(p2^2,n),x^r-1,x),n):
 if B2[i]=1 then p:=p*p2: fi:
 od:
 for a from 1 to trunc(2*sqrt(phi(r))*log[2](n)) do
  if modp(subs(A=a,p),n) <> 0 then return(false): fi:
 od:

 #PART 6
 return(true):

end:



# Our old AKS, including improvements from Josh

#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: