#C17.txt, March 27, 2025. Shor's algorithm
Help17:=proc(): print(`Ord(x,n) , FindFact1(n), FindFact(n,K), Factorize(n,K) , ModExp(x,r,n) `): end:

#Ord(x,n): inputs a large pos. integer n and x<n m such that gcd(x,n) outputs (by brute force)
#the order of x mod n, i.e. the smallest pos. integer r such that x^r=1
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(n): tries to find a non-trivial factor of n or returns FAIL
FindFact1:=proc(n) local x,r:
x:=rand(2..n-1)():
r:=Ord(x,n):

if r=FAIL then
  RETURN(r):
fi:

if r mod 2=1 then
  RETURN(FAIL):
fi:

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

#FindFact(n,K): tries to find a factor of n using FindFact1(n) K times or returns FAIL (with prob.
#less than 1/2^K
FindFact:=proc(n,K) local i,p:
for i from 1 to K do
 p:=FindFact1(n):
 if p<>FAIL then
   RETURN(p):
 fi:
od:
FAIL:
end:

#Factorize(n,K): inputs a pos. integer n and outputs the list of prime factors
Factorize:=proc(n,K) local L,p:
if isprime(n) then
  RETURN([n]):
fi:

p:=FindFact(n,K):
if p=FAIL then
 RETURN(FAIL):
fi:
sort([ op(Factorize(p,K)), op(Factorize(n/p,K))    ]):
end:

#ModExp(x,r,n): x^r mod n the fast way where r is given in reverse binary [2^0,2^1,2^2,...]
ModExp:=proc(x,r,n) local k,T,p,i:
k:=nops(r)+1:
T[0]:=1:
T[1]:=x:
for i from 1 to k-1 do
 T[i+1]:=T[i]^2 mod n: 
od:
p:=1:
for i from 0 to nops(r)-1 do
if r[i+1]=1 then
  p:=p*T[i+1] mod n:
fi:
od:
p:
end: