#OK to post homework
#Ramesh Balaji,3/10/2024,hw15

# Part 1
# Note A(t)/B(t) = sum(Ci(t)/(1-tyi), i=1..n).
# Note that (1-tyj)A(t)/B(t) = sum((1-tyj)Ci(t)/(1-tyi)) for any yi
# Then we have that if t=1/yi then (1-tyj)A(t)/B(t) = C(t) so note that
# lim t=1/yi (1-tyj)A(t)/B(t) = C(t) so C(t) = -yiA(t)/B'(t) by L'hopitals so we
# are done.

# Part 2
#corrected after class
#OurPF(R,t): inputs a rational function of t outputs the
#partial fraction decomopition 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:


#aSr(a): solves the system of Ramanujan using his method with partial fractions
#outputs the rational function whose partial fraction would solve the system
#modified after class
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:

d:=aSr([2,3,16,31,103,235,674,1669,4526,11595]);
for i from 0 to 9 do
    print(evalf(d[iquo(i,5) + 1][irem(i,5) + 1]))
od:


# Part 3
qS:=proc(x,y, q) local n,i,r:
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) mod q,r=0..2*n-1)]:
end:

qOurPF:=proc(R,t,q) local n,Top,Bot,zer,Y,i:
Top:=numer(R):
Bot:=denom(R):

zer := [];
    print(q);
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(1/zer[i] mod q,i=1..nops(zer))]):
[[seq( (-Y[i]*subs(t=1/Y[i],Top)/subs(t=1/Y[i],diff(Bot,t))) mod q,i=1..n)],Y]:
end:

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) :
    print(Top, Bot);

#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)):
    print(f);
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: