The Museum of HP CalculatorsHP Articles Forum
[Return to the Index ]
[ Previous | Next ]
A "Cadillac " Quadratic Solver for the hp 33s and HP35s
Posted by Palmer O. Hanson, Jr. on 29 Apr 2008, 11:40 a.m.
NOTE: The program for the hp33s was previously published as Article 593 on 10 March 2006. The program sequence had errors at A0008, B0022 and E0009. I was unable to correct those errors because I did not have a valid password. This article repeats the hp 33s material with the errors corrected and adds a similar program for the HP35s.
Thhe Polynomial Root Finder program in pages 15-20 through 15-32 of the hp 33s user's guide offers a quadratic solution as an option. That routine is a "garden variety" solution which is adequate for most quadratic problems; however, users may occasionally encounter quadratic problems which cannot be correctly solved with that routine. An example is a problem for which the discriminant is erroneously found to be zero because the two parts of the discriminant have more significant digits than the number of significant digits carried by the machine on which the quadratic solver is installed.
A quadratic solver which appears on pages 208-211 of the HP-15C Advanced Functions Handbook uses a so-called "tricky property of the SUM + and SUM - commands" in a way that is equivalent to carrying twenty significant digits in the calculation of the discrimiinant. As I was working my way through the hp 33s user's guide I found the following statement on page 1-19: "... During some complicated internal calculations, the calculator uses 15-digit precision for immediate results. ..." That sounded a lot like the "tricky properties". I did some tests that convinced me that the "tricky properties" were available in the hp 33s. I then decided to try to translate the HP-15C program for use with the hp 33s. Steps 011 through 025 of the HP-15C program are
011 CLEAR SUM
012 RCL 8
013 STO 7
014 RCL / I
015 RND
016 RCL I
017 SUM -
018 RCL 9
019 x <> 7
020 x <> y
021 RCL 8
022 SUM -
023 R Down
024 SUM -
025 RCL 7
and so on where there is an additional SUM - at step 052 and there are no SUM + commands in the loop. That means that the HP-15C is able to accept inputs to the statistics registers which would cause the number of entries (n) to become negative. The statistics routines in the HP-41, the HP-32S and several other HP machines can also do that. The statistics routine in the hp 33s cannot. After a little thought it occurred to me that hp 33s routine may not care how the value for the number of entries (n) had been accumulated in register 28. To prove that I did a CL SUM and followed it by three SUM + commands where 2 was in the y register and 3 was in the x register. I recalled the means and saw the values 3 and 2 as expected. Then I stored 28 (the number of the statistics register which holds the value of n) in I, entered a 3 in the display, and did a STO+ (i). I recalled the means and saw the values 1.5 and 1.0 . To make a long story short I was able to successfully translate the HP-15C program for use on the hp 33s by storing a four for n in statistics register 28 after the CL SUM command and before the SUM - commands. The program which follows implements the HP-15C program on the hp 33s. Actually, it is not a direct translation from the HP-15C program but rather is a translation from my translation of the HP-15C program for use on the HP-41 as described in my submission "Update: 13 Digit Precision on the HP-15C" of 23 July 2003 which appears in Archive 13 of the Forum.
The HP-33s and HP-35s programs below accept the coefficients A, B, and C of the quadratic Ax^2 + Bx + C = 0. For complex roots Flag zero is set and the imaginary part is in the upper part of the display. The entered value for B is divided by -2 to convert it for use by the solution routine which is for the modified quadratic Ax^2 - 2Bx + C = 0. This is similar to the quadratic solution method on pages 110-111 of the HP-27 Owner's Handbook. The value of A does not have to be a value of one as with the quadratic routine from the hp 33s Users Guide. Some test cases follow the program listing . You will find that the quadratic program from the hp 33s user's guide fails to get the correct answer for several of the cases. The translated program gets the correct answer in all eight cases.
HP-33S Program Listing
H0001 LBL H
H0002 INPUT A
H0003 STO D
H0004 INPUT B
H0005 -2
H0006 /
H0007 STO F
H0008 STO H
H0009 INPUT C
H0010 STO G
H0011 STO I
H0012 CF 0
H0013 SCI 2
A0001 LBL A
A0002 CL SUM Clear statistics registers
A0003 28
A0004 STO i
A0005 4
A0006 STO(i)
A0007 33
A0008 STO i
A0009 RCL H
A0010 STO(i)
A0011 RCL D
A0012 /
A0013 RND
A0014 RCL D
A0015 SUM -
A0016 RCL I
A0017 SUM xy Recall statistics product from register 33
A0018 x<>y
A0019 STO(i)
A0020 R DOWN Roll down
A0021 x<>y
A0022 RCL H
A0023 SUM -
A0024 R DOWN
A0025 SUM -
A0026 SUM xy
A0027 ABS
A0028 RCL I
A0029 ABS
A0030 x<=y? Is x less than or equal to y
A0031 GTO B
A0032 ENTER
A0033 R UP Roll up
A0034 STO H
A0035 SUM xy
A0036 STO I
A0037 ABS
A0038 1 E 20
A0039 X
A0040 RCL G
A0041 ABS
A0042 x<=y?
A0043 GTO A
B0001 LBL B
B0002 FIX 9
B0003 RCL H
B0004 x^2
B0005 STO(i)
B0006 RCL D
B0007 RCL I
B0008 SUM -
B0009 RCL F
B0010 SUM xy
B0011 x<0?
B0012 GTO E
B0013 SQR x Square root of x
B0014 RCL F
B0015 SGN
B0016 X
B0017 +
B0018 STO J
B0019 x<>y
B0020 /
B0021 RCL J
B0022 x<>0? x is not equal to zero
B0023 GTO F
B0024 CLx
B0025 STOP
F0001 LBL F
F0002 RCL G
F0003 x<>y
F0004 /
F0005 STOP
E0001 LBL E
E0002 SF 0
E0003 +/-
E0004 SQR x
E0005 RCL D
E0006 /
E0007 x<>y
E0008 R UP
E0009 /
E0010 STOP
HP35S Program Listing
The hp 33s program listing above requires a number of changes for use on the HP35s; for example, to account for the different assignment of memory locations for the statistics registers A0003 must be changed from +28 to -27 and A0007 must be changed from +33 to -32. Changes are also required to account for the difference between i and (i) with the hp 33s and I and (I) with the HP35s. It would also be necessary to assign a different variable name wherever I is used in the HP 33s program. Finally, the hp 33s version uses five labels because direct address transfers are not available. The result of all of this is that I rewrote the hp 33s program for the HP35s using direct addressing to save labels, and in the process made the other necessary changes.
Q001 LBL Q
Q002 SF 10
Q003 EQN Ax2 + Bx + C
Q004 CF 10
Q005 INPUT A
Q006 STO D
Q007 INPUT B
Q008 -2
Q009 /
Q010 STO F
Q011 STO H
Q012 INPUT C
Q013 STO G
Q014 STO K
Q015 CF 0
Q016 SCI 2
Q017 CL SUM Clear statistics registers
Q018 -27
Q019 STO I
Q020 4
Q021 STO (I)
Q022 -32
Q023 STO I
Q024 RCL H
Q025 STO (I)
Q026 RCL D
Q027 /
Q028 RND
Q029 RCL D
Q030 SUM -
Q031 RCL K
Q032 SUM xy Recall statistics product from register 33
Q033 x<>y
Q034 STO (I)
Q035 R DOWN Roll down
Q036 x<>y
Q037 RCL H
Q038 SUM -
Q039 R DOWN
Q040 SUM -
Q041 SUM xy
Q042 ABS
Q043 RCL K
Q044 ABS
Q045 x<=y? Is x less than or equal to y
Q046 GTO Q059
Q047 ENTER
Q048 R UP Roll up
Q049 STO H
Q050 SUM xy
Q051 STO K
Q052 ABS
Q053 1 E 20
Q054 X
Q055 RCL G
Q056 ABS
Q057 x<=y?
Q058 GTO Q017
Q059 FIX 9
Q060 RCL H
Q061 x^2
Q062 STO (I)
Q063 RCL D
Q064 RCL K
Q065 SUM -
Q066 RCL F
Q067 SUM xy
Q068 x<0?
Q069 GTO Q087
Q070 SQR x Square root of x
Q071 RCL F
Q072 SGN
Q073 X
Q074 +
Q075 STO L
Q076 x<>y
Q077 /
Q078 RCL L
Q079 x<>0? x is not equal to zero
Q080 GTO Q083
Q081 CLx
Q082 STOP
Q083 RCL G
Q084 x<>y
Q085 /
Q086 STOP
Q087 SF 0
Q088 +/-
Q089 SQR x
Q090 RCL D
Q091 /
Q092 x<>y
Q093 R UP
Q094 /
Q095 STOP
TEST CASES
The HP-15C Advanced Functions Handbook describes the program as "like Grandmother's expensive chinaware, reserved for special occasions." I suspect that may be because it uses more memory and much more time to run. For the test cases below the hp 33s never uses more than a second or two to reach a solution. The HP35s runs somewhat longer.
Cases 1 through 4 are equivalent to cases 3 through 6 on page 207 of the HP-15C Advanced Functions Handbook . Cases 5 and 6 and are equivalent to the problems defined on page 211 of HP-15C Advanced Functions Handbook. Cases 7 and 8 are two additional problems that I devised which require a 15 digit capability if the discriminant is to be calculated correctly.
Case 1:
a = 1E-13 b = -2 c = 1
R1 = 2E13 R2 = 0.5
Case 2:
a = 654,323 b = -1,308,644 c = 654,321
R1 = 1 R2 = 0.999996943...
Case 3:
a = 11,713 b = -1,470,492 c = 46,152,709
Re = 62.77179203... ; Im = 8.5375E-05
Case 4:
a = 80,841 b = -1,975,288 c = 12,066,163
Re = 12.21711755... ; Im = 0.001374514...
Case 5:
a = 4,877,361,379 b = -9,754,525,226 c = 4,877,163,849
Re = 0.999979750 ; Im = 2.8995463E-10
Case 6:
a = 1 b = -222,223 c = 12,193,329,370
R1 = 123,458 ; R2 = 98765
Case 7:
a = 11,111,119 b = -22,222,222 c = 11,111,103
R1 = 1 ; R2 = 0.999998560001
Case 8::
a = 8,441,600 b = -22,222,222 c = 14,624,809
Re = 1.31623282316 Im = ± 1.05290400129E-6
Go back to the main exhibit hall