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