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 X<Y? |
77 GTO 00
78 RCL 06 79 RCL 07 80 XEQ "DIV" 81 GTO 04 82 LBL 00 83 RCL 06 84 0 85 LBL 04 86 STO 15 87 X<>Y 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.