ModExp:=proc(x,

#Ord(x,n): the order of x mod n, the smallest r such that x^r=1 mod n
Ord:=proc(x,n) local p,r:
if igcd(x,n)<>1 then
  RETURN(FAIL):
fi:
p:=x:
for r from 1 while p mod n<>1 do 
p:=p*x mod n:
od:
r:
end:


FindFact1:=proc(n) local x,r,f:
x:=rand(4..n-1)():
if igcd(x,n)<>1 then
  RETURN(FAIL):
fi:


r:=Ord(x,n):
if r mod 2=1 then
  RETURN(FAIL):
fi:
if x^(r/2) mod n =-1 then
  RETURN(FAIL):
fi:
igcd(x^(r/2)-1,n):
end:

FindFact:=proc(n) local f,i:
for i from 1 to 10000 do
f:=FindFact1(n):
if f<>FAIL then
  RETURN(f):
fi:
od:
FAIL:
end:
FindFact(9*7*11);

                               7

Factorize:=proc(n) local f:
if isprime(n) then
 RETURN([n]):
fi:
f:=FindFact(n):
[op(Factorize(f)), op(Factorize(n/f))]:
end: