#OK to post homework
#Lucy Martinez, 03-10-2024, Assignment 15


# Problem 2: Use the corrected procedure aSr(a) in C15.txt to check that Ramanujan's example
# on the second page of his seminal paper is correct.
# [Warning: Maple is not so good with simplifying expressions with radicals,
#  you can either do "evalf", or to check that two algebraic expressions,
#  expressed in radicals, a and b, represent the same number, do : simplify(a-b)
# Answer:
# L:=evalf(aSr([2,3,16,31,103,235,674,1669,4526,11595]));
# outputs
# L := [[-0.6000000000, 2.023606792, 1.576393202, 1.288854381,-2.288854383],[-1., 2.618033984, 0.3819660114, -1.618033991, 0.6180339890]]

# The first list, L[1], is [-0.6000000000, 2.023606792, 1.576393202, 1.288854381,-2.288854383]
# which coincides with the [x,y,z,u,v] list of the paper which is:

# L1:=evalf([-3/5,(18+sqrt(5))/10,(18-sqrt(5))/10,-(8+sqrt(5))/(2*sqrt(5)),(8-sqrt(5))/(2*sqrt(5))]);
# which outputs
# L1 := [-0.6000000000, 2.023606798, 1.576393202, -2.288854382, 1.288854382]

# The second list, L[2] is [-1., 2.618033984, 0.3819660114, -1.618033991, 0.6180339890]
# which coincides with the [p,q,r,s,t] list of the paper which is:

# L2:=evalf([-1,(3+sqrt(5))/2,(3-sqrt(5))/2,(sqrt(5)-1)/2,-(sqrt(5)+1)/2]); which outputs
# L2 := [-1., 2.618033988, 0.381966012, 0.6180339880, -1.618033988]
# NOTE: the last two elements of L1 and L2 are switched when comparing it to L



# Problem 3: Write a finite field analog of procedures, S(x,y), OurPF(R,t), and aSr(a),
# call them qS(x,y,q), qOurPF(R,t,q), and qaSr(a,q) where everyhing is done modulo q
# Test it with q=11 and the example on p. 134 of the book, in other words, 
# find qaSr([2,8,4,5,3,2],11);
# [Hint: to find the roots of the denominator just do it by brute force,
#  collecting all i such that B(i) mod q=0. To get the reciprocal of a mod q do a&^(-1) mod q]

# ANSWER: Executing qaSr([2,8,4,5,3,2],11); outputs
# [[4, 2, 7], [3, 5, 9]]
# which coincides with the book :)

## HERE are the procedures:

# qS(x,y,q): Maps a pair of lists of numbers of the same length, say n, and outputs
# a list of numbers of length 2n:
# [x[1]+...+x[n], x[1]*y[1]+...+x[n]*y[n], ... , x[1]*y[1]^(2*n+1)+...+x[n]*y[n]^(2*n+1) ]
# modulo q
qS:=proc(x,y,q) local i,r,n:
if not (type(x,list) and type(y,list) and nops(x)=nops(y)) then
 return(FAIL):
fi:
n:=nops(x):

[seq(add(x[i]*y[i]^r mod q,i=1..n),r=0..2*n-1)]:
end:


#qOurPF(R,t,q): Inputs a rational function of t, and outputs
# the partial fraction decomposition over the complex numbers of the form
# [a1,r1],...,[an,rn] where 1/r1,..., 1/rn are the complex roots of the bottom
qOurPF:=proc(R,t,q) local n,Top,Bot,zer,Y,i,d:
Top:=numer(R):
Bot:=denom(R):

zer:=[]:
for i from 0 to q-1 do
 if subs(t=i,Bot) mod q = 0 then
  zer:=[op(zer),i]:
 fi:
od:
n:=nops(zer): 
Y:=sort([seq( zer[i]&^(-1) mod q,i=1..nops(zer))]):
[[seq(-Y[i]*subs(t= Y[i]&^(-1) mod q ,Top)/subs(t=Y[i]&^(-1) mod q,diff(Bot,t)) mod q,i=1..n)],Y]:
end:


# aSr(a,q): solves the system of Ramanujan using his method with partial fractions
# outputs the rational function  whose partial fractions would solve the system 
# all done modulo q
qaSr:=proc(a,q) local n,A,B,i,Top,Bot,f,eq,var,t:
if not( type(a,list) and nops(a) mod 2=0) then
 RETURN(FAIL):
fi:
n:=nops(a)/2:
Top:=add(A[i+1]*t^i,i=0..n-1):
Bot:=1+add(B[i]*t^i,i=1..n):
f:=expand(Top-Bot*add(a[i]*t^(i-1),i=1..2*n)) mod q:
eq:={seq(coeff(f,t,i),i=0..2*n-1)}:
var:={seq(A[i],i=1..n),seq(B[i],i=1..n)}:
var:=solve(eq,var):
qOurPF(normal(subs(var,Top)/subs(var,Bot)),t,q):
end:



##################################################### from class:
Help15:=proc():
print(`S(x,y), aS(a), aSr(a),OurPF(R,t)`):
end:


# S(x,y): Maps a pair of lists of numbers of the same length, say n, and outputs
# a list of numbers of length 2n:
# [x[1]+...+x[n], x[1]*y[1]+...+x[n]*y[n], ... , x[1]*y[1]^(2*n+1)+...+x[n]*y[n]^(2*n+1) ]
S:=proc(x,y) local i,r,n:
if not (type(x,list) and type(y,list) and nops(x)=nops(y)) then
 return(FAIL):
fi:
n:=nops(x):
[seq(add(x[i]*y[i]^r,i=1..n),r=0..2*n-1)]:
end:


# antiS
# aS(a): reverse S(x,y) using Maple's solve command
# inputs a list of length 2*n and outputs lists x,y of length n, such that S(x,y)=a
aS:=proc(a) local x,y,X,Y,i,eq,var,r,n:
if not (type(a,list) and nops(a) mod 2 = 0) then
 return(FAIL):
fi:
n:=nops(a)/2:
eq:={seq(add(x[i]*y[i]^r,i=1..n)=a[r+1],r=0..2*n-1)}:
var:={seq(x[i],i=1..n),seq(y[i],i=1..n)}:
var:=solve(eq,var):
[[seq(subs(var,x[i]),i=1..n)],[seq(subs(var,y[i]),i=1..n)]]:
end:


# aSr(a): solves the system of Ramanujan using his method with partial fractions
# outputs the rational function  whose partial fractions would solve the system 
aSr:=proc(a) local n,A,B,i,Top,Bot,f,eq,var,t:
if not( type(a,list) and nops(a) mod 2=0) then
 RETURN(FAIL):
fi:
n:=nops(a)/2:

Top:=add(A[i+1]*t^i,i=0..n-1):

Bot:=1+add(B[i]*t^i,i=1..n):

#Top/Bot=add(a[i]*t^(i-1),i=1..2*n):

f:=expand(Top-Bot*add(a[i]*t^(i-1),i=1..2*n)):
eq:={seq(coeff(f,t,i),i=0..2*n-1)}:
var:={seq(A[i],i=1..n),seq(B[i],i=1..n)}:
var:=solve(eq,var):
OurPF(normal(subs(var,Top)/subs(var,Bot)),t):
end:



# our own partial fraction decomposition
#OurPF(R,t): Inputs a rational function of t, and outputs
# the partial fraction decomposition over the complex numbers of the form
# [a1,r1],...,[an,rn] where 1/r1,..., 1/rn are the complex roots of the bottom
OurPF:=proc(R,t) local n,Top,Bot,zer,Y,i:
Top:=numer(R):
Bot:=denom(R):
zer:=[solve(Bot,t)]:
n:=nops(zer):
Y:=sort([seq(1/zer[i],i=1..nops(zer))]):
[[seq(-Y[i]*subs(t=1/Y[i],Top)/subs(t=1/Y[i],diff(Bot,t)),i=1..n)],Y]:
end: