#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 10
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.
Help := proc(): print(`Help(), GuessAlg(L,P,x, MaxDegP), GuessAlg1(L,P,x, DegP), GuessAlg11(L,P,x, DegP, Degx), DiagToSeq(f,x,y,m), DiagToAlg(R,x,y,P,t,maxdegP), GFtoSeq(f,q,m), AlgToSeq(F, P, x, K, hint)`): end:


## Stuff from c10.txt:

end:
#C10.txt: starting the Alg. Formal Power Series
#1+1!*x+2!*x^2+3!*x^3+ ......

#GuessAlg(L,P,x,MaxDegP): Inputs a sequence of numbers L
#and tries to guess  an alg. eq. of the from
#F(P,x)=0. where P is the symbol for the generating
#function of L, P(x)=add(L[i+1]*x^i,i=0..infinity)
#of max. degree MaxDegP
GuessAlg:=proc(L,P,x, MaxDegP) local g,DegP:

for DegP from 1 to MaxDegP do
 g:=GuessAlg1(L,P,x,DegP):
 if g<>FAIL then
     RETURN(g):
 fi:
od:

FAIL:

end:


#GuessAlg1(L,P,x,DegP): Inputs a sequence of numbers L
#and tries to guess  an alg. eq. of the from
#F(P,x)=0. where P is the symbol for the generating
#function of L, P(x)=add(L[i+1]*x^i,i=0..infinity)
#of degree DegP in P
GuessAlg1:=proc(L,P,x, DegP) local g,Degx:



for Degx from 0 while (1+Degx)*(1+DegP)<nops(L)-6 do

g:=GuessAlg11(L,P,x, DegP,Degx):

if g<>FAIL then
  RETURN(g):
fi:

od:

FAIL:

end:



#GuessAlg11(L,P,x,DegP, Degx): Inputs a sequence of numbers L
#and tries to guess  an alg. eq. of the from
#F(P,x)=0. where P is the symbol for the generating
#function of L, P(x)=add(L[i+1]*x^i,i=0..infinity)
#of degree DegP in P and Degx in x
GuessAlg11:=proc(L,P,x, DegP, Degx) local g, F,a, i,j,
var,P1,F1,eq, var1,v:

F:=add(add(a[i,j]*x^i*P^j,i=0..Degx),j=0..DegP):

var:={seq(seq(a[i,j],i=0..Degx),j=0..DegP)}:

P1:=add(L[i]*x^(i-1),i=1..nops(L)):

F1:=expand(subs(P=P1,F)):

eq:={seq( coeff(F1,x,i)=0, i=0..nops(L)-1)}:

var1:=solve(eq,var):

F:=subs(var1,F):

if F=0 then
 RETURN(FAIL):
else
F:=subs({seq(v=1, v  in var)},F):

F:=add(factor(coeff(F,P,i))*P^i,i=0..degree(F,P)):
RETURN(F):
fi:



end:

#############
# Problem 2 #
#############

DiagToSeq := proc(f,x,y,m):
    [seq(coeff(coeff(mtaylor(f, [x,y], 2*m+2), x, i), y, i), i=0..m)]:
end:

DiagToAlg := proc(R,x,y,P,t,maxdegP) local L:
    L := DiagToSeq(R, x, y, maxdegP*4+12):
    return GuessAlg(L, P, t, maxdegP):
end:

## Answers to the questions:

# 1. F(x, P) = 1 + (x-1)P
# 2. F(x, P) = 1/4 + (x-1/4)P^2
# 3. F(x, P) = -1/9 + (1/9 + 2/9 x + x^2)P^2
# 4. F(x, P) = 1 + 2(4t-3)(4t^2-8t+1)P^2 + (8t+5)(4t^2-8t+1)^2P^4

#############
# Problem 3 #
#############

# From HW 5:
GFtoSeq := proc(f,q,m):
    [seq(coeff(taylor(f, q, m+1), q, i), i=0..m)]:
end:


# In some cases there may be more than one generating function that satisfies the
# equation. The "hint" argument allows you to specify the first few
# terms of the sequence (as a list) to ensure that you get the right
# sequence.

# N.B.: I did not use the hint. The hint made the problem much harder :(

AlgToSeq := proc(F, P, x, K, hint) local a, f, i, j, poly, eqns, solns, start:
    poly := add(a[i]*x^i, i=0..degree(F, x)):
    eqns := {seq(coeff(subs(P=poly, F), x, i)=0 , i=0..degree(F, x))}:
    if (hint <> NULL) then:
        eqns := eqns union {seq(a[i] = hint[i+1], i=0..nops(hint)-1)}:
    fi:
    solns := solve(eqns):
    poly := subs(solns, poly);
    f := poly:


    start := degree(poly):


    for i from start+1 to K do:
        poly := f+a[i]*x^i:
        eqns := {seq(coeff(subs(P=poly, F), x, j)=0, j=0..i)}:
        solns := solve(eqns):
        f := subs(solns, poly):
    od:
    return [seq(coeff(f,x,i), i=0..K)]:
end:

#############
# Problem 5 #
#############

# I'm not sure if I have really solved the problem in its full
# generality (because I don't understand what the full generality is), but
# at least it correctly solves the special case specified in the
# problem.

LIF := proc (F, z, P, n) local offset, G, phi:
    offset := -(coeff(expand(F), z, 0)-P):
    if (P in convert(offset, set)) then:
        return FAIL:
    fi:
    G := subs(P=P-offset, F):
    phi := -(G-P)/z:
    return (subs(P=z, phi))^n/n/z^(n-1):
end: