# Ok to post homework
# Lucy Martinez, 04-10-2025, Assignment 20



######################
#Problem 1: Write a procedure CFx(x,k) that
# inputs any positive number (given in decimals)
# and outputs the first k terms of its infinite continued fraction
CFx:=proc(x,k) local a,r,x1:
if k=0 then
 return([]):
fi:
a:=trunc(x):
r:=x-a:
if r=0 then
 return([a]):
else 
return([a,op(CFx(1/r,k-1))]):
fi:
end:



#Problem 2: Write a procedure GuessPCF(x)
# that inputs a quadratic irrationality x
# (e.g. sqrt(n) where n is not a square-free)
# and outputs two lists [L1,L2]
# such that the continued fraction of x is [L1,op(L2)infinity]
# Find GuessPCF(n) for all non-square integers beween 2 and 100


Digits:=100:

GuessPCF:=proc(x) local L1,bool,L,L2,p,i,k:
L1:=CFx(x,1):
for k from 30 to 50 do
L:=[op(2..-1,CFx(x,k))]:
if not FindPer(L)=FAIL then
 return([L1, FindPer(L)]):
fi:
od:
end:


#FindPer(L): finds the smallest period or return FAIL
FindPer:=proc(L) local p:
for p from 1 to trunc(nops(L)/2) do
if IsPer1(L,p) then
   RETURN([op(1..p,L)]):
fi:
od:
FAIL:
end:


#IsPer1(L,p): Is the list L periodic of period p
IsPer1:=proc(L,p) local i:
for i from 1 to nops(L)-p do
 if L[i]<>L[i+p] then
  RETURN(false):
fi:
od:
true:
end:


