#Nathan Fox
#Homework 6
#I give permission for this work to be posted online

#Read procedures from class
read(`C6.txt`):

#Linear Algebra
with(LinearAlgebra):

#Help procedure
Help:=proc():
    print(` H(m,r,a,n) , h(m,r,a,n) `):
    print(` GuessRatioRatio(a,n,r,d:=20) `):
    print(` ProveConj(Exp, a, n, r) , MyGuessRF1(L,d1,d2,n) `):
    print(` MyGuessRF(L, n) , GuessDoubleRF1(L,d1,d2,m,n) `):
    print(` GuessDoubleRF(L,m,n) `):
end:

##PROBLEM 1##

#Generalized Hilbert Matrix
#m=dimension
#r=offset
#a=expression in n
#n=variable
H:=proc(m, r, a, n) local i, j:
    return [seq([seq(subs(n=r+i+j-1, a), j=1..m)], i=1..m)]:
end:

#Determinant of H(m,r,a,n)
h:=proc(m, r, a, n)
    return Determinant(H(m, r, a, n)):
end:

#Guess Ratio Ratio
#inputs the expression a,
#the discrete variables n and r,
#and outputs a guessed expression,
#as a rational function of n, for 
#(h(n+2,r)/h(n+1,r))/(h(n+1,r)/h(n,r))
#
#NOTE: returns fail if can't find a rational expression of low
#enough degree
GuessRatioRatio:=proc(a,n,r,d:=20) local m, L:
    L:=Ratio1(Ratio1([seq(h(m, r, a, n), m=1..(d+2))])):
    return factor(GuessRF(L, n)):
end:

##PROBLEM 2##

#Proves (rigorously!, by induction on n), that the output of
#GuessRatioRatio(a,n,r) is indeed correct
#using recurrence
#h(n,r)=(h(n-1,r)*h(n-1,r+2)-h(n-1,r+1)^2)/h(n-2,r+2)
#Returns true if proved, false if not
ProveConj:=proc(Exp, a, n, r) local L, m, N, R, i, Dodgson, Dex:
    Dex:=proc(NN, RR)
        return subs({n=NN, r=RR}, Exp):
    end:
    Dodgson:=proc(NN, RR)
        return (Dex(NN-1, RR) * Dex(NN-1, RR+2)
            - Dex(NN-1, RR+1)^2) / Dex(NN-2, RR+2):
    end:
    L:=Ratio1(Ratio1([seq(h(m, r, a, n), m=1..6)])):
    #Base cases
    if {evalb(subs(n=i, Exp) <> L[i])} <> {true} then
        print(`Please check your base cases and try again.`);
        return false:
    fi:
    #Induction
    return evalb(simplify((Dex(N+2, R)/Dex(N+1, R))/
        (Dex(N+1, R)/Dex(N, R)))  = simplify((Dodgson(N+2, R)/
        Dodgson(N+1, R))/(Dodgson(N+1, R)/Dodgson(N, R)))):
end:
#Note: This doesn't work, and I don't see how to write Dodgson's
#rule in terms of the ratio ratios, which is what would be required
#to fix it

##PROBLEM 3##

#I cannot do this problem in full, as I cannot get problem 2
#to work.  I get  (n + 2) (n + 1) (n + r + 2) (n + r + 1)
#                  ------------------------------------------
#                                                           2
#                  (2 n + 4 + r) (2 + 2 n + r) (2 n + 3 + r)
#as the rational function, though.

##PROBLEM 4##

#GuessRF1(L, d1, d2, n): inputs a list of numbers (or expressions)
#L, non-neg. integers d1 and d2, and a (discrete) variable name
#n and tries to guess a rational function of n of top degree
#d1 and bottom degree d2, such that L[i]-R(i)=0 for all i from
#1 to nops(L)
MyGuessRF1:=proc(L,d1,d2,n) local R, a, b, i, j, eq, var:
    if nops(L) <= 2*max(d1, d2)+6 then
        #print(`Make list bigger`):
        return FAIL:
    fi:

    R:=(n^d1+add(a[i]*n^i, i=0..d1-1))/add(b[j]*n^j, j=0..d2):
    var:={seq(a[i], i=0..d1-1), seq(b[j], j=0..d2)}:
    eq:={seq(numer(L[i]-subs(n=i,R))=0, i=1..nops(L))}:
    var:=solve(eq, var):
    if var = NULL then
        return FAIL:
    fi:
    R:=subs(var, R):
    if {seq( numer(L[i]-subs(n=i,R)), i=1..nops(L))}<>{0} then
        #print(`It worked up to `, nops(L), ` terms, but that's life`):
        return FAIL:
    fi:
    return R:
end:

#MyGuessRF(L,n): inputs a list of numbers (or expressions)
#L, and a (discrete) variable name
#n and tries to guess a rational function of n of degree
#<=(nops(L)-6)/2 such that L[i]-R(i)=0 for all i from 1 to nops(L)
MyGuessRF:=proc(L, n) local d1, d2, guess:
    for d1 from 0 to (nops(L)-6)/2 do
        for d2 from 0 to (nops(L)-6)/2 do
            guess:=MyGuessRF1(L, d1, d2, n):
            if guess <> FAIL then
                return guess:
            fi:
        od:
    od:
    return FAIL:
end:

##PROBLEM 5##

#GuessDoubleRF1(L, d1, d2, m, n): inputs a list of numbers
#(or expressions)
#L, non-neg. integers d1 and d2, and (discrete) variable names
#m and n and tries to guess a rational function of m and n of
#top degree d1 and bottom degree d2, such that L[i]-R(i)=0 for
#all i from 1 to nops(L)
GuessDoubleRF1:=proc(L,d1,d2,m,n)
local R, a, b, i, j, eq, var, ndeg:
    if nops(L) <= 2*max(d1, d2)+6 then
        return FAIL:
    fi:

    #Need to try all top degree leading coefficients, since
    #some might be 0
    for ndeg from 0 to d1 do
        R:=(m^(d1-ndeg)*n^ndeg+add(a[i][d1]*m^(d1-i)*n^i,
            i=ndeg+1..d1)+add(add(a[i][j]*m^(j-i)*n^i, i=0..j),
            j=0..d1-1)) / add(add(b[i][j]*m^(j-i)*n^i, i=0..j),
            j=0..d2):
        var:={seq(a[i][d1], i=ndeg+1..d1), seq(seq(a[i][j],
            i=0..j), j=0..d1-1), seq(seq(b[i][j], i=0..j),
            j=0..d2)}:
        eq:={seq(seq(numer(L[i][j]-subs({m=i,n=j},R))=0,
            j=1..nops(L[i])), i=1..nops(L))}:
        var:=solve(eq, var):
        if var = NULL then
            next:
        fi:
        R:=subs(var, R):
        return R:
    end:
    return FAIL:
end:

#GuessDoubleRF(L,m,n): inputs a list of numbers (or expressions)
#L, and (discrete) variable names
#m and n and tries to guess a rational function of m and n of
#degree<=(nops(L)-6)/2 such that L[i]-R(i)=0 for all i from 1 to nops(L)
GuessDoubleRF:=proc(L, m, n) local d1, d2, guess:
    for d1 from 0 to (nops(L)-6)/2 do
        for d2 from 0 to (nops(L)-6)/2 do
            guess:=GuessDoubleRF1(L, d1, d2, m, n):
            if guess <> FAIL then
                return guess:
            fi:
        od:
    od:
    return FAIL:
end: