# OK to post homework
# Aurora Hiveley, 4/10/25, Assignment 21

Help := proc(): print(`CFx(x,k), GuessPCF(x)`): end:

# "turn it up to 11", aka we need more precision for this to work
Digits := 500:

# 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,L,y:
if evalf(x) <= 0 then 
    RETURN(FAIL):
fi:

L := []:
y := evalf(x):

while nops(L) < k and y<>infinity do 
    a := trunc(y):
    L := [op(L),a]:

    y := 1/(y-a):
od:

L:
end:



# GuessPCF(x): 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].

# example (from pablo): sqrt(2) has a continued fraction [1,2,2,2,2,...] so it has a repeating part. 
# Other quadratic irrationalities will also have a repeating part. So in the case of sqrt(2) for example you would return:  [[1],[2]]

GuessPCF:=proc(x) local L,L2,i,j,R,r,T:
if evalf(x) <= 0 then 
    RETURN(FAIL):
fi:

L := CFx(x,101):   # still pretty fast
L2 := []:   # initialize empty repeated section
i := 1:

# check if it repeats
while nops(L2)=0 and i<50 do 
  R := L[2..2+(i-1)]:   # candidate repeated section, has length i 
  r := floor(100/i):  # how many repeated sections fit into L? 
  T := [seq(op(R), j=1..r)]:  # test section 

  if L[2..i*r+1]=T then 
    L2 := R:
  fi:

  i++:

od:

# if no repeating section found from first 100 entries, return an error
if nops(L2)=0 then 
  L2 := [`error`]:
fi:


[[L[1]],L2]:
end:



# Find GuessPCF(n) for all non-square integers beween 2 and 100
# for n from 2 to 99 do 
#   if not type(sqrt(n),integer) then 
#     print(GuessPCF(sqrt(n)));
#   fi:
# od:

#[[1], [2]]
#[[1], [1, 2]]
#[[2], [4]]
#[[2], [2, 4]]
#[[2], [1, 1, 1, 4]]
#[[2], [1, 4]]
#[[3], [6]]
#[[3], [3, 6]]
#[[3], [2, 6]]
#[[3], [1, 1, 1, 1, 6]]
#[[3], [1, 2, 1, 6]]
#[[3], [1, 6]]
#[[4], [8]]
#[[4], [4, 8]]
#[[4], [2, 1, 3, 1, 2, 8]]
#[[4], [2, 8]]
#[[4], [1, 1, 2, 1, 1, 8]]
#[[4], [1, 2, 4, 2, 1, 8]]
#[[4], [1, 3, 1, 8]]
#[[4], [1, 8]]
#[[5], [10]]
#[[5], [5, 10]]
#[[5], [3, 2, 3, 10]]
#[[5], [2, 1, 1, 2, 10]]
#[[5], [2, 10]]
#[[5], [1, 1, 3, 5, 3, 1, 1, 10]]
#[[5], [1, 1, 1, 10]]
#[[5], [1, 2, 1, 10]]
#[[5], [1, 4, 1, 10]]
#[[5], [1, 10]]
#[[6], [12]]
#[[6], [6, 12]]
#[[6], [4, 12]]
#[[6], [3, 12]]
#[[6], [2, 2, 12]]
#[[6], [2, 12]]
#[[6], [1, 1, 3, 1, 5, 1, 3, 1, 1, 12]]
#[[6], [1, 1, 1, 2, 1, 1, 1, 12]]
#[[6], [1, 2, 2, 2, 1, 12]]
#[[6], [1, 3, 1, 1, 2, 6, 2, 1, 1, 3, 1, 12]]
#[[6], [1, 5, 1, 12]]
#[[6], [1, 12]]
#[[7], [14]]
#[[7], [7, 14]]
#[[7], [4, 1, 2, 1, 4, 14]]
#[[7], [3, 1, 1, 3, 14]]
#[[7], [2, 1, 6, 1, 2, 14]]
#[[7], [2, 2, 2, 14]]
#[[7], [2, 14]]
#[[7], [1, 1, 4, 1, 1, 14]]
#[[7], [1, 1, 1, 1, 1, 1, 14]]
#[[7], [1, 2, 7, 2, 1, 14]]
#[[7], [1, 2, 1, 14]]
#[[7], [1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14]]
#[[7], [1, 6, 1, 14]]
#[[7], [1, 14]]
#[[8], [16]]
#[[8], [8, 16]]
#[[8], [5, 2, 1, 1, 7, 1, 1, 2, 5, 16]]
#[[8], [4, 16]]
#[[8], [3, 3, 1, 4, 1, 3, 3, 16]]
#[[8], [2, 1, 2, 1, 2, 16]]
#[[8], [2, 2, 1, 7, 1, 2, 2, 16]]
#[[8], [2, 16]]
#[[8], [1, 1, 5, 5, 1, 1, 16]]
#[[8], [1, 1, 1, 1, 16]]
#[[8], [1, 1, 1, 16]]
#[[8], [1, 2, 1, 1, 5, 4, 5, 1, 1, 2, 1, 16]]
#[[8], [1, 3, 2, 3, 1, 16]]
#[[8], [1, 4, 1, 16]]
#[[8], [1, 7, 1, 16]]
#[[8], [1, 16]]
#[[9], [18]]
#[[9], [9, 18]]
#[[9], [6, 18]]
#[[9], [4, 1, 1, 4, 18]]
#[[9], [3, 1, 1, 1, 8, 1, 1, 1, 3, 18]]
#[[9], [3, 18]]
#[[9], [2, 1, 1, 1, 2, 18]]
#[[9], [2, 3, 3, 2, 18]]
#[[9], [2, 18]]
#[[9], [1, 1, 5, 1, 5, 1, 1, 18]]
#[[9], [1, 1, 2, 4, 2, 1, 1, 18]]
#[[9], [1, 1, 1, 4, 6, 4, 1, 1, 1, 18]]
#[[9], [1, 2, 3, 1, 1, 5, 1, 8, 1, 5, 1, 1, 3, 2, 1, 18]]
#[[9], [1, 2, 1, 18]]
#[[9], [1, 3, 1, 18]]
#[[9], [1, 5, 1, 1, 1, 1, 1, 1, 5, 1, 18]]
#[[9], [1, 8, 1, 18]]
#[[9], [1, 18]]





### sequence of periods

# S := []:
# for n from 2 to 99 do 
#   if not type(sqrt(n),integer) then 
#     S := [op(S),GuessPCF(sqrt(n))[2]];
#   fi:
# od:

# seq(nops(s),s in S);
# 1, 2, 1, 2, 4, 2, 1, 2, 2, 5, 4, 2, 1, 2, 6, 2, 6, 6, 4, 2, 1, 2, 
# 4, 5, 2, 8, 4, 4, 4, 2, 1, 2, 2, 2, 3, 2, 10, 8, 6, 12, 4, 2, 
# 1, 2, 6, 5, 6, 4, 2, 6, 7, 6, 4, 11, 4, 2, 1, 2, 10, 2, 8, 6, 
# 8, 2, 7, 5, 4, 12, 6, 4, 4, 2, 1, 2, 2, 5, 10, 2, 6, 5, 2, 8, 
# 8, 10, 16, 4, 4, 11, 4, 2

# this is A013943 in the OEIS 
# longest period occurs for sqrt(94)





# find the best rational approximations to e and Pi with denominator less than 10000000

# f := EvalCF(CFx(Pi,1)):
# k := 1:
# while denom(f) < 10000000 do 
#   f := EvalCF(CFx(Pi,k)):
#   k++:
# od:
# f;

# 80143857/25510582 for Pi 
# 28245729/10391023 for e 










### copied from C21.txt 

#C21.txt, April 10, 2025
Help21:=proc(): print(`Eq510(q,n,c,x), Eq510(256,33,10,5) , Histogram510(q,n,x), BA1(f,r) , BA(f) , CF(f) `):
CVG(f):
end:
read `C17.txt`:
Digits:=20:

Eq510:=proc(q,n,c,x) local r:
r:=Ord(x,n):
evalf(abs(1/r*int(exp(2*Pi*I*mods(r*c,q)/r*u),u=0..1))^2):
end:

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 hope[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) f=a1+1/(a2+ 1/(a3+ 1/....)
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):
RETURN([a,op(CF(1/(f-a)))]):
fi:
end:

#EvalCF(L): inputs a list of pos. integers and outputs the fraction whose cont. fraction it is
#(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:

#From C20.txt 
#C20.txt, April 7, 2025
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), Eq59(q,n,k,c,x), Histogram56(q,k,n,x)  `):end:
read `C17.txt`:

#Eq55(q,k,n,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,k,n,c,x): the value of Eq. (5.5) in Shor's paper
Eq56:=proc(q,n,k,c,x) local r,b:
r:=Ord(x,n):
if k>=r then
  RETURN(FAIL):
fi:

evalf(abs(add(exp(2*Pi*I*(b*r+k)*c/q),b=0..trunc((q-k-1)/r))/q)^2):
end:

#Eq56(q,k,n,c,x): the value of Eq. (5.6) in Shor's paper
Eq56:=proc(q,n,k,c,x) local r,b:
r:=Ord(x,n):
if k>=r then
  RETURN(FAIL):
fi:

evalf(abs(add(exp(2*Pi*I*(b*r+k)*c/q),b=0..trunc((q-k-1)/r))/q)^2):
end:

#Eq57(q,k,n,c,x): the value of Eq. (5.7) in Shor's paper
Eq57:=proc(q,n,k,c,x) local r,b:
r:=Ord(x,n):
if k>=r then
  RETURN(FAIL):
fi:

#evalf(abs(add(exp(2*Pi*I*(b*({r*c}_q)/q),b=0..trunc((q-k-1)/r))/q)^2):
evalf(abs(add(exp(2*Pi*I*b*(mods(r*c,q))/q),b=0..trunc((q-k-1)/r))/q)^2):
end:

#Eq58(q,k,n,c,x): the value of Eq. (5.8) in Shor's paper (w/o the error term)
Eq58:=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*b*(mods(r*c,q))/q),b=0..trunc((q-k-1)/r))/q)^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*(mods(r*c,q))/r),u=0..r/q*trunc((q-k-1)/r))/q)^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:

#end C20.txt