# OK to post homework
# Aurora Hiveley, 3/28/25, Assignment 17

Help := proc(): print(`MyOrd(x,n)`): end:

with(numtheory):


# if gcd(x,n), the order of x mod n always exists, since, thanks to Euler, for every such x
# 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. 

# MyOrd(x,n): (i) computes phi(n) (ii) finds its set of divisors, (iii) sorts them 
# (iv) keeps trying x^r mod n until it gets 1.

MyOrd := proc(x,n) local L,r:
if igcd(x,n)<>1 then
 RETURN(FAIL):
fi:

L := Factorize(phi(n),100):      # automatically sorted

for r in L do 
    if x^r mod n = 1 then 
        RETURN(r):
    fi:
od:

FAIL:
end:

# Take random large x and n and compare the time it takes.
# Why is even this more (possibly more) clever way of computing the order of x mod n 
# not crack RSA with a classical computer?

# generally, MyOrd is faster. 
# i assume this still does not work to crack RSA since as n gets very, very large 
# there are far too many factors of phi(n) to check




# A weighted directed graph with n vertices can be described as an n by n matrix with 
# non-negative entries (for the sake of simplicity, integers). 

# RWG(n,K): generates such a matrix (as a list-of-lists) where the diagonal is 0 
# (the distance from a vertex to itself is 0)
# note the matrix need not be symmetric since graph is directed 

RWG := proc(n,K) local ra,T,i,j:
ra := rand(0..K): 

for i from 1 to n do
    for j from 1 to n do 
        if i = j then 
            T[i,i] := 0:        # diagonal is all zeros
        else 
            T[i,j] := ra():
        fi:
    od:
od:

# reformat as list of lists 
[seq( [seq( T[i,j] , j=1..n)], i=1..n )]:
end:



# Implement Dijksra's algorithm that inputs a weighted directed graph (given as matrix) 
# and two vertices i, and j, and outputs (1) the length of the shortest path from i to j 
# (2) the actual path.

Dijkstra := proc(M,i,j) local P,U,D,n,k,l,u,elig:
n := nops(M):
P := [[]$n]: # paths from i to each vertex
U := {seq(1..n)}: # unvisited, step 1
D := [infinity$n]: # distances, step 2

# adjustments for starting vertex
D[i] := 0: # make distance to starting vertex 0

# identify current node, called k 
while nops(U)<>0 and min(D)<>infinity do 
    # determine unvisited at shortest distance -- step 3
    elig := [seq(D[u], u in U)]:
    member(min(elig),elig,'k'):
    k := U[k]:

    if k = j then 
        RETURN([D[j],P[j]]):
    fi:

    for l in U do         # step 4
        if M[k][l]<>0 and D[k] + M[k][l] < D[l] then   # update path and dist if shorter path found
            D[l] := D[k]+M[k][l]: 
            P[l] := [op(P[l]),k]:
        fi:
    od:

    U := U minus {k}:   # step 5

od:

[D[j],P[j]]:    # step 6
end:



# Give a few examples by using RWG(n,10);


## example 1
# R := RWG(5,10);
# output: R := [[0, 8, 3, 4, 5], [2, 0, 6, 9, 8], [2, 4, 0, 7, 9], [5, 9, 8, 0, 5], [3, 3, 4, 1, 0]]
# Dijkstra(R,2,1);
# output: [2, [2]]

## example 2
# R := RWG(3,10);
# output: R := [[0, 10, 10], [2, 0, 1], [10, 10, 0]]
# Dijkstra(R,1,2);
# output: [10, [1]]

## example 3 
# R := RWG(10,10);
# output: R := [[0, 2, 5, 2, 4, 2, 0, 3, 5, 3], [6, 0, 10, 2, 7, 10, 6, 1, 6, 6], [4, 7, 0, 1, 7, 3, 0, 1, 9, 7], [0, 7, 3, 0, 8, 6, 9, 5, 2, 5], [6, 9, 3, 3, 0, 8, 0, 5, 8, 5], [7, 4, 4, 6, 8, 0, 5, 9, 3, 7], [9, 3, 6, 4, 5, 3, 0, 4, 4, 7], [2, 10, 3, 10, 6, 10, 3, 0, 7, 10], [10, 7, 1, 5, 0, 5, 10, 10, 0, 10], [2, 9, 6, 0, 8, 5, 1, 3, 6, 0]]
# Dijkstra(R,5,9);
# output: [5, [5, 4]]





















### copied from C17.txt 

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