#!/usr/local/bin/maple

# Nathaniel Shar
# HW 3
# Experimental Mathematics

Help := proc():
    print(`IsConst(L); Di(L); IsPoly(L); GuessPoly(L,x), GuessRat1(L, x, d1, d2), GuessRat(L, x), FF(x, k), BetterGuessPoly(L, x)`);
end:


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

# From class 3: checks if the list L is constant

IsConst := proc(L):
    evalb(nops(convert(L, set)) <= 1);
end:

# From class 3: returns the list of differences.

Di := proc(L):
    [seq(L[i] - L[i-1], i=2..nops(L))];
end:

# From class 3: checks if the list L appears to be polynomial.

IsPoly := proc(L) local i, list:
    if IsConst(L) then return 0: end:
    list := L:
    for i from 1 to nops(L)-5 do:
        if IsConst(Di(list)) then:
            return i:
        else:
            list := Di(list):
        fi:
    od:
    return FAIL:
end:

# GuessPoly: This version only forms the first d+1 equations. Please
# note that the call to IsPoly ensures that the sequence is in fact a
# polynomial of degree d, so there is no need to check this a second time by
# plugging the other values into the polynomial.
# However, it is done anyway here because it is specified in the assignment.

GuessPoly := proc(L,x) local d, eqns, P, i, a, guess:
    d := IsPoly(L):
    if d = 0 then return L[1]: end:
    if d = FAIL then return FAIL: end:
    P := add(a[i]*x^i, i=0..d);
    eqns := [seq(subs(x=i, P) = L[i], i=1..d+1)]:
    guess := sort(eval(P, solve(eqns))):
    for i from d+2 to nops(L) do:
        if subs(x=i, guess) <> L[i] then return FAIL: fi:
    od:
    return guess:
end:

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

# GuessRat1: Tries to find a rational function R such that R(i) = L[i]
# such that the degree of the numerator is d1, the degree of the
# denominator is d2, and the denominator is monic.

GuessRat1 := proc(L, x, d1, d2) local eqns, solns, R, a, b:
    if nops(L) <= d1+d2+4 then return FAIL: fi:
    R := add(a[i]*x^i, i=0..d1)/(add(b[i]*x^i, i=0..d2-1) + x^d2):
    eqns := [seq(subs(x=i, R) = L[i], i=1..nops(L))]:
    solns := solve(eqns):
    if solns = NULL then return FAIL: fi:
    eval(R, solns):
end:

#############
# Problem 4 #
#############

# GuessRat: Tries to find a rational function  R of smallest degree
# such that R(i) = L[i] and the denominator of R is monic.

# Warning: takes a lot longer (asymptotically) than GuessPoly!

GuessRat := proc(L, x) local sum, d1, guess:
    if nops(L) <= 3 then return FAIL: fi:
    for sum from 0 to nops(L)-4 do:
        for d1 from 0 to sum do:
            guess := GuessRat1(L, x, d1, sum-d1):
            if guess <> FAIL then return guess: fi:
        od:
    od:
    FAIL:
end:

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

# Falling factorial function: returns x to the k falling.
# k must be nonnegative.

FF := proc(x, k):
    product(x-i, i=0..k-1):
end:

#############
# Problem 6 #
#############

# First note that P_k(x+1) = (x+1)P_{k-1}(x). Meanwhile, P_k(x) =
# (x-k+1)P_{k-1}(x). Subtracting, we get
#            P_k(x+1) - P_k(x) = kP_{k-1}(x),                          (1)
# as desired.

# Given a polynomial P(x), we can therefore write
#  P(x) = a_0 + a_1P_1(x)/1! + a_2P_2(x)/2! + ... + a_nP_n(x)/n!,      (2)
# where a_i = \Delta^i P(0). To see this,  we first show that any
# polynomial of degree d can be written in the form
# c_0 + c_1P_1(x) + c_2P_2(x) + ... + c_dP_d(x).
# The proof is by induction on d. This is trivial when d = 0 or d = 1.
# For d >= 2, first note that
#            P_d(x) = x^d + k_{d-1}x^{d-1} + ... + k_0
# for some constants k_i. Now, by the induction hypothesis, let
# k_{d-1}x^{d-1} + ... + k_1x + k_0 be written in the form
# c_0 + c_1P_1(x) + ... + c_{n-1}P_{n-1}(x). Then
# x^d = - (c_0 + c_1P_1(x) + ... + c_{n-1}P_{n-1}(x)) + P_d(x).
# Thus we have written x^d in the desired form, and by the induction
# hypothesis we can also write x^{d-i} in the same form, so we can
# write any polynomial of degree d.

# Now, put P(x) = c_0 + c_1P_1(x) + ... + c_dP_d(x), and take the kth
# difference of P. By (1),
#   \Delta^k P(x) = sum_{i=0}^{d-k} c_{k+i}P_i(k+i)P_i(x)
# In particular, the constant term is k!c_{k+i}, so
#                \Delta^i P(0) = k!c_k.
# Solving for c_k, we get c_k = \Delta^i P(0)/k! = a_k/k!, and this
# completes the proof of (2).

# Of course a similar formula can be written centered at (x-1).

#############
# Problem 7 #
#############

# Guess polynomials using the calculus of finite differences.

BetterGuessPoly := proc(L, x) local i, d, coeffs, diff:
    d := IsPoly(L):
    if d = FAIL then return FAIL: fi:
    diff := L[1..d+1]:
    coeffs := []:
    for i from 0 to d do:
        coeffs := [op(coeffs), diff[1]/factorial(i)]:
        diff := Di(diff):
    od:
    return sort(expand(add(coeffs[i]*FF((x-1), i-1), i=1..d+1))):
end:

# Timing results show that for polynomials of large order,
# BetterGuessPoly is not actually better. For example,
# "GuessPoly([seq(i^361 + 17*i^71 + i^59 + 3, i=1..1000)], x)"
# took 18.6 seconds, while
# "BetterGuessPoly([seq(i^361 + 17*i^71 + i^59 + 3, i=1..1000)], x)"
# took 100.5 seconds.
# HOWEVER, almost all of this time was spent in the last step
# (expanding the falling factorials).