Multivariate Polynomials for the HP-41
Overview
1°) Two Variables
1-a) Evaluation of Polynomial Functions
1-b) Degree of a Polynomial
1-c) Product of 2 Polynomials
2°) N Variables
2-a) Evaluation of Polynomial Functions
2-b) Product of 2 Polynomials
2-c) Deleting the Zeros
2-d) Simplifications
2-e) Partial Derivatives
1°) Two Variables
-In this paragraph, the coefficients a , b , c , d , e , f , g , h , i , j , .... of a polynomial p(x,y) are to be stored into consecutive registers as follows:
p(x,y) = a + b x + c y + c x2 + d x.y + e y2 + f x3 + g x2y + h x.y2 + i y3 + j x4 + k x3y + l x2y2 + m x.y3 + n y4 +......
-Thus, we start with the constant term, then the monomials of degree
1 , the monomials of degree 2 , ... etc ...
-The polynomial is identified by the control number of its block of
registers: bbb.eee that begins with Rbb and ends with
Ree.
-All the coefficients must be stored ( including the zeros )
-So, we use:
1 register for a polynomial of degree
0
3 registers for a polynomial of degree
1
6 registers for a polynomial of degree
2
10 registers for a polynomial of degree 3
and, more generally, (p+1)(p+2)/2 registers for a polynomial
of degree p
-Inversely, the degree of the polynomial = (1/2)(-3+(8n+1)1/2)
where n is the number of coefficients in the block of registers ( cf §1-b)
).
a) Evaluation of Polynomial
Functions
Data Registers: R00 = bbb.eee ..... 1+eee.eee ( Registers Rbb thru Ree are to be initialized before executing "PXY" )
• Rbb ...... • Ree = coefficients of the polynomial
Flags: /
Subroutines: /
-Synthetic registers M N O P may be replaced by R01 R02
R03 R04
-In this case, delete line 42 and replace line 02 by 0
STO 04 RDN
| 01 LBL "PXY"
02 CLA 03 STO M 04 RDN 05 STO N 06 X<>Y 07 STO 00 08 RCL IND X 09 STO O |
10 ISG 00
11 FS? 30 12 GTO 03 13 LBL 01 14 RCL P 15 INT 16 E3 17 / 18 STO P |
19 SIGN
20 RCL IND 00 21 LBL 02 22 ISG 00 23 RCL M 24 * 25 X<>Y 26 RCL N 27 * |
28 STO Z
29 RCL IND 00 30 * 31 + 32 ISG P 33 GTO 02 34 ST+ O 35 ISG 00 36 GTO 01 |
37 LBL 03
38 RCL M 39 SIGN 40 X<> N 41 RCL O 42 CLA 43 END |
( 75 bytes )
| STACK | INPUTS | OUTPUTS |
| Z | bbb.eee | / |
| Y | y | y |
| X | x | p(x,y) |
| L | / | x |
with bbb > 000
Example: p(x,y) = 1 + 2 x + 3 y + 3 x2 + 2 x.y - 4 y2 + 2 x3 + 4 x2y - x.y2 + 3 y3 >>>> Evaluate p(3,4)
-Suppose we choose the first register = R11
-Store 1 2
3 3 2
-4 2
4 -1 3
into R11 R12 R13 R14
R15 R16 R17 R18 R19 R20 respectively
( control number = 11.020 )
11.020 ENTER^
4
ENTER^
3
XEQ "PXY" >>>> p(3;4) = 348
( in 4 seconds )
b) Degree of a Polynomial
-This short routine calculates the degree of a polynomial from its control
number and it produces DATA ERROR ( line 20 ) if this control number is
invalid.
-It is used to compute the product of 2 polynomials in paragraph 1-c)
Data Registers: /
Flags: /
Subroutines: /
| 01 LBL "DEGXY"
02 INT 03 LASTX 04 FRC 05 E3 06 * 07 1 08 + 09 X<>Y 10 - 11 8 12 * 13 1 14 + 15 SQRT 16 3 17 - 18 2 19 / 20 FACT 21 X<> L 22 END |
( 34 bytes )
Examples:
11.020 XEQ "DEGXY" yields
3
1.028
R/S
yields 6
1.029
R/S
yields DATA ERROR 1.029 is not a valid control-number
c) Product of 2 Polynomials
-Storing the coefficients as explained above allows to write a relatively
fast program to multiplicate 2 polynomials of 2 variables:
Data Registers: R00 to R07: temp ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PROXY" )
• Rbb ...... • Ree = coefficients of the
1st polynomial bb > 07
• Rbb' ..... • Ree' = coefficients of the 2nd polynomial
bb' > 07
Flags: /
Subroutine: "DEGXY"
-If you don't have an HP-41CX, replace lines 27-28-29 by
SIGN
CLX
STO 00
LBL 00
STO IND L
ISG L
GTO 00
| 01 LBL "PROXY"
02 STO 03 03 RDN 04 STO 02 05 X<>Y 06 STO 01 07 XEQ "DEGXY" 08 1 09 + 10 X<>Y 11 XEQ "DEGXY" 12 ST+ Z 13 + 14 1 15 + |
16 *
17 2 18 / 19 RCL 03 20 + 21 DSE X 22 E3 23 / 24 RCL 03 25 + 26 STO 03 27 CLRGX 28 CLX 29 STO 00 30 LBL 01 |
31 RCL 03
32 RCL 00 33 + 34 STO 04 35 LASTX 36 8 37 * 38 1 39 ST+ 00 40 + 41 SQRT 42 1 43 - 44 2 45 / |
46 INT
47 STO 05 48 CLX 49 STO 06 50 RCL 02 51 STO 07 52 LBL 02 53 RCL IND 01 54 RCL IND 07 55 * 56 ST+ IND 04 57 ISG 06 58 GTO 03 59 RCL 05 60 ST+ 04 |
61 RCL 06
62 INT 63 E3 64 / 65 STO 06 66 LBL 03 67 ISG 04 68 CLX 69 ISG 07 70 GTO 02 71 ISG 01 72 GTO 01 73 RCL 03 74 END |
( 113 bytes )
| STACK | INPUTS | OUTPUTS |
| Z | bbb.eee | / |
| Y | bbb.eee' | / |
| X | bbb'' | bbb.eee'' |
with bbb and bbb' > 007
Example:
p(x,y) = 1 + 2 x + 3 y + 3 x2 + 2 x.y
- 4 y2 + 2 x3 + 4 x2y - x.y2
+ 3 y3
q(x,y) = 2 + 3 x + 4 y + 5 x2 + 6 x.y
+ 7 y2
-If we choose Rbb = R11 for p(x,y) and Rbb' = R21 for q(x,y)
-Store 1 2
3 3 2
-4 2
4 -1 3
into R11 R12 R13 R14
R15 R16 R17 R18 R19 R20 respectively
( control number = 11.020 )
and 2
3 4
5 6
7
into R21 R22 R23
R24 R25 R26 respectively
( control number = 21.026 )
-We can take for instance Rbb" = R31
11.020 ENTER^
21.026 ENTER^
31
XEQ "PROXY" >>>> bbb.eee'' = 31.051
( in 47 seconds ) and we have the coefficients of the product
in registers R31 thru R51, namely:
p(x,y) q(x,y) = 2 + 7 x + 10 y + 17 x2 + 27 x.y
+ 11 y2 + 23 x3 + 53 x2y + 26 x.y2
+ 11 y3 + 21 x4 + 48 x3y + 26 x2y2
- 5 x.y3 -16 y4
+ 10 x5 + 32 x4y + 33 x3y2
+ 37 x2y3 + 11 x.y4 + 21 y5
2°) N Variables
-In the following programs, each monomial µ xa
yb zc td ue is stored
into 2 consecutive registers
-The coefficient µ is stored in one data register and the
exponents are coded in the next register under the form aa.bbccddee
-For instance, to store 12 x2 y3 z4 t15 u6 in registers R14-R15 ,
12
STO 14
2.03041506 STO 15
-Thus, we assume that all the ( partial ) degrees are smaller than 99
and that the number of variables does not exceed 5.
-However, if we are sure that all the degrees will be smaller than
10, we could deal with polynomials of up to 10 variables ( x2
y3 z4 t5 u6 v7 is
coded 2.34567 )
-On the other hand, if you want to handle degrees up to 999, the maximum
number of variables will be reduced to 3 and each exponent will be coded
with 3 digits.
( x21 y34 z471 will be coded
21.034471 )
a) Evaluation of Polynomial
Functions
Data Registers: R00: temp ( Registers R01 thru Rnn & Rbb thru Ree are to be initialized before executing "PEVAL" )
• R01 = x • R02 = y ......
• Rbb ...... • Ree = coefficients of the polynomial
Flags: /
Subroutines: /
-If you don't want to use the synthetic registers M N O , replace
them by R06 R07 R08
-In this case, delete line 32 and replace line 01 by 0
STO 06 RDN
-Replace line 21 by 10 if you code the exponents with one
digit - provided they are smaller than 10
or by E3 if you code the exponents with
3 digits - provided they are smaller than 1000
| 01 LBL "PEVAL"
02 CLA 03 STO 00 04 STO O 05 LBL 01 06 CLX 07 STO N |
08 RCL IND 00
09 ISG 00 10 RCL IND 00 11 LBL 02 12 ISG N 13 CLX 14 FRC |
15 RCL IND N
16 LASTX 17 INT 18 Y^X 19 ST* Z 20 CLX 21 E2 |
22 *
23 X#0? 24 GTO 02 25 X<>Y 26 ST+ M 27 ISG 00 28 GTO 01 |
29 RCL O
30 SIGN 31 X<> M 32 CLA 33 END |
( 58 bytes )
| STACK | INPUTS | OUTPUTS |
| X | bbb.eee | p(x,y,z,...) |
| L | / | bbb.eee |
with bbb > nnn ( n = number of variables )
Example: p(x,y,z) = 3 x y4 z7 + 4 x2 y2 z3 - 7 y5 z + 2 Evaluate p(2,3,4)
-If you choose to store the first coefficient into R11
3 STO 11 1.0407
STO 12
4 STO 13 2.0203
STO 14
-7 STO 15 0.0501 STO 16
2 STO 17
0 STO 18
( control number = 11.018 )
2 STO 01
3 STO 02
4 STO 03
11.018 XEQ
"PEVAL" >>>> p(2,3,4) = 7965038 ( in 8 seconds
)
Note:
-There will be a DATA ERROR line 18 if a variable and a corresponding
exponent equal zero because the HP-41 refuses to compute 00
-So, you can add after line 17 X#0? GTO 02
SIGN STO Y LBL 02 and you'll get for
instance the correct result p(0,3,4) = -6802
-You can also replace line 18 by a M-Code routine Y^X that does produce 00 = 1 for example:
098 "X"
01E "^"
019 "Y"
0F8 READ3(X)
10E A=C ALL
0B8 READ2(Y)
355 ?NCXQ
this subroutine checks A and C for alphadata
050 14D5
executes A<>C and SETDEC
128 WRIT4(L) saves X in L-register
2EE ?C#0? ALL
02F JC+ 05
04E C=0 ALL C
35C PT= 12
=
050 LD@PT-1 1
023 JNC+04
044 CLRF 4
C
045 ?NCXQ
=
06C 1B11
A^C
0E8 WRIT3(X)
078 READ1(Z)
0A8 WRIT2(Y)
046 C=0 S&X
270 RAMSLCT
038 READDATA
068 WRIT1(Z)
3E0 RTN
XEQ "Y^X" will work like the built-in function Y^X but it yields 00 = 1
-However, this M-Code routine does not check for overflows and underflows.
b) Product of 2 Polynomials
Data Registers: R00 to R05: temp ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PMUL" )
• Rbb ...... • Ree = coefficients of the 1st polynomial
bb > 05
• Rbb' ..... • Ree' = coefficients of the 2nd polynomial
bb' > 05
Flag: F07 is cleared
Subroutine: "DEL0" ( cf next paragraph
)
| 01 LBL "PMUL"
02 STO 03 03 DSE X 04 E-3 05 STO 05 06 ST+ 05 07 * 08 ST+ 03 09 RDN 10 STO 02 11 STO 04 12 X<>Y 13 STO 01 14 SIGN |
15 2 E-5
16 ST+ 01 17 + 18 STO 00 19 SF 07 20 LBL 01 21 RCL IND 01 22 RCL IND 04 23 * 24 RCL 00 25 RCL 01 26 + 27 ISG 04 28 RDN |
29 RCL IND T
30 RCL IND 04 31 + 32 RCL 00 33 RCL 03 34 + 35 X<>Y 36 FS?C 07 37 GTO 03 38 LBL 02 39 RCL IND Y 40 X=Y? 41 GTO 04 42 RDN |
43 ISG Y
44 GTO 02 45 LBL 03 46 STO IND Y 47 SIGN 48 - 49 X<>Y 50 STO IND Y 51 RCL 05 52 ST+ 03 53 GTO 05 54 LBL 04 55 RDN 56 SIGN |
57 -
58 X<>Y 59 ST+ IND Y 60 LBL 05 61 ISG 04 62 GTO 01 63 RCL 02 64 STO 04 65 ISG 01 66 GTO 01 67 RCL 03 68 XEQ "DEL0" 69 END |
( 113 bytes )
| STACK | INPUTS | OUTPUTS |
| Z | bbb.eee | / |
| Y | bbb.eee' | / |
| X | bbb'' | bbb.eee'' |
with bbb and bbb' > 005
Example: p(x,y,z) = 2 x + 3 x.y2z + 4 x3y.z4 q(x,y,z) = 4 - 6 y2z + x2y.z4
-With bb = 07 & bb' = 13
2 STO 07
3 STO 09
4 STO 11 control
number
4 STO 13 -6
STO 15 1
STO 17 control number
1 STO 08 1.0201
STO 10 3.0104 STO 12
= 7.012
0 STO 14 0.0201 STO 16
2.0104 STO 18
= 13.018
-If you choose bb'' = 21
7.012 ENTER^
13.018 ENTER^
21 XEQ
"PMUL" >>>> 21.030 ( execution time
= 15 seconds ) and we have:
R21 = 8
R23 = 4
R25 = 18
R27 = -18 R29
= -21
R22 = 1
R24 = 5.0208 R26 = 3.0104
R28 = 1.0402 R30 = 3.0305
-So, p(x,y,z) q(x,y,z) = 8 x + 4 x5y2z8 + 18 x3y.z4 - 18 x.y4z2 - 21 x3y3z5
Notes:
-Each new monomial is compared to the other ones to collect like terms
as soon as possible ( lines 38 to 44 & 54 to 59 )
-It reduces the number of registers during the calculations.
-A shorter program may be written if the HP-41 simplifies the polynomial
at the end
but sometimes, ( too? ) many registers are used, especially
if there are many like-terms.
-"SIM" is listed in paragraph 2-d)
| 01 LBL "PMUL2"
02 STO 00 03 STO 03 04 RDN 05 STO 02 06 STO 04 07 X<>Y 08 STO 01 |
09 LBL 01
10 RCL IND 01 11 RCL IND 04 12 * 13 STO IND 03 14 1 15 ST+ 03 16 ST+ 04 |
17 RCL 01
18 + 19 RCL IND X 20 RCL IND 04 21 + 22 STO IND 03 23 ISG 03 24 CLX |
25 ISG 04
26 GTO 01 27 RCL 02 28 STO 04 29 ISG 01 30 ISG 01 31 GTO 01 32 RCL 03 |
33 DSE X
34 E3 35 / 36 RCL 00 37 + 38 XEQ "SIM" 39 END |
( 69 bytes )
| STACK | INPUTS | OUTPUTS |
| Z | bbb.eee | / |
| Y | bbb.eee' | / |
| X | bbb'' | bbb.eee'' |
with bbb and bbb' > 004
-The same example is solved in 21 seconds ( instead of 15 )
-Note that the monomials are stored in a different order.
-But if you are sure that there is no like-term in the product, delete
line 38 and "PMUL2" will be much faster than "PMUL"
c) Deleting the Zeros
-This subroutine is called by "PMUL" and "SIM" to eliminate the monomials
that equal zero.
Data Registers: R00 to R02: temp ( Registers Rbb thru Ree are to be initialized before executing "DEL0" )
• Rbb ...... • Ree = coefficients of the polynomial
bb > 03
Flags: /
Subroutines: /
| 01 LBL "DEL0"
02 STO 00 03 STO 01 04 FRC 05 E3 06 * 07 STO 02 08 2 E-5 |
09 ST+ 00
10 LBL 01 11 ISG 00 12 FS? 30 13 GTO 02 14 RCL IND 00 15 X#0? 16 GTO 01 |
17 RCL 02
18 DSE X 19 RCL IND X 20 STO IND 00 21 RCL 00 22 1 23 + 24 RCL IND 02 |
25 STO IND Y
26 2 27 ST- 00 28 ST- 02 29 E3 30 / 31 ST- 00 32 ST- 01 |
33 GTO 01
34 LBL 02 35 RCL IND 01 36 X=0? 37 STO IND 02 38 RCL 01 39 END |
( 71 bytes )
| STACK | INPUTS | OUTPUTS |
| X | bbb.eee | bbb.eee' |
with bbb > 003
- bbb.eee = bbb.eee' if all the coefficients
are different from zero
-If all the coefficients equal zero eee' = bbb+1 &
Rbb = Ree' = 0
d) Simplifications
-This subroutine collects like terms to simplify a polynomial expression.
Data Registers: R00 to R04: temp ( Registers Rbb thru Ree are to be initialized before executing "SIM" )
• Rbb ...... • Ree = coefficients of the polynomial
bb > 04
Flags: /
Subroutine: "DEL0" ( cf § 2-c)
)
-Line 51 may be replaced by GTO "DEL0"
| 01 LBL "SIM"
02 STO 01 03 1.00002 04 + 05 STO 02 06 RCL 01 07 FRC 08 E3 09 * 10 STO 00 11 LBL 01 |
12 RCL 02
13 STO 03 14 RCL IND X 15 LBL 02 16 ISG 03 17 FS? 30 18 GTO 03 19 RCL IND 03 20 X<>Y 21 X#Y? 22 GTO 02 |
23 STO 04
24 RCL IND 00 25 STO IND 03 26 RCL 02 27 1 28 ST- 00 29 ST- 03 30 - 31 RCL IND 00 32 X<> IND 03 33 ST+ IND Y |
34 2 E-3
35 ST- 02 36 ST- 03 37 DSE 00 38 SIGN 39 ST- 03 40 RCL 04 41 GTO 02 42 LBL 03 43 ISG 02 44 GTO 01 |
45 RCL 00
46 E3 47 / 48 RCL 01 49 INT 50 + 51 XEQ "DEL0" 52 END |
( 96 bytes )
| STACK | INPUTS | OUTPUTS |
| X | bbb.eee | bbb.eee' |
with bbb > 004
Example: p(x,y) = 2 x + 3 x2y + 4 x.y + 5 x2y2 + 8 x.y + 7 x2y
-If we store this polynomial into R05 R06 .....
2 STO 05
3 STO 07 4
STO 09 5 STO 11
8 STO 13 7
STO 15
1 STO 06 2.01
STO 08 1.01 STO 10 2.02
STO 12 1.01 STO 14
2.01 STO 16
5.016 XEQ "SIM" >>>> 5.012 ( in 8 seconds ) and we find
R05 = 2 R07 = 10
R09 = 12 R11 = 5
R06 = 1 R08 = 2.01
R10 = 1.01 R12 = 2.02
-So p(x,y) = 2 x + 10 x2y + 12 x.y + 5 x2y2
Note:
-This program may be used to add 2 polynomials:
-Store the 2 polynomials into consecutive registers so that their control
numbers (bbb.eee)1 and (bbb.eee)2 verify
bbb2 = 1+eee1
-Then, key in bbb1.eee2
XEQ "SIM"
e) Partial Derivatives
-This program calculates ¶ a+b+c+...p
/ ¶xa¶yb¶zc....
[ d a+b+c+...p / dxadybdzc....
if the symbol ¶ doesn't appear correctly
]
Data Registers: R00 to R04: temp ( Registers Rbb thru Ree are to be initialized before executing "PDRV" )
• Rbb ...... • Ree = coefficients of the polynomial
bb > 04
Flags: /
Subroutine: "SIM" ( cf § 2-d)
)
-Replace line 31 by 10 if you code the exponents with one
digit - provided they are smaller than 10
or by E3 if you code the exponents with
3 digits - provided they are smaller than 1000
| 01 LBL "PDRV"
02 STO 00 03 X<>Y 04 STO 01 05 STO 02 06 ISG 02 07 LBL 01 08 RCL IND 02 09 RCL 00 |
10 ST- IND 02
11 DSE 02 12 LBL 02 13 RCL Y 14 INT 15 ST- Z 16 RCL Y 17 INT 18 ST- Z |
19 X=0?
20 GTO 04 21 X<>Y 22 LBL 03 23 ST* IND 02 24 DSE X 25 "" 26 DSE Y 27 GTO 03 |
28 LBL 04
29 RDN 30 CLX 31 E2 32 ST* Z 33 * 34 X#0? 35 GTO 02 36 2 |
37 ST+ 02
38 ISG 02 39 GTO 01 40 RCL 01 41 XEQ "SIM" 42 END |
( 74 bytes )
| STACK | INPUTS | OUTPUTS |
| Y | bbb.eee | / |
| X | aa.bbccddee | bbb.eee' |
with bbb > 004
Example: p(x,y,z) = 2 x3 y4 z5 + x2 y z3 + x y z3 + 3 x12 y7 z10 + 4 Calculate ¶6p / ¶x2¶y ¶z3 ( aa.bbcc = 2.0103 )
-If p(x,y,z) is stored in R05 R06 ....
2
STO 05 1
STO 07 1
STO 09 3
STO 11 4 STO 13
3.0405 STO 06 2.0103
STO 08 1.0103 STO 10
12.0710 STO 12 0 STO 14
control number = 5.014
5.014 ENTER^
2.0103 XEQ "PDRV" >>>>
5.010 ( in 24 seconds ) and we have
R05 = 2880 R07 = 12
R09 = 1995840
R06 = 1.0302 R08 = 0
R10 = 10.0607
-Thus, ¶6p / ¶x2¶y ¶z3 = 2880 x y3 z2 + 12 + 1995840 x10 y6 z7
Note:
-This routine replaces the polynomial by its derivative.
-So you'll have to store it in a second block of registers if you want
to save it.