(35S) - Quadratic root finder -- high performance

[Douglas De Boer]

In the attached PDF file I offer an HP 35s program to find the roots of a quadratic equation. The program requires the calculator to be in RPN mode. Place the coefficients of the quadratic on the stack and execute the program. The answers are displayed in the x and y registers. For example, to use the program to find the roots of 3x^2 + 2x -1 = 0 press [3] [ENTER] [2] [ENTER] [1] [+/-] [XEQ] [Q] [ENTER]. After execution of this example, y = 0.33333... and x = -1.

I realize there are much shorter quadratic root finder programs available for the HP 35s. The goal of this program is to handle all types of cases with accuracy, both real and complex, including cases that challenge the accuracy and range limits of the calculator. It also leaves the "last x" register unaffected and has an undo feature to remove the result from the stack and return the original coefficients to the stack. To undo, press [XEQ] [Q] [0] [0] [3].

The program uses storage registers H, O, and U as temporary storage. It leaves the "user flags" unaffected. (It manipulates flag 10 but will leave the flag as it was when the program started.)

Try some of the examples shown in the attached PDF file on other quadratic root finder programs to see if they perform as accurately as this program does.

This program uses algorithms adapted from this journal article:
R.D. Middlebrook, "Methods of Design-Oriented Analysis: The Quadratic Equation Revisited," in Frontiers in Education, 1992 Proceedings. Twenty-Second Annual Conference, pages 95-101, published by IEEE.

The heart of these methods is also contained in this PDF:
http://www.rdmiddlebrook.com/D_OA_Rules&...0Roots.pdf

[Jeff_Kearns]

(05-30-2014 02:29 AM)Douglas De Boer Wrote:  The goal of this program is to handle all types of cases with accuracy, both real and complex, including cases that challenge the accuracy and range limits of the calculator.

I haven't entered the 123 line program yet, however it would be interesting if you could compare its performance with this 95 line "Cadillac" Quadratic Solver by the late Palmer O. Hanson Jr. for the HP-35s.

It solves the following 8 cases 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

Regards,

Jeff K

[Douglas De Boer]

Quote:Case 1:

a = 1E-13 b = -2 c = 1

R1 = 2E13 R2 = 0.5

Jeff, Thanks for pointing out this "Cadillac" program. The first case is solved incorrectly by my program--an outright bug. My program gives R2 = 0. I'll look into it and fix that when I find time.

Most of the other cases show some discrepancies compared to my program. These look like matters of numerical accuracy. Since the "Cadillac" program has been around for a while, I think it is fair to trust that one rather than mine at this time. I'll keep looking at this, mainly since it is an interesting puzzle for me.

Thanks for your comments!

Regards,

--DDB

[Jeff_Kearns]

(05-30-2014 04:15 AM)Douglas De Boer Wrote:  Thanks for your comments!
--DDB

DDB,

The Quadratic Solution is also an (amateur) interest of mine, especially as it pertains to the complex plane. I highly recommend the following math videos called Calculus Revisited where this impressive lecturer Herbert Gross explains complex numbers in a very interesting manner.

A great thread on this forum related to which quadratic solution should be used can be found here. There are also some neat (and short) Quadratic Equation solvers for the 42s, 32sii and 33s on this site if you do a search. There is also a short HP-15C routine by Gerson Barbosa and Thomas Klemm:

001- LBL B
002- ENTER
003- R^
004- /
005- R^
006- LSTx
007- /
008- 2
009- CHS
010- /
011- ENTER
012- ENTER
013- x^2
014- R^
015- -
016- SQRT
017- -
018- x<>y
019- LSTx
020- +
021- RTN

Usage: To solve for 'x' in Ax^2+Bx+C=0

Set Flag 8

1) Enter coefficients: A [ENTER] B [ENTER] C
2) GSB B
3) Roots appear in X and Y registers
4) If there are complex roots (and Flag 8 is not set), result is Error 0
5) Set Flag 8 and rerun
6) First root is in X register. Real part and imaginary parts can be seen by toggling Re<>Im. To see second root, press X<>Y and toggle Re<>Im

From the HP-15C Advanced Functions Handbook, page 208:

"Programs do exist which, while carrying only 10 significant digits during arithmetic, will calculate the roots of any quadratic correctly to at least nine significant digits regardless of how nearly coincident those roots may be. All such programs calculate
d = b2 - ac by some trick tantamount to carrying 20 significant digits whenever b2 and ac nearly cancel, so those programs are a lot longer and slower than the simple subroutine provided above."

The manual includes an 82 line program that accomplishes just that, likened to "Grandmother's expensive chinaware, reserved for special occasions, leaving subroutine "A" for everyday use."

Long live the HP-15C!

Jeff K

[Douglas De Boer]

Jeff,

The bug in my program is that in ax^2 + bx +c = 0, the b must be non-negative in Middlebrook's technique, which I ignored. There is no loss of generality because if b < 0 then just multiply the whole equation through by -1. That does not change the roots. Taking this into account improves the results considerably on my program.

But there are more issues, which I cannot exactly put my finger on at this moment. Middlebook's technique still seems to suffer from an inaccurate calculation of the discriminant in some cases. I'll keep thinking about it. And I think the next thing I'd like to do is try some of the other programs myself, including the "Cadillac" program.

Thanks for all the enlightening help.

Regards,

--DDB

[W.B.Grant]

Thank you both. This has been most interesting and very informative.
Best,
WBG

[Thomas Klemm]

(05-30-2014 06:01 PM)Douglas De Boer Wrote:  And I think the next thing I'd like to do is try some of the other programs myself, including the "Cadillac" program.

For those interested on how the "Cadillac" solver works:
Solve the real quadratic equation \(c−2bz+az^2=0\) for real or complex roots.


Cheers
Thomas