Help:=proc(): print(`RE(a,n,p) , REd(a,n,p)`): 
print(`IsComposite(n,a), IsPrime(n,k) `):
end:

#RE(a,n,p): a^n mod p 
RE:=proc(a,n,p) local x:
option remember:

if n=0 then
 1:
elif n=1 then
  a mod p:
elif n mod 2 =0 then
 RE(a,n/2,p)^2 mod p
else
 RE(a,n-1,p)*a mod p:
fi:


end:

REd:=proc(a,n,p): a^n mod p: end:
REd1:=proc(a,n,p): a&^n mod p: end:

#A definite test to prove that n is composite
#using a as a trial witness. Returns
#true if n passed the test, FAIL if the test fails
IsComposite:=proc(n,a) local d, s,n1,r:
n1:=n-1:

for s from 1 while n1 mod 2 =0 do n1:=n1/2: od:
d:=n1:
s:=s-1:

if {evalb(RE(a,d,n)<>1), 
seq(evalb(RE(a,d*2^r,n)<>n-1),r=0..s-1)}={true} 
 then
   true:
else
  false:
fi:

end:

#IsPrime(n,k): The Miller-Rabin famous probabilstic
#primality test
#Inputs an integer n and a pos. integer k
#and performs IsComposite(n,a) k times for
#k random a's. If they all say FAIL then you
#return true (with high prob.) otherwise false
IsPrime:=proc(n,k) local a,i,ra:
ra:=rand(2..n-1):

for i from 1 to k do
 a:=ra():
 if IsComposite(n,a) then
    RETURN(false):
 fi:
od:

true:

end: