# OK to post homework
# Aurora Hiveley, 4/8/25, Assignment 20

Help := proc(): print(`Eq510(q,k,n,c,x)`): end:

with(combinat):


### project: 

i'll be working with pablo, lucy, and omar on implementing gabow's algorithm



### Eq510 

#Eq510(q,k,n,c,x): the value of Eq. (5.10) in Shor's paper (w/o the error term)
Eq510:=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)*u)/r),u=0..1)/r)^2): # leaving abs and square in according to dr. z's comments from class
end:


### approximation for probability of observing [c,k]

# evalc(abs(evalc(int(exp(2*Pi*I*x*u),u=0..1)))^2):
# simplify(%):

## output: 
# (1 - cos(2*Pi*x))/(2*Pi^2*x^2)




### histograms 

# Histogram56(128,86,4,3);
# Histogram56(256,86,4,3);

# Histogram56(512,86,4,3);
# Histogram56(128,86,4,9);
# Histogram56(256,172,4,3);
# Histogram56(128,86,16,9);
















### copied (& modified) 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), Ord(x,n)  `):end:
# read `C17.txt`:

#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. Try:
#Ord(2,101);
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:


#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.8) 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..63)]);
end: