Part Frac Expan

# Partial Fraction Expansion for the HP-41

Overview

-Given a rational function   R(x) = P(x) / Q(x)        with  Q(x) = [ q1(x) ]µ1 ..................  [ qn(x) ]µn     and  gcd( qi , qj ) = 1  for all i # j ,
this program returns the partial fraction expansion:

R(x) = E(x) +  p1,1(x) / [ q1(x) ]µ1 + p1,2(x) / [ q1(x) ]µ1-1 + ........................ + p1,µ1(x) / q1(x)                   where  deg pi,k < deg qi
+ ..........................................................................................................................

+ pn,1(x) / [ qn(x) ]µn + pn,2(x) / [ qn(x) ]µn-1 + ........................ + pn,µn(x) / qn(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) / { [ qi(x) ]µi q(x) }   where  q(x) = ¶ j#i [ qj(x) ]µj            ( ¶ = product )

since A0 = q and A1 = qiµi  are relatively prime,  there are polynomials U and V that verify   A0 U + A1 V = p    ( Bezout's identity )

-Successive Euclidean divisions produce                                                           We also define:

A0  = A1 Q1 + A2         deg A2 < deg A1                                                       U0 = 1 , U1 = 0  and  Uk+1 = Uk-1 - Qk Uk
A1  = A2 Q2 + A3         deg A3 < deg A2
............................................................

Am-1 = Am Qm + Am+1    deg Am+1 < deg Am   and  Am+1 = 0

-The next to last remainder Am = gcd( q , qiµi )  is a constant polynomial since  q & qiµi are relatively prime.
-One of the polynomials U is then obtained by  U = p Um / Am
-Finally,  pi,1 , pi,2 , ............. , pi,µi  are the remainders  R1 R2 ..... Rµi  in the successive divisions by qi

U  = B1 qi + R1  ,   B1 = B2 qi + R2  ,   B2 = B3 qi + R3  ,  .......................

- V could be computed in a similar way ( V0 = 0 , V1 = 1  and  Vk+1 = Vk-1 - Qk Vk  ,  V = p Vm / Am  ),
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)
•  Rbb1 .... •  Ree1 = coeff of q1(x)  .......................  •  Rbbn .... •  Reen = coeff of qn(x)    ( denominators )
•  Rbb" = bb1.ee1 ...................  •  Ree" = bbn.een       ( control numbers )
•  Rbb' = µ1 ...........................  •  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

 01  LBL "PFE"   02  STO 04   03  STO 05   04  STO 07   05  STO 14   06  FRC   07   E3   08  *   09  1   10  +   11  STO 08   12  RCL 14   13  INT   14  -   15  ST+ 08   16  .1   17  %   18  +   19  ST- 05   20  ST+ 14   21  X<>Y   22  STO 06   23  RCL 05   24  STO 11   25  RCL 14   26  STO 10   27  STO 12   28  RCL 08   29  STO 13   30  LBL 01   31  RCL IND 07   32  STO 09   33  RCL IND 11   34  LBL 02   35  RCL 08   36  XEQ 04   37  DSE 09   38  X=0?   39  GTO 00   40  RCL IND 11   41  RCL Z   42  XEQ "PRO"   43  GTO 02   44  LBL 00   45  STO IND 12   46  FRC   47   E3   48  *   49  1   50  + 51  STO 08   52  ISG 07   53  ISG 11   54  ISG 12   55  GTO 01   56  RCL 08   57  RCL IND 10   58  LBL 03   59  ISG 10   60  X=0?   61  GTO 00   62  RCL IND 10   63  RCL Z   64  XEQ "PRO"   65  RCL 08   66  XEQ 04   67  GTO 03   68  LBL 04   69  STO T   70  SIGN   71  LBL 08   72  CLX   73  RCL IND Y   74  STO IND L   75  ISG L   76  CLX   77  ISG Y   78  GTO 08   79  R^   80  LASTX   81  STO Z   82  DSE X   83   E3   84  /   85  +   86  RTN   87  LBL 00   88  STO 07   89  FRC   90  RCL 06   91  FRC   92   E3   93  ST* Z   94  *   95  RCL 06   96  INT   97  -   98  X<>Y   99  STO 08 100  RCL 07 101  INT 102  - 103  X>Y? 104  GTO 00 105  RCL 06 106  RCL 07 107  XEQ "DIV" 108  GTO 04 109  LBL 00 110  RCL 06 111  0 112  LBL 04 113  STO 15 114  X<>Y 115  STO 06 116  SF 01 117  ISG 08 118  LBL 05 119  RCL 07 120  RCL 08 121  XEQ 04 122  RCL IND 14 123  XEQ "DIV" 124  STO 09 125  FRC 126   E3 127  * 128  STO 10 129  RCL 09 130  INT 131  - 132  ISG 10 133  CLX 134  RCL IND 14 135  FRC 136   E3 137  * 138  RCL IND 14 139  INT 140  - 141  RCL 10 142  1.001 143  * 144  STO 11 145  LASTX 146  + 147  STO 12 148  CLX 149  STO IND 12 150  SIGN 151  STO IND 11 152  RDN 153  X<=Y? 154  GTO 06 155  RCL 09 156  X<> IND 14 157  STO 09 158  CLX 159  STO IND 11 160  SIGN 161  STO IND 12 162  LBL 06 163  RCL 09 164  RCL IND 14 165  XEQ "DIV" 166  STO 01 167  CLX 168  RCL IND 14 169  STO 09 170  SIGN 171  - 172  ISG X 173  X=0? 174  GTO 00 175  XEQ "DEL" 176  STO IND 14 177  RCL IND X 178  X=0? 179  GTO 00 180  RCL 12 181  RCL 01 182  RCL 12 183  FRC 184   E3 185  * 186  1 187  + 188  XEQ "PRO" 189  RCL 11 190  X<>Y 191  ENTER^ 192  FRC 193   E3 194  * 195  1 196  + 197  XEQ "ADD" 198  X<> 12 199  RCL 10 200  XEQ 04 201  STO 11 202  RCL 12 203  RCL Z 204  XEQ 04 205  STO 12 206  GTO 06 207  LBL 00 208  RCL 12 209  RCL 09 210  ISG X 211  ACOS 212  1 213  - 214  XEQ "DIV" 215  RCL 10 216  XEQ 04 217  RCL 06 218  RCL Z 219  XEQ "PRO" 220  LBL 07 221  STO 01 222  FRC 223   E3 224  STO 00 225  * 226  RCL 01 227  INT 228  - 229  STO 02 230  RCL IND 05 231  FRC 232  RCL 00 233  * 234  RCL IND 05 235  INT 236  - 237  X<=Y? 238  GTO 00 239  RCL 08 240  STO Z 241  + 242  STO 03 243  RCL 00 244  / 245  + 246  X<> 01 247  INT 248  RCL 03 249  RCL 02 250  - 251  STO 03 252  RCL 00 253  / 254  + 255  RCL 02 256  1 257  ST- 03 258  + 259  RCL 00 260  X^2 261  / 262  + 263  REGMOVE 264  RCL 03 265  RCL 00 266  / 267  RCL 08 268  + 269  CLRGX 270  LBL 00 271  RCL 01 272  RCL IND 05 273  XEQ "DIV" 274  STO 11 275  CLX 276  RCL 13 277  XEQ 04 278  X<>Y 279  STO 13 280  RCL 11 281  DSE IND 04 282  GTO 07 283  ISG 04 284  ISG 05 285  ISG 14 286  GTO 05 287  CF 01 288  RCL 06 289  RCL 14 290  INT 291  RCL 13 292  DSE X 293   E3 294  / 295  + 296  RCL 15 297  END

( 452 bytes / SIZE max )

 STACK INPUTS OUTPUTS Z / bbb.eee3 Y bbb.eee bbb.eee2 X bbb'.eee' bbb.eee1

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.eee1     is the control number of E(x)                                         Exception:   if  E(x) = 0 ,  X-output = 0.
bbb.eee2     is the control number of the new numerators
bbb.eee3     is the control number of the expanded p(x)

Example1:       R(x) = P(x)/Q(x) = ( 6 x5 - 19 x4 +20 x3 - 7 x2 + 7 x + 10 ) / [ ( 2 x2 + x + 1 ).( x - 2 )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     p1,1(x)  = 2 x + 3
R35 = 4                      so     p2,1(x) = 4
R36 = 5                      so     p2,2(x) = 5

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

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

-SIZE 054 is enough for this example.

Example2:       R(x) = P(x)/Q(x) = x5 /( 3 x2 + 1 )2

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 x2 + 1 )2 - 2 (x/9) / ( 3 x2 + 1 )

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

-SIZE at least 039.

Another checking:

-With P(x) = ( 192 x11 + 552 x10 + 1427 x9 + 2825 x8 + 4694 x7 + 6247 x6 + 6745 x5 + 6172 x4 + 4071 x3 + 2329 x2 + 829 x + 125 )
and  Q(x) = ( 3 x + 1 )2 ( x2 + x + 1 )3 ( 2 x2 - x + 3 )2

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

-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 qj = 1    the numerators are constants
by groups of 2 numbers if  deg qj = 2   the numerators are polynomials of degree 1
by groups of 3 numbers if  deg qj = 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 qj(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  qi & qj are relatively prime for all i # j

-For instance:

1 / [ ( x3 + x + 1 )2 ( x5 + x + 1 )  ] = ( 7 x2/3 - 5 x/3 + 11/3 ) / ( x3 + x + 1 )2 + ( -11 x2 + 25 x/3 - 52/3 ) / ( x3 + x + 1 )
+ ( 11 x4 - 25 x3/3 + 19 x2/3 - 5 x + 44/3 ) / ( x5 + x + 1 )

-Lines 210 to 213:    ISG X   ACOS   1   -     will display a DATA ERROR message if the qi-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

 01  LBL "PFE"   02  STO 04   03  STO 05   04  STO 07   05  STO 14   06  DEG F   07  LASTX   08  STO 08   09  X<>Y   10  X+1   11  ST+ 08   12  .1   13  %   14  +   15  ST- 05   16  ST+ 14   17  R^   18  STO 06   19  RCL 05   20  STO 11   21  RCL 14   22  STO 10   23  STO 12   24  RCL 08   25  STO 13   26  LBL 01   27  RCL IND 07   28  STO 09   29  RCL IND 11   30  LBL 02   31  RCL 08   32  XEQ 04   33  DSE 09   34  X=0?   35  GTO 00   36  RCL IND 11   37  RCL Z   38  XEQ "PRO" 39  GTO 02   40  LBL 00   41  STO IND 12   42  DEG F   43  LASTX   44  STO 08   45  ISG 07   46  ISG 11   47  ISG 12   48  GTO 01   49  RCL 08   50  RCL IND 10   51  LBL 03   52  ISG 10   53  X=0?   54  GTO 00   55  RCL IND 10   56  RCL Z   57  XEQ "PRO"   58  RCL 08   59  XEQ 04   60  GTO 03   61  LBL 04   62  LCO   63  ENTER^   64  DEG F   65  X<> L   66  X<>Y   67  RTN   68  LBL 00   69  STO 07   70  DEG F   71  LASTX   72  STO 08   73  CLX   74  RCL 06   75  DEG F   76  XY   88  STO 06   89  SF 01   90  LBL 05   91  RCL 07   92  RCL 08   93  XEQ 04   94  RCL IND 14   95  XEQ "DIV"   96  STO 09   97  DEG F   98  LASTX   99  STO 10 100  CLX 101  RCL IND 14 102  DEG F 103  RCL 10 104  1.001 105  * 106  STO 11 107  LASTX 108  + 109  STO 12 110  CLX 111  STO IND 12 112  SIGN 113  STO IND 11 114  RDN 115  X<=Y? 116  GTO 06 117  RCL 09 118  X<> IND 14 119  STO 09 120  CLX 121  STO IND 11 122  SIGN 123  STO IND 12 124  LBL 06 125  RCL 09 126  RCL IND 14 127  XEQ "DIV" 128  STO 01 129  CLX 130  RCL IND 14 131  STO 09 132  SIGN 133  - 134  ISG X 135  X=0? 136  GTO 00 137  XEQ "DEL" 138  STO IND 14 139  RCL IND X 140  X=0? 141  GTO 00 142  RCL 12 143  DEG F 144  RCL 12 145  RCL 01 146  LASTX 147  XEQ "PRO" 148  RCL 11 149  X<>Y 150  ENTER^ 151  DEG F 152  X<> L 153  XEQ "ADD" 154  X<> 12 155  RCL 10 156  LCO 157  STO 11 158  DEG F 159  RCL 12 160  LASTX 161  LCO 162  STO 12 163  GTO 06 164  LBL 00 165  RCL 12 166  RCL 09 167  ISG X 168  ACOS 169  X-1 170  XEQ "DIV" 171  RCL 10 172  XEQ 04 173  RCL 06 174  RCL Z 175  XEQ "PRO" 176  LBL 07 177  STO 01 178  DEG F 179  STO 02 180  RCL IND 05 181  DEG F 182  X<=Y? 183  GTO 00 184  RCL 08 185  STO Z 186  + 187  STO 03 188  X/E3 189  + 190  X<> 01 191  RCL 03 192  RCL 02 193  - 194  LCO 195  INT 196  X-1 197  X/E3 198  RCL 08 199  + 200  CLRGX 201  LBL 00 202  RCL 01 203  RCL IND 05 204  XEQ "DIV" 205  STO 11 206  CLX 207  RCL 13 208  XEQ 04 209  X<>Y 210  STO 13 211  RCL 11 212  DSE IND 04 213  GTO 07 214  ISG 04 215  ISG 05 216  ISG 14 217  GTO 05 218  CF 01 219  RCL 06 220  RCL 14 221  INT 222  RCL 13 223  X-1 224  X/E3 225  + 226  RCL 15 227  END

( 376 bytes / SIZE max )

 STACK INPUTS OUTPUTS Z / bbb.eee3 Y bbb.eee bbb.eee2 X bbb'.eee' bbb.eee1

>>>>  The instructions are identical.