# Ok to post homework
# Lucy Martinez, 03-30-2025, Assignment 17

with(combinat):
with(NumberTheory):

######################
#Problem 1: Tijil Kiran noticed that if gcd(x,n)=1, then
# the order of x mod n always exists, since, thanks to Euler,
# for every such x relatively prime to n, x^phi(n) mod n =1.
# Hence the order must be a divisor of phi(n).
# The Maple package numtheory has both phi and divisors.
# Write a procedure MyOrd(x,n) that
# (i) computes phi(n)
# (ii) finds its set of divisors,
# (iii) sorts them
# (iv) keeps trying x^r mod n until it gets 1.

# Take random large x and n and compare the time it takes.
# Some Data:

# x:=rand(1000000..10000000)():
# n:=rand(1000000..10000000)():
# time(MyOrd(x,n))=0.125
# time(Ord(x,n))=0.609

# x:=rand(500..10000000)():
# n:=rand(500..10000000)():
# time(MyOrd(x,n))= 0.156
# time(Ord(x,n))=2.140

MyOrd:=proc(x,n) local ph,D,r:
ph:=phi(n):
D:=Divisors(ph):
for r in D do
 if x^r mod n = 1 then
  return r:
 fi:
od:
end:


#Problem 2: A weighted directed graph with n vertices
# can be described a pair [n,E] (like the one generated by procedure RG(n,p)
# of C2.txt), together with a table T defined on E
# where for any member e of E T[e] is the distance.
# Write a procedure, RWG(n,p,K) whose output is a triple [n,E,op(T)]
# where T is the table of distances,
# where the distances are random integers from 1 to K.

# For example the graph with vertices {1,2,3} and set of edges {{1,2},{1,3}}
# and the distance between 1 and 2 is 5 and between 1 and 3 is 2 is expressed as:
# [3, {{1, 2}, {1, 3}}, table([{1, 3} = 2, {1, 2} = 5])]
RWG:=proc(n,p,K) local ra, G,E,e,T:
ra:=rand(1..K):
G:=RG(n,p):
E:=G[2]:
for e in E do
T[e]:=ra():
od:
[n,E,op(T)]:
end:



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


#C2.txt: Jan. 27, 2025 Exp Math (Dr. Z.)
Help2:=proc():
print(`LC(p), RG(n,p), Cliques(G,k)`):
end:


#LC(p): inputs a rational number between 0 and 1 and outputs true with prob. p
LC:=proc(p) local a,b,ra:
if not (type(p,fraction) and p>= 0 and p<=1) then
 return fail:
fi:

a:=numer(p):
b:=denom(p):

ra:=rand(1..b)():

if ra<=a then
 true:
else
 false:
fi:
end:


RG:=proc(n,p) local E,i,j:
E:={}:
for i from 1 to n do
 for j from i+1 to n do
  if LC(p) then
   E:=E union {{i,j}}:
  fi:
 od:
od:
[n,E]:
end: