Partial Fraction Expansion for the HP-41

Overview
 

-Given a rational function   R(x) = P(x) / Q(x)        with  Q(x) = [ q₁(x) ]^µ1 ..................  [ q_n(x) ]^µn     and  gcd( q_i , q_j ) = 1  for all i # j ,
  this program returns the partial fraction expansion:

         R(x) = E(x) +  p_1,1(x) / [ q₁(x) ]^µ1 + p_1,2(x) / [ q₁(x) ]^µ1-1 + ........................ + p_1,µ1(x) / q₁(x)                   where  deg p_i,k < deg q_i
                            + ..........................................................................................................................

                            + p_n,1(x) / [ q_n(x) ]^µn + p_n,2(x) / [ q_n(x) ]^µn-1 + ........................ + p_n,µn(x) / q_n(x)

  E(x) is the quotient in the Euclidean division  P(x) = E(x) Q(x) + p(x)      p(x) is the remainder.

-The Euclidean algorithm is used to perform the decomposition as follows:
-We have to expand  p(x) / { [ q_i(x) ]^µi q(x) }   where  q(x) = ¶ _j#i [ q_j(x) ]^µj            ( ¶ = product )

  since A₀ = q and A₁ = q_i^µi  are relatively prime,  there are polynomials U and V that verify   A₀ U + A₁ V = p    ( Bezout's identity )

-Successive Euclidean divisions produce                                                           We also define:

   A₀  = A₁ Q₁ + A₂         deg A₂ < deg A₁                                                       U₀ = 1 , U₁ = 0  and  U_k+1 = U_k-1 - Q_k U_k
   A₁  = A₂ Q₂ + A₃         deg A₃ < deg A₂
   ............................................................

  A_m-1 = A_m Q_m + A_m+1    deg A_m+1 < deg A_m   and  A_m+1 = 0

-The next to last remainder A_m = gcd( q , q_i^µi )  is a constant polynomial since  q & q_i^µi are relatively prime.
-One of the polynomials U is then obtained by  U = p U_m / A_m
-Finally,  p_i,1 , p_i,2 , ............. , p_i,µi  are the remainders  R₁ R₂ ..... R_µi  in the successive divisions by q_i

   U  = B₁ q_i + R₁  ,   B₁ = B₂ q_i + R₂  ,   B₂ = B₃ q_i + R₃  ,  .......................

- V could be computed in a similar way ( V₀ = 0 , V₁ = 1  and  V_k+1 = V_k-1 - Q_k V_k  ,  V = p V_m / A_m  ),
   but it's not useful here, so "PFE" does not calculate this polynomial.

-The method is repeated for i = 1 , 2 , ......... , n
 

Program Listing#1  ( without M-Code )
 
 

Data Registers:              R00 thru R15: temp                             ( All the "• Registers" are to be initialized before executing "PFE" )

                                      •  Rbb .................   •  Ree = coefficients of the numerator P(x)
                                      •  Rbb₁ .... •  Ree₁ = coeff of q₁(x)  .......................  •  Rbb_n .... •  Ree_n = coeff of q_n(x)    ( denominators )
                                      •  Rbb" = bb₁.ee₁ ...................  •  Ree" = bb_n.ee_n       ( control numbers )
                                      •  Rbb' = µ₁ ...........................  •  Ree' = µ_n                     ( exponents )

                                         Ree'+1  ...... temp                  --- This program may use a lot of data registers ---

Flag:    F01 is cleared at the end
Subroutines:  "PRO"  "DIV"  ADD"  "DEL"  ( cf "Polynomials for the HP-41" )
 

-Line 286 is a three-byte GTO 05
 
 

 
    ( 452 bytes / SIZE max )
 
 

   where

     bbb.eee   is the control number of the numerator   P(x)
     bbb'.eee'  is the control number of the exponents  µ_j

 >>>>    eee' +1 = the first free-register - All the registers Rnn  for  nn >  eee'   must be free.                                          bbb > 15

    bbb.eee₁     is the control number of E(x)                                         Exception:   if  E(x) = 0 ,  X-output = 0.
    bbb.eee₂     is the control number of the new numerators
    bbb.eee₃     is the control number of the expanded p(x)

Example1:       R(x) = P(x)/Q(x) = ( 6 x⁵ - 19 x⁴ +20 x³ - 7 x² + 7 x + 10 ) / [ ( 2 x² + x + 1 ).( x - 2 )² ]
 

  Store   [ 6 -19  20 -7  7  10 ]   into   R16  R17  R18  R19  R20  R21      control number  = 16.021

    and      [ 2  1  1 ]  into  R22  R23 R24  &  [ 1  -2 ]  into  R25 R26   control numbers = 22.024  &  25.026

         store these 2 control numbers into R27 & R28   ( 22.024  STO 27 ,  25.026  STO 28 )

  and store the exponents  1 and 2  into  R29 & R30   ( control number = 29.030 )

     16.021   ENTER^
     29.030   XEQ "PFE"  >>>>   16.017                  --- Execution time = 1mn55s ---
                                       RDN    33.036
                                       RDN    18.021

   R16 = 3 ,  R17 = 1       so      E(x) = 3 x + 1

         R33 = 2   R34 = 3      so     p_1,1(x)  = 2 x + 3
         R35 = 4                      so     p_2,1(x) = 4
         R36 = 5                      so     p_2,2(x) = 5

-Finally,     R(x) = 3 x + 1 + ( 2 x + 3 ) / ( 2 x² + x + 1 ) + 4 / ( x - 2 )² + 5 / ( x - 2 )

-We have also  [ R18  R19  R20  R21 ]  = [ 12  -12  -5  6 ]    whence    R(x) = 3 x + 1 + ( 12 x³ - 12 x² - 5 x + 6 ) / [ ( 2 x² + x + 1 ).( x - 2 )² ]

-SIZE 054 is enough for this example.
 

Example2:       R(x) = P(x)/Q(x) = x⁵ /( 3 x² + 1 )²
 

   Store  [ 1  0  0  0  0  0 ]  into  R16  R17  R18  R19  R20  R21   control number = 16.021

        and    [ 3  0  1 ]  into  R22  R23  R24   control number = 22.024  store it into R25  and store the exponent 2 into R26

  Then,

       16.021   ENTER^
       26.026   XEQ "PFE"  >>>>   16.017                  --- Execution time = 41s ---
                                         RDN    28.031
                                         RDN    18.021

   R16 = 1/9 ,  R17 = 0       therefore    E(x) = x/9 + 0

   [ R28  R29  R30  R31 ] = [ 1/9  0  -2/9  0 ]

-Whence:     R(x) = x/9 + (x/9) / ( 3 x² + 1 )² - 2 (x/9) / ( 3 x² + 1 )

   [ R18  R19  R20  R21 ] = [ -2/3  0  -1/9  0 ]    thus    R(x) = x/9 + ( -2x³/3 - x/9 ) / ( 3 x² + 1 )²

-SIZE at least 039.
 

Another checking:
 

-With P(x) = ( 192 x¹¹ + 552 x¹⁰ + 1427 x⁹ + 2825 x⁸ + 4694 x⁷ + 6247 x⁶ + 6745 x⁵ + 6172 x⁴ + 4071 x³ + 2329 x² + 829 x + 125 )
  and  Q(x) = ( 3 x + 1 )² ( x² + x + 1 )³ ( 2 x² - x + 3 )²

  It yields   R(x) = P(x) / Q(x) = 2 / ( 3 x + 1 )² + 1 / ( 3 x + 1 ) + ( 3 x + 4 ) / ( x² + x + 1 )³ + ( 5 x + 6 ) / ( x² + x + 1 ) + ( 7 x + 8 ) / ( 2 x² - x + 3 )²

-If you store the coefficients and control numbers as above into contiguous registers from R16,

   16.027  ENTER^
   39.041  XEQ "PFE"  >>>>    0                           --- execution time = 8mn58s ---    SIZE at least 104
                                     RDN   45.056
                                     RDN   16.027
 

      R45 = 2                           R46 ~ 1
      R47 ~ 3    R48 ~  4         R49 = 0   R50 ~ 0         R51 ~ 5   R52 ~ 6
      R53 ~ 7    R54 ~  8         R55 ~ 0   R56 ~ 0

-The contents of these registers are to be read

        by groups of 1  number if  deg q_j = 1    the numerators are constants
        by groups of 2 numbers if  deg q_j = 2   the numerators are polynomials of degree 1
        by groups of 3 numbers if  deg q_j = 3   the numerators are polynomials of degree 2 ,  and so on ....

-Roundoff-errors are more important in this example:  for some coefficients,  only 5 or 6 decimals are exact.
 

Notes:

-Usually,  the denominators q_j(x) are polynomials of degree 1 or polynomials of degree 2 with a negative discriminant.
-However, this program will also work to get a - partial - decomposition in other cases, provided  q_i & q_j are relatively prime for all i # j

-For instance:

      1 / [ ( x³ + x + 1 )² ( x⁵ + x + 1 )  ] = ( 7 x²/3 - 5 x/3 + 11/3 ) / ( x³ + x + 1 )² + ( -11 x² + 25 x/3 - 52/3 ) / ( x³ + x + 1 )
                                                            + ( 11 x⁴ - 25 x³/3 + 19 x²/3 - 5 x + 44/3 ) / ( x⁵ + x + 1 )
 

-Lines 210 to 213:    ISG X   ACOS   1   -     will display a DATA ERROR message if the q_i-polynomials are not relatively prime.
-Otherwise, they could be deleted.
 

Program Listing#2
 

-If you can use the M-Code routines   LCO   ( cf "Miscellaneous Short Routines" )
  and  DEG F   X/E3   X+1   X-1  ( cf "A few M-Code Routines" ),  the program will be both shorter & faster.

-Line 217 is a three-byte GTO 05
 
 

 
    ( 376 bytes / SIZE max )
 
 

 
 >>>>  The instructions are identical.