### Homework 3 - January 26, 2012 - Patrick Devlin ##
# I have no preferences about keeping my work private #
### Programs from class: ###
Help:=proc() : print ( `Di(L), IsConst(L), IsPol(L), GuessPol(L,x), GuessRat1(L,x,d1,d2), GuessRat(L,x), FF(x,k), BetterGuessPol(L,x), TimeTrial(L,x)`): end proc:
# Di(L): the sequence of differences of consecutive entries of L
# For example, Di([1,4,9]) should be [3,5]
Di:=proc(L) local i:
return [seq(L[i] - L[i-1], i=2..nops(L))]:
end proc:
#IsConst(L): is the sequence of numbers L a constant sequence
IsConst:=proc(L):
return evalb(nops(convert(L,set)) <= 1):
end proc:
#IsPol(L): this returns the smallest possible degree of a polynomial that fits the data
IsPol:=proc(L) local i, L1, foundPoly, tolerance:
if (nops(L)<1) then return FAIL: fi:
L1:=L:
tolerance:=nops(L):
for i from -1 to tolerance do
if (L1=[0$nops(L1)]) then return i: fi:
L1:=Di(L1):
od:
return FAIL:
end proc:
### Problem 1 ###
# Skipped because it is for novices #
### Problem 2 ###
#GuessPol(L,x): this returns a guess of a polynomial, p(x), that fits the first couple terms of L
GuessPol:=proc(L,x) local p, i, a, d:
d:=IsPol(L):
if(d=FAIL) then return FAIL: fi:
if(d<1) then return L[1]: fi:
p:=0:
for i from 0 to d do
p:=p+a[i] * x^i:
od:
p:=sort(eval(p,solve([seq(eval(p,x=i)=L[i],i=1..min(d+1,nops(L)))])),x): # This won't fail since it only fits the polynomial to the first d+1 points
for i from 1 to nops(L) do: # this tests the polynomial at all points in L to see if it in fact fits
if (eval(p,x=i) <> L[i]) then return FAIL: fi:
od:
return p:
end proc:
### Problem 3 ###
# takes a list L of data, an indeterminant x, and degrees d1 and d2, and tries to fit a rational
# polynomial p(x)/q(x) to the data in L with deg(p)=d1, deg(q)=d2. It returns FAIL on failure.
GuessRat1:=proc(L,x,d1,d2) local r,p, q, i, a, b, solns:
if(d1+d2+4 > nops(L)) then return FAIL: fi:
if( (d1 < 0) or (d2 < 0)) then return FAIL: fi:
p:=0:
q:=0:
for i from 0 to d1 do # form p(x)
p:= p+a[i] * x^i:
od:
for i from 0 to d2 do # form q(x)
q:=q+b[i] * x^i:
od:
try # try to solve equations for p and q (try statement for good measure)
solns:=solve( [seq(eval(p,x=i)=L[i]*eval(q,x=i),i=1..nops(L))] , {seq(a[i],i=0..d1),seq(b[i],i=0..d2)}):
catch:
return FAIL:
end try:
if((solns = NULL) or (simplify(eval(q,solns))=0)) then return FAIL: fi: # Checks to see if solns went ok
try
r:=simplify(eval(p/q,solns)): # tries to form what would be the output, r = p/q
catch:
return FAIL:
end try:
for i from 1 to nops(L) do
if(eval(r,x=i) <> L[i]) then return FAIL: fi: # checks the output against data for good measure
od:
return r:
end proc:
### Problem 4 ###
# GuessRat(L,x) takes a list L and indeterminant x and tries to spit out a rational polynomial
# r(x) = p(x)/q(x) that fits L with deg(p) + deg(q) minimal. If it can't find
# such a poly with deg(p) + deg(q) + 4 <= nops(L), it returns FAIL
GuessRat:=proc(L,x) local degSum, d1, r:
for degSum from 0 to nops(L) do # note, we need not worry about degSum getting too big since GuessRat1 does that already
for d1 from 0 to degSum do
r:=GuessRat1(L,x,d1,degSum-d1):
if (r <> FAIL) then return r: fi:
od:
od:
return FAIL:
end proc: ## Since GuessRat1 already simplifies the answer, I bet GuessRat1(L,x) will always either return fail or something unique anyway
### Problem 5 ###
# FF(x,k) returns the falling factorial of x for k nonnegative integer,
# FF(x,0) -> 1, FF(x,1) -> x, FF(x,2) -> x(x-1), FF(x,3) -> x(x-1)(x-2), etc.
FF:=proc(x,k) local p, i:
p:=1:
for i from 0 to k-1 do:
p:=p*(x-i):
od:
return p:
end proc:
### Problem 6 ###
# To prove P_k (x+1) - P_k (x) = k P_{k-1} (x), simply divide both sides by k! and it reduces to Pascal's identity for binomial coefficients.
# We may in fact assume x is a positive integer greater than k because we know two polynomials, f and g, are equal for
# all values of x if and only if they are equal for all sufficiently large integer inputs
# (This is because the polynomial f(x) - g(x) would have infinitely
# many zeros, implying it is the zero polynomial)
# Now let F(x) be any entire complex-valued function (any function that is
# analytic on all of C) [in particular, F(x) may be a polynomial].
# Then since F is entire, we know that, F(x) may be written uniquely as
# F(x) = a_0 + a_1 x^1 + a_2 x^2 + a_3 x^3 + ..., where a[0], a[1], ... are in C
#
# Let C[[x]] denote the ring of formal power series over C.
# Then C[[x]] may be viewed as a countably infinite vector space over C,
# with one basis B:={1, x, x^2, x^3, ...}.
# However, the set B':={1, x, x(x-1), x(x-1)(x-2), ...} = {P_0 (x), P_1 (x), P_2 (x), ...}
# is also a basis of C[[x]] over C. Thus, since the analytic function F(x)
# has a unique expansion in the basis B, it also has a unique expansion in B'.
#
# Therefore, the function F(x) may be uniquely written in the form
# F(x) = b_0 P_0 (x) + b_1 P_1 (x) + b_2 P_2 (x) + ..., and this
# series expansion also converges for all x in C.
#
# Now we can recover the coefficients b_0 by taking discrete derivatives.
# In particular, we have F(0) = b_0 + b_1 * 0 + b_2 * 0 + ... = b_0.
# Moreover, if Di is the forward difference operator, Di(F) := F(x+1) - F(x),
# and Di^k is the k^th iterated composition of Di, then we have
# Di^k (F)(x) = sum( b_n Di^k P_{n} (x), n=0..infinity)
# = sum( b_{n} n(n-1)(n-2)...(n-k+1) P_{n-k} (x), n=0..infinity)
# = sum( b_{n+k} (n+k)(n+k-1)...(n+1) P_{n} (x), n=0..infinity)
# Therefore, we have [Di^k (F) (0) ] / k! = b_{k}.
# (Note we may push the Di^k through the infinite summation because F(x) is entire)
#
# Therefore, we may recover the coefficients b_k by b_k = [Di^k (F)(0)] / k!,
# and in particular, if F(x) is a polynomial, this will recover the entire
# polynomial F(x)
### Problem 7 ###
## (See programs below)
# It seems like BetterGuessPol is faster than GuessPol on the surface level.
# However, if you require for instance that you expand the output of
# BetterGuessPol or you check the output against the input list
# (which GuessPol does unnecessarily, so it's "only fair")
# This time improvement is much less noticeable.
# BetterGuessPol(L,x) returns an expression for the polynomial of minimal degree that fits the data L.
BetterGuessPol:=proc(L,x) local i, p, L2,ff:
p:=0:
L2:=L:
ff:=1:
for i from 1 to nops(L) do
if (L2=[0$nops(L2)]) then return p: fi: # if L1 is all 0's, then we're done
p:=p+ff*L2[1]:
L2:=Di(L2):
ff:=ff*(x-i)/i: # note this is NOT ff:=ff*(x-i+1)/i as expected since we need to compensate
# for the fact that we can only access "L2[1]" and not "L2[0]" as our formula requires
od:
return FAIL:
end proc:
# TimeTrial(L,x) times the performance of GuessPol and BetterGuessPol on the input L
TimeTrial:=proc(L,x) local t0, timeG, timeBG,p,i:
t0:=time():
GuessPol(L,x):
timeG:=time()-t0:
t0:=time():
p:=BetterGuessPol(L,x):
sort(expand(p),x): # expand
for i from 1 to nops(L) do # check to see that the program did as it should have, which is "only fair" since GuessPol does
if(eval(p,x=i) <> L[i]) then print(`Error!`): fi:
od:
timeBG:=time()-t0:
printf(`\nGuessPol took %e seconds, and BetterGuessPol took %e seconds.\n`, timeG, timeBG):
end proc: