%% Below, I've used Maple, treating polynomials as /expressions/.  I
%% could, instead, have treated polynomials as functions.


{{{Maple}}} Wed Feb 10 21:56:13 EST 2010

    |\^/|     Maple 10 (SUN SPARC SOLARIS)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2005
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.


%% Here is our numerator.

> Numxm :=

          10    9      8      7       6       5      4       3       2
    29 + x   + x  + 7 x  + 2 x  + 13 x  - 17 x  + 5 x  - 34 x  + 14 x  - 12 x ;


%% And here is our denominator.

                      8    7      6      5    4       3      2
>         Dxm := 8 + x  + x  + 4 x  - 2 x  + x  - 11 x  + 2 x  - 4 x ;

%% Let's compute the quotient and remainder polynomials:

> Quox := quo(Numxm, Dxm, x, 'Remx');   Remx  := Remx ;

                                         2
                                Quox := x  + 3

                                      7      3
                             Remx := x  + 3 x  + 5

%% Since I know how to integrate a poly, Quox, we'll now deal with
%% integrating the rational fnc
%% 		Remx / Dxm .

%% The form I want to write Remx/Dxm  in is:

> parfracexp := A1/(x-1) +  A2/(x-1)^2 + ((B1*x)+C1)/(x^2  + x + 2)  + ((B2*x)+C2)/(x^2  + x + 2)^2  + ((B3*x)+C3)/(x^2  + x + 2)^3 ;


                 A1        A2      B1 x + C1      B2 x + C2       B3 x + C3
  parfracexp := ----- + -------- + ---------- + ------------- + -------------
                x - 1          2    2             2         2     2         3
                        (x - 1)    x  + x + 2   (x  + x + 2)    (x  + x + 2)


%% So we have a set, {A1,A2, B1,C1, B2,C2, B3,C3}, of eight unknown numbers.
%%
%%   In order to solve for them, note that the difference
%% 	[parfracexp * Dxm] - Remx
%% must be the zero-poly; which is what I called "zip" in class and in the
%% "Primer on Polynomials".

%%   So we will simply set each coeff of   [parfracexp*Dxm] - Remx   equal to
%% zero, then solve the resulting system of eight eqns; this will give us
%% values for {A1,A2, B1,C1, B2,C2, B3,C3}.


> Diffx := collect(simplify((parfracexp * Dxm) - Remx), x);

                        7                     6                              5
Diffx := (B1 + A1 - 1) x  + (A2 + 2 A1 + C1) x  + (2 B1 + 6 A1 + 3 A2 + B2) x

                                              4
     + (4 A1 - 4 B1 + 2 C1 + C2 + 9 A2 - B2) x

                                                       3
     + (-4 C1 + 13 A2 + B3 - 3 + B1 + 5 A1 - C2 + B2) x

                                                            2
     + (-6 A1 + C3 + 18 A2 - 4 B1 + C1 - 3 B2 - 2 B3 + C2) x

     + (-3 C2 - 4 C1 + 12 A2 + 2 B2 - 2 C3 + B3 - 4 A1 + 4 B1) x - 8 A1 + 4 C1

     + 8 A2 - 5 + C3 + 2 C2

%% The following step is a mathematically-unnecessary step... but, for some
%% reason, the Maple `collect' cmd does not collect the constant term as a
%% multiplier of x^0,  So, below, I multiply Diffx by x, and collect, so that
%% NOW, the eight coeffs of x^8, x^7, ..., x^1 are the coeffs I want.

> DiffxTimesx := collect(simplify(x*Diffx), x) ;

                              8                     7
DiffxTimesx := (B1 + A1 - 1) x  + (A2 + 2 A1 + C1) x

                                  6                                          5
     + (2 B1 + 6 A1 + 3 A2 + B2) x  + (4 A1 - 4 B1 + 2 C1 + C2 + 9 A2 - B2) x

                                                       4
     + (-4 C1 + 13 A2 + B3 - 3 + B1 + 5 A1 - C2 + B2) x

                                                            3
     + (-6 A1 + C3 + 18 A2 - 4 B1 + C1 - 3 B2 - 2 B3 + C2) x

                                                                2
     + (-3 C2 - 4 C1 + 12 A2 + 2 B2 - 2 C3 + B3 - 4 A1 + 4 B1) x

     + (-8 A1 + 4 C1 + 8 A2 - 5 + C3 + 2 C2) x

%% So there, above, are the eight coeffs I want to set to zero and then solve
%% for the eight unknows.


%%   The following Maple snippet will extract the coeffs, create eqns setting
%% each coeff equal to zero, and put them in a Maple "set", which is indicated by
%% the braces.

> Eqns := { seq( op(1, op(n,DiffxTimesx)) = 0, n=1..8 ) } ;

Eqns := {B1 + A1 - 1 = 0, A2 + 2 A1 + C1 = 0, 2 B1 + 6 A1 + 3 A2 + B2 = 0,

    4 A1 - 4 B1 + 2 C1 + C2 + 9 A2 - B2 = 0,

    -4 C1 + 13 A2 + B3 - 3 + B1 + 5 A1 - C2 + B2 = 0,

    -6 A1 + C3 + 18 A2 - 4 B1 + C1 - 3 B2 - 2 B3 + C2 = 0,

    -3 C2 - 4 C1 + 12 A2 + 2 B2 - 2 C3 + B3 - 4 A1 + 4 B1 = 0,

    -8 A1 + 4 C1 + 8 A2 - 5 + C3 + 2 C2 = 0}

%% Now we solve this system of linear eqns.  We have eight equations in eight
%% unknowns, and it will have a unique soln.

> solve( Eqns, {A1,A2, B1,C1, B2,C2, B3,C3} );


       273                  59                           -69       -1
{ B1 = ---, A2 = 9/64, B3 = --, C3 = 9/8, C2 = 9/8, B2 = ---, C1 = ---,
       256                  16                           32        128

         -17
    A1 = ---   }
         256


%% As a partial check, let's ask Maple to directly compute the partial-fraction
%% decomposition of Remx/Dxm.


> convert( Remx / Dxm , parfrac, x) ;


   18 + 59 x          -2 + 273 x           9            17
---------------- + ---------------- + ----------- - -----------
     2         3         2                      2   256 (x - 1)
16 (x  + x + 2)    256 (x  + x + 2)   64 (x - 1)

          36 - 69 x
     + ----------------
            2         2
       32 (x  + x + 2)


%%%%%%% Let's see if our `convert' produced the same answer %%%%%%%%%%%

  The immediately-above expression indeed gives the same values as what
solving "by hand" produced.  F'irinstance, the expression

              17
        - -----------
          256 (x - 1)

from the convert(...) can be rewritten as

            -17 / 256
           ----------       .
            (x - 1)

And indeed, the A1 that we solved for, up above, is


         -17
    A1 = ---  .
         256

%% Pretty Nifty, eh? %%