#The continued fractions are as follows
# where the first element is n for sqrt(n)
# 2, [[1], [2]]
# 3, [[1], [1, 2]]
# 5, [[2], [4]]
# 6, [[2], [2, 4]]
# 7, [[2], [1, 1, 1, 4]]
# 8, [[2], [1, 4]]
# 10, [[3], [6]]
# 11, [[3], [3, 6]]
# 12, [[3], [2, 6]]
# 13, [[3], [1, 1, 1, 1, 6]]
# 14, [[3], [1, 2, 1, 6]]
# 15, [[3], [1, 6]]
# 17, [[4], [8]]
# 18, [[4], [4, 8]]
# 19, [[4], [2, 1, 3, 1, 2, 8]]
# 20, [[4], [2, 8]]
# 21, [[4], [1, 1, 2, 1, 1, 8]]
# 22, [[4], [1, 2, 4, 2, 1, 8]]
# 23, [[4], [1, 3, 1, 8]]
# 24, [[4], [1, 8]]
# 26, [[5], [10]]
# 27, [[5], [5, 10]]
# 28, [[5], [3, 2, 3, 10]]
# 29, [[5], [2, 1, 1, 2, 10]]
# 30, [[5], [2, 10]]
# 31, [[5], [1, 1, 3, 5, 3, 1, 1, 10]]
# 32, [[5], [1, 1, 1, 10]]
# 33, [[5], [1, 2, 1, 10]]
# 34, [[5], [1, 4, 1, 10]]
# 35, [[5], [1, 10]]
# 37, [[6], [12]]
# 38, [[6], [6, 12]]
# 39, [[6], [4, 12]]
# 40, [[6], [3, 12]]
# 41, [[6], [2, 2, 12]]
# 42, [[6], [2, 12]]
# 43, [[6], [1, 1, 3, 1, 5, 1, 3, 1, 1, 12]]
# 44, [[6], [1, 1, 1, 2, 1, 1, 1, 12]]
# 45, [[6], [1, 2, 2, 2, 1, 12]]
# 46, [[6], [1, 3, 1, 1, 2, 6, 2, 1, 1, 3, 1, 12]]
# 47, [[6], [1, 5, 1, 12]]
# 48, [[6], [1, 12]]
# 50, [[7], [14]]
# 51, [[7], [7, 14]]
# 52, [[7], [4, 1, 2, 1, 4, 14]]
# 53, [[7], [3, 1, 1, 3, 14]]
# 54, [[7], [2, 1, 6, 1, 2, 14]]
# 55, [[7], [2, 2, 2, 14]]
# 56, [[7], [2, 14]]
# 57, [[7], [1, 1, 4, 1, 1, 14]]
# 58, [[7], [1, 1, 1, 1, 1, 1, 14]]
# 59, [[7], [1, 2, 7, 2, 1, 14]]
# 60, [[7], [1, 2, 1, 14]]
# 61, [[7], [1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14]]
# 62, [[7], [1, 6, 1, 14]]
# 63, [[7], [1, 14]]
# 65, [[8], [16]]
# 66, [[8], [8, 16]]
# 67, [[8], [5, 2, 1, 1, 7, 1, 1, 2, 5, 16]]
# 68, [[8], [4, 16]]
# 69, [[8], [3, 3, 1, 4, 1, 3, 3, 16]]
# 70, [[8], [2, 1, 2, 1, 2, 16]]
# 71, [[8], [2, 2, 1, 7, 1, 2, 2, 16]]
# 72, [[8], [2, 16]]
# 73, [[8], [1, 1, 5, 5, 1, 1, 16]]
# 74, [[8], [1, 1, 1, 1, 16]]
# 75, [[8], [1, 1, 1, 16]]
# 76, [[8], [1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16]]
# 77, [[8], [1, 3, 2, 3, 1, 16]]
# 78, [[8], [1, 4, 1, 16]]
# 79, [[8], [1, 7, 1, 16]]
# 80, [[8], [1, 16]]
# 82, [[9], [18]]
# 83, [[9], [9, 18]]
# 84, [[9], [6, 18]]
# 85, [[9], [4, 1, 1, 4, 18]]
# 86, [[9], [3, 1, 1, 1, 8, 1, 1, 1, 3, 18]]
# 87, [[9], [3, 18]]
# 88, [[9], [2, 1, 1, 1, 2, 18]]
# 89, [[9], [2, 3, 3, 2, 18]]
# 90, [[9], [2, 18]]
# 91, [[9], [1, 1, 5, 1, 5, 1, 1, 18]]
# 92, [[9], [1, 1, 2, 4, 2, 1, 1, 18]]
# 93, [[9], [1, 1, 1, 4, 6, 4, 1, 1, 1, 18]]
# 94, [[9], [1, 2, 3, 1, 1, 5, 1, 8, 1, 5, 1, 1, 3, 2, 1, 18]]
# 95, [[9], [1, 2, 1, 18]]
# 96, [[9], [1, 3, 1, 18]]
# 97, [[9], [1, 5, 1, 1, 1, 1, 1, 1, 5, 1, 18]]
# 98, [[9], [1, 8, 1, 18]]
# 99, [[9], [1, 18]]



#Problem 3: Is the sequence of lengths of their periods in the OEIS?
# Which is the largest period

# Answer: Yes, it is :) here is the A no:
#  A013943
#  So far, I saw it might be 146. There are 10,000 elements
#  in this list

#Problem 4: Find the best rational approximations to e and Pi
# with denominator less than 10000000




####################################from previous classes:
#C21.txt: April 10, 2025
#Prime Factorization, Section 5 of Peter W. Shor's paper
Help21:=proc():
print(`Eq510(q,n,c,x), Histogram510(q,n,x),BA1(f,r),BA(f)`):
print(`CF(f),EvalCF(L), CVG(f)`):
end:



#Eq510(q,n,c,x): the value of Eq. (5.10) in Shor's paper
Eq510:=proc(q,n,c,x) local r:
r:=Ord(x,n):
evalf(abs(int(exp(2*Pi*I*(mods(r*c,q))/r*u),u=0..1)/r)^2):
end:

#Histogram510(q,n,x): the plot of [c,Eq10(q,n,c,x)] for c from 1 to q-1
Histogram510:=proc(q,n,x) local c:
plot([seq([c,Eq510(q,n,c,x)],c=1..q-1)]);
end:


#BA1(f,r): the closest fraction d/r to f, where r is fixed
BA1:=proc(f,r) local b:
b:=round(f*r)/r:
[b, evalf(abs(b-f)*denom(f))]:
end:


#BA(f): the best fraction d/r where r is between 1 and sqrt(denom(f))
BA:=proc(f) local r,cha,rec,hope:
cha:=BA1(f,1)[1]:
rec:=BA1(f,1)[2]:
for r from 2 to trunc(sqrt(denom(f))) do
 hope:=BA1(f,r):
 if hope2[2]<rec then
  cha:=hope[1]:
  rec:=hope[2]:
 fi:
od:
cha,rec:
end:


#CF(f): inputs an integer or fraction outputs the list consisting of
# the continued-fraction
#11/4=2+3/4=2+1/(4/3)=2+1/(1+1/3) 
CF:=proc(f) local a:
if not (type(f,rational) and f>0) then
 return(FAIL):
fi:

if type(f,integer) then
 return([f]):
else
a:=trunc(f):
[a,op(CF(1/(f-a)))]:
fi:
end:


#EvalCF(L): inputs a list of pos. integers and outputs the fraction 
# whose continued fraction it is
# the reverse of CF(f)
EvalCF:=proc(L):
if nops(L)=1 then
 L[1]:
else
 L[1]+1/EvalCF([op(2..nops(L),L)]):
fi:
end:


#CVG(f): inputs a rational number f and outputs its list of CONVERGENTS
CVG:=proc(f) local L,i:
L:=CF(f):
[seq(EvalCF([op(1..i,L)]),i=1..nops(L))]:
end:


###############################

#C20.txt: April 07, 2025
#Prime Factorization, Section 5 of Peter W. Shor's paper
#https://sites.math.rutgers.edu/~zeilberg/EM25/shor.pdf

Help20:=proc():
print(`Eq55(q,n,k,c,x), Eq56(q,n,k,c,x), Eq57(q,n,k,c,x), Eq58(q,n,k,c,x),`):
print(`Eq59(q,k,n,c,x), Histogram56(q,k,n,c)`):
end:


#Eq55(q,n,k,c,x): the value of Eq. (5.5) in Shor's paper
Eq55:=proc(q,n,k,c,x) local a,su:
su:=0:
for a from 0 to q-1 do
 if x^a-x^k mod n=0 then
  su:=evalc(su+exp(2*Pi*I*a*c/q)):
 fi:
od:
evalf((abs(su/q))^2):
end:


#Eq56(q,n,k,c,x): the value of Eq. (5.6) in Shor's paper
Eq56:=proc(q,n,k,c,x) local b,r:
r:=Ord(x,n):
evalf(abs(1/q*add(exp(2*Pi*I*(b*r+k)*c/q), b=0..trunc((q-k-1)/r)))^2):
end:


#Eq57(q,n,k,c,x): the value of Eq. (5.7) in Shor's paper
Eq57:=proc(q,n,k,c,x) local b,rc,r:
r:=Ord(x,n):
# mods(r*c,q) is the same as computing r*c mod q
evalf(abs(1/q*add(exp(2*Pi*I*b*mods(r*c,q)/q), b=0..trunc((q-k-1)/r)))^2):
end:



#Eq58(q,n,k,c,x): the integral of Eq. (5.8) in Shor's paper
Eq58:=proc(q,n,k,c,x) local b,rc,r:
r:=Ord(x,n):
# mods(r*c,q) is the same as computing r*c mod q
evalf(abs(1/q*int(exp(2*Pi*I*b*mods(r*c,q)/q), b=0..trunc((q-k-1)/r)))^2):
end:

#Eq59(q,k,n,c,x): the value of Eq. (5.9) in Shor's paper (w/o the error term)
Eq59:=proc(q,n,k,c,x) local r,b:
r:=Ord(x,n):
if k>=r then
  RETURN(FAIL):
fi:
evalf(abs(int(exp(2*Pi*I*u*(mods(r*c,q))/r),u=0..r/q*trunc((q-k-1)/r))/r)^2):
end:


#Added after class
#Histogram56(q,k,n,c): the plot of [c,Eq56(q,k,n,c,x)] for c from 1 to q-1
Histogram56:=proc(q,k,n,x) local c:
plot([seq([c,Eq56(q,k,n,c,x)],c=1..q-1)]);
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: