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

# Nathaniel Shar
# HW 9
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

Help := proc(): print(`GuessRec1(L, order), GuessRec(L), RecToSeq(Seq, n), T(n, x), GFToSeq(f,q,m), SeqToGF(L, q)`): end:

# Stuff from Class 6 (somewhat modified):


# Guesses homogeneous linear recurrence of the specified order with
# constant coefficients for the sequence L.
# L(n) + c_1L(n-1) + ... + c_rL(n-r) = 0.
# Output:
# [[c_1, c_2, ..., c_r], InitialConditions]
GuessRec1 := proc(L, order) local eqns, c, vars:
    if (nops(L) < 2*order+4) then return FAIL: fi:
    eqns := [seq(L[n] + sum(c[i]*L[n-i], i=1..order), n=order+1..nops(L))]:
    vars := solve(eqns, {seq(c[i], i=1..order)}):
    if (vars = NULL) then return FAIL: fi:
    [subs(vars, [seq(c[i], i=1..order)]), L[1..order]]:
end:


# Guesses homogeneous linear recurrence of the smallest possible order with
# constant coefficients for the sequence L.
# L(n) + c_1L(n-1) + ... + c_rL(n-r) = 0.
# Output:
# [[c_1, c_2, ..., c_r], InitialConditions]


GuessRec := proc(L) local i, answer:
    for i from 1 to floor((nops(L)-4)/2) do:
        answer := GuessRec1(L, i):
        if answer <> FAIL then: break: fi:
    od:
    return answer:
end:

# Stuff from HW 6:

RecToSeq := proc(Seq, n) local ini, co, out, i, j:
    ini := Seq[2]:
    co := Seq[1]:
    out := ini:
    for i from nops(ini)+1 to n do:
        out := [op(out), add(-co[j]*out[i-j], j=1..nops(co))]:
    od:
    return out:
end:


#############
# Problem 1 #
#############

T := proc(n, x) option remember:
    if n = 1 then:
        return 1:
    elif n = 2 then:
        return 4*x^2:
    else:
        (-3 + 4*x^2 + T(n-2, x)*T(n-1, x) - 12*T(n-1, x)*x^2 - 8*T(n-2, x)*x^2 + 4*T(n-2, x)^2*x^2 - 16*T(n-2, x)*T(n-1, x)*x^4 + 8*T(n-2, x)*T(n-1, x)*x^2 + T(n-1, x)^2 + T(n-2, x) + 2*T(n-1, x))/(-T(n-2, x) - T(n-1, x) - 1);
    fi:
end:


#############

# From HW 5:

GFToSeq := proc(f,q,m):
    [seq(coeff(taylor(f, q, m+1), q, i), i=0..m)]:
end:

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

SeqToGF := proc(L, q) local C, den, num, L2, term:
    C := GuessRec(L):
    num := 0:
    if C = FAIL then return C: fi:
    den := 1 + add(C[1][i]*q^i, i=1..nops(C[1])):
    L2 := convert(expand(add(L[i]*q^(i-1), i=1..nops(C[1])) * den), list):
    for term in L2 do:
        if degree(term, q) < nops(C[1]) then:
            num := num + term
        fi:
    od:
    return num/den;
end:

#######################
# It appears to work...