### 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: