The Museum of HP Calculators

HP Articles Forum

[Return to the Index ]
[ Previous | Next ]

Eighth Order Linear Equation Solver for the hp 33s

Posted by Palmer O. Hanson, Jr. on 6 Oct 2006, 12:24 a.m.

I started with a TI-59 linear equation solution program by Henrik Ohlsson which appeared in the 81-2 issue of the Swedish newsletter Programbiten. Henrik's program expanded the capability of the TI-59 from order 8 to order 16.

I translated Henrik's program for the Durabrand 828 as a means of demonstrating the capability of that little low cost machine. I was able to fit a sixth order linear equation solution into the Durabrand. The translation was essentially a straightforward, brute force modification from AOS to EOS.

I translated the Durabrand program into a sixth order linear equation solution for the hp 33s. The resulting hp 33s program was essentially a brute force translation of Henrik's program from the TI-59 for the Durabrand and from the Durabrand for the hp 33s. The resulting sixth order linear equation solution was documented as Article 676. The hp 33s program wasn't pretty but it was functional. Only the 26 data registers were used. Register arithmetic was used throughout (no use of the stack in RPN and no use of parentheses or equal signs in ALG) so the the program as written would run in either RPN or ALG mode. That wasn't necessary but it was something that I had wanted to demonstrate.

Next, I found a way to use the statistics registers to expand the capability of the hp 33s program to solve a seventh order set of linear equations. Only limited changes were required to make the expansion. The changes were posted in a September 17 message "Expansion to seventh order" in the Forum.

Finally, I expanded the solution to solve an eighth order set of linear equations. To do that Henrik's method required 33 data registers. The combination of the data registers and the statistics registers of the hp 33s yields only 32 registers. I found a modification to Henrik's method which eliminated the need for one register. The sixth and seventh order solutions had used data registers Y and Z as scratch pads in order to implement the register arithmetic type solution. I had to free up those registers to yield the 32 registers needed for my modification to Henrik's method. To do that I changed substantial parts of the program to use the RPN stack. The resulting program follows.

I won't be doing an expansion to ninth order. Henrik's method would require five more data registers to do that.

Program Listing:

N0001   LBL N
N0002   CLVARS
N0003   CLSum
N0004   INPUT N
N0005   STO A
N0006   0
N0007   STO N
N0008   8
N0009   STO C
O0001   LBL O
O0002   0
O0003   STO G
O0004   RCL A
O0005   STO B
O0006   1
O0007   STO+ B
O0008   8
O0009   STO F
O0010   RCL D
O0011   IP
O0012   STO E
O0013   STO D
O0014   1
O0015   STO+ D
O0016   RCL E
O0017   x=0?
O0018   GTO C
A0001   LBL A
A0002   XEQ S
A0003   +/-
B0001   LBL B
B0002   ENTER
B0003   ENTER
B0004   RCL F
B0005   XEQ T
B0006   R Down
B0007   RCL (i)
B0008   x
B0009   RCL C   
B0010   XEQ T
B0011   R Down
B0012   STO+ (i)
B0013   R Down
B0014   R Down
B0015   1
B0016   STO+ C
B0017   STO+ F
B0018   R Down
B0019   DSE B
B0020   GTO B
B0021   RCL A
B0022   STO B
B0023   STO- C
B0024   1
B0025   STO+ B
B0026   STO- C
B0027   DSE E
B0028   GTO A
C0001   LBL C
C0002   XEQ S
C0003   STO E
C0004   RCL C
C0005   XEQ T
C0006   RCL (i)
C0007   STO+ E
C0008   1
C0009   STO- B
D0001   LBL D
D0002   XEQ S
D0003   ENTER
D0004   1
D0005   STO+ C
D0006   R Down
D0007   RCL C
D0008   XEQ T
D0009   R Down
D0010   RCL (i)
D0011   +
D0012   0
D0013   STO (i)
D0014   1
D0015   STO- C
D0016   R Down
D0017   R Down
D0018   RCL C
D0019   XEQ T
D0020   R Down
D0021   STO (i)
D0022   RCL E
D0023   STO/ (i)
D0024   1
D0025   STO+ C
D0026   DSE B
D0027   GTO D
D0028   RCL D
D0029   IP
D0030   STO E
D0031   DSE E
D0032   GTO I
D0033   GTO H
I0001   LBL I
I0002   7
I0003   STO C
I0004   RCL A
I0005   STO B
E0001   LBL E
E0002   1
E0003   STO+ C
E0004   STO+ G
E0005   RCL C
E0006   XEQ T
E0007   0
E0008   x<> (i)
F0001   LBL F
F0002   ENTER
F0003   ENTER
F0004   1
F0005   STO+ C
F0006   R Down
F0007   RCL F
F0008   XEQ T
F0009   R Down
F0010   RCL (i)
F0011   +/-
F0012   x
F0013   RCL C
F0014   XEQ T
F0015   R Down
F0016   RCL (i)
F0017   +
F0018   0
F0019   STO (i)
F0020   R Down
F0021   RCL G
F0022   STO- C
F0023   R Down
F0024   RCL C
F0025   XEQ T
F0026   R Down
F0027   STO (i)
F0028   RCL G
F0029   STO+ C
F0030   R Down
F0031   1
F0032   STO+ F
F0033   R Down
F0034   R Down
F0035   DSE B
F0036   GTO F
F0037   RCL A
F0038   STO- F
F0039   STO B
F0040   DSE E
F0041   GTO E
F0042   RCL G
F0043   STO- C
F0044   1
F0045   STO+ C
G0001   LBL G
G0002   RCL F
G0003   XEQ T
G0004   RCL (i)
G0005   RCL C
G0006   XEQ T
G0007   R Down
G0008   STO (i)
G0009   RCL F
G0010   XEQ T
G0011   0
G0012   STO (i)
G0013   1
G0014   STO+ C
G0015   STO+ F
G0016   DSE B
G0017   GTO G
H0001   LBL H
H0002   DSE A
H0003   GTO O
H0004   7
H0005   STO i
H0006   RCL D
H0007   IP
H0008   STO+ i
H0009   1000
H0010   STO/ i
H0011   8
H0012   STO+ i
J0001   LBL J
J0002   RCL (i)
J0003   STOP
J0004   ISG i
J0005   GTO J
J0006   0
J0007   1/x
J0008   STOP
S0001   LBL S
S0002   0.01
S0003   STO+ D
S0004   RCL D
S0005   FIX 2
S0006   STOP
S0007   FIX 9
S0008   RTN
T0001   LBL T
T0002   STO i
T0003   27
T0004   -
T0005   x<0?
T0006   RTN
T0007   R Down
T0008   1
T0009   STO+ i
T0019   RTN
To use the program press XEQ N and see the prompt N? in the display. Enter the order of the solution desired (n = 2 to 6) and press R/S. 1.01 will appear in the display indicating that the machine is ready to accept the A(1,1) element of the matrix Enter the A(1,1) element and press R/S. 1.02 will appear in the display indicating that the machine is ready to accept the A(1,2) element of the matrix. Procced as before to enter the A(1,2) element and see the prompt for the A(1.3) element, and so on. When the column element of the prompt is n+1 enter the corresponding element of the vector. When the last element of the vector has been entered the machine will stop with 8.000000000 in the Y register display and the first element of the solution, X(1), in the X register display. Press R/S to see X(2), and so on. When all of elements of the solution vector have been displayed the message DIVIDE BY 0 will be displayed. The elements of the solution can be recalled starting with register H followed by subsequent registers as needed by the order of the system.. .

If you run the problem with the matrix as a 8x8 sub-Hilbert and the vector all ones you will get the following results:

   Exact                  hp 33s                 Relative Error

-72 -79.85564194 0.10674 2520 2752.76793665 0.09237 -27720 -29976.5053369 0.08140 138600 148685.022829 0.07276 -360360 -384064.400797 0.06578 504504 534783.269832 0.06002 -360360 -380245.814081 0.05518 102960 108218.029350 0.05107

for a mean relative error of 7.316E-02. Hendrik's program running on the TI-59 yields a mean relative error of 3.331E-02. The ML-02 program running on the TI-59 yields a mean relative error of 2.212E-03. The matrix solution in the HP-28 yield a mean relative error of 1.801E-03. .

Some notes on the program:

1. With Henrik's program on the TI-59 the calculations after the entry of an element required a significant amount of time so that the user was always waiting on the machine. For example, in a fifth order solution when the vector element is entered in response to the prompt 3.06 the TI-59 will run for nine seconds in fast mode or fifteen seconds in normal mode before the next prompt (4.01) appears. The much higher calculating speed of the hp 33s yields much shorter wait times.

2. Henrik used an INV Dsz command in his TI-59 program. I had to synthesize the equivalent in the hp 33s which yielded the curious code at steps D0021 through D0033..

3. The T routine provides the methodology for easily adding the statistics registers immediately after the data registers, essenially skipping the i register.

4. The STO/ (i) command at step D0023 will cause a divide by zero error for some systems. Consider the following problem:

   2   -1   -3     1        1

-2 1 1 3 1

2 -1 0 -1 1

-3 2 -1 0 1

where the solution 5.0, 8.25, 0.5, 0.75 can be readily obtained with the TI-59 ML-02 program, with the H-41 Math Pac and with others. If you attempt to solve the problem with this hp 33s program you will get a "Divide by 0" error after the A(2,3) element is entered. If you start the program over again and enter the first row and then the third row you will again get the"Divide by 0" error after the A(2,3) element is entered. However, if you start the program one more time and enter the first row, followed by the fourth row, the second row and the third row in that order you will see the correct solution. The point of all this is that if you get the "Divide by 0" error after the entry of an element you must try another arrangement of the rows.


[ Return to the Message Index ]

Go back to the main exhibit hall