#March 7, 2024 C15.txt
Help:=proc(): print(`S(x,y), aS(a), aSr(a,t) , OurPF(R,t) `):end:

#S(x,y): maps a pair of lists of numbers of the same length, n, say,
#and outputs a list of numbers of length 2*n
#[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 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,i=1..n),r=0..2*n-1)]:
end:

#aS(s): 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 n,x,y,X,Y,i,eq,var,r:
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)[1]:
[[seq(subs(var,x[i]),i=1..n)],[seq(subs(var,y[i]),i=1..n)]]:
end:

#aSr(a,t): solves the system of Ramanujan using his method with partial fractions
#outputs the rational function whose partial fraction would solve the system
aSr:=proc(a,t) local n,A,B,i,Top,Bot,f,eq,var:
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):
factor(subs(var,Top)/subs(var,Bot)):
end:

#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,X,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))]):

#R(t)= x[1]/(1-y[1]*t)+...+x[n]/(1-y[n]*t):
#x[1]/(1-y1*t)+ x[2]/(1-y[2]*t)= Top(t)/Bot(t)
#x(1)*(1-y[2]*t)+x[2](1-y[1]*t))/(1-y[1]*t)*(1-y[2]*t)):
#Top(t)=x[1]*(1-y[2]*t)+x[2](1-y[1]*t))+x[3]/(1-
[[seq(subs(t=1/Y[i],diff(Top,t)),i=1..n)],Y]:
end: