Which quadratic solution should we use?



#8

A "more efficient and accurate" method for the solution of quadratic equations appeared in the HP-67 Standard Pac. A similar method which "avoids destructive cancellation" appeared in the HP-15C Advanced Functions Handbook. Those "better" methods have not been consistently used in later H-P applications.

Is there a reason for not using the "more efficient and accurate" methods?

Is there an earlier reference on the "more efficient and accurate" methods than that in the HP-67 Standard Pac?


Avoiding Destructive Cancellation in the Solution of a Quadratic Equation

Page 31 of Maurice D. Weir's Calculus by Calculator (Prentice-Hall, 1982) gives the typical algebra textbook solution for a quadratic equation but then states:

"... The program QUADS calculates the real roots (if they exist) for any specified quadratic equation. The formula used is a modification of Equation 1.11 [the textbook form] to reduce roundoff error."

Analysis of the QUADS program reveals that the solution sequence is:

1. Calculate the determinant D = B2 - 4AC

2. Calculate the value E = (sqr(D) + abs(B)/2A

3. Change the sign of E if B<0

4. Calculate the real roots as C/AE and E.

A more thorough treatment of the subject with a similar solution method appears in the appendix "Accuracy of Numerical Calculations" to the HP-15C Advanced Functions Handbookalso published in 1982. The discussion on page 191 presents the textbook formula and an algebraically equvalent calculation sequence which avoids "destructive cancellation". (Readers of that reference should note that for some reason the authors chose to define the quadratic as c - 2bz + az2 . The more common definition is used here.) The discussion states

"... Such a program will be listed later (page 205) and must be used in those instances, common in engineering, when the smaller root y is needed accurately despite the fact that the quadratic's other unwanted root is relatively large."

A demonstration problem is provided. For a quadratic equation of the form ax2 + bx + c and a = 1E-13, b = -2 and c = 1 the exact solution is r1 = 2E13 and r2 = 0.5. Implementation of the textbook formula on most HP and TI machines will yield r1 = 2E13 and r2 = 0. .

The rest of this submission presents all the gory details which led up to this submission.


Some H-P History on a Better Quadratic Solution Methoid

This history is limited to the hardware and documentation available from my own collection.

The method was not included in the quadratic solution in the HP-65 Standard Pac printed in 1974 or in the quadratic solution in the HP-27 Owner's Handbook printed in 1976.

Pages 198 through 203 of the HP-67 Owner's Handbook and Programming Guide (Revision C 12/76) presented a quadratic solver program which used the textbook formula. The end of the discussion includes the following statement "For a more efficient and accurate method of finding the roots of a quadratic equation, see the Polynomial Evaluation program in your HP-67 Standard Pac."

The Polynomial Evaluation program on pages 9-01 though 9-04 of the HP-67 Standard Pac (1976) uses a routine which addresses destructive cancellation. It solves the test case successfully.

Pages 61 through 65 of the HP-33E/33C Owner's Handbook and Programming Guide (April 1979) presents a quadratic solver program which uses the textbook formula. It can not be expected to and does not solve the test case successfully.

Pages 179 through 183 of my copy of the HP-41C Owner's Handbook and Programming Guide (Revision B 8/79) present a quadratic solver program which uses the textbook formula. The program cannot solve the test case satisfactorily. A note at the bottom of page 179 does offer the following caution "Some values of a, b, and c may result in misleading answers because their solutions require greater than twelve digits of accuracy."

The "Polynomial Solutions/Evaluation" program of the HP-41 Math Pac (11/79) includes a quadratic solver which uses the textbook formula. It cannot solve the test case satisfactorily.

Pages 206 through 210 of the HP-11C Owner's Handbook and Problem Solving Guide (Revision B 9/81) describes how to develop a quadratic solver using the textbook formula.

A thorough discussion of the method for avoiding destructive cancellation appears on pages 191 and 205-211 of the HP-15C Advanced Functions Handbook which was published in 1982. A program for the HP-15C is included which will solve the test case satisfactorily.

The QUAD option of the SOLVE menu of the HP-28S (1988) does not solve the test problem correctly.

Pages 172 through 174 of my copy of the HP-42S RPN Scientific Calculator Owner's Manual (10/88) present a demonstration of the conversion of a quadratic solver program from the HP-41CV Owner's Manual for use on the HP-42S. The solution uses the textbook formula and will not solve the test case satisfactorily..


Some TI History on A Better Quadratic Solution Method

The "Solution of Quadratic Equation" in the SR-56 Application Library (1976) and page IV-82 of Personal Programming (the Owner's Manual for the TI-58/59, 1977) only present quadratic solver programs using the textbook formula which will not solve the test case satidfactorily..

A program for the TI-59 by Bill Skillman which uses a method which will solve the test case satisfactorily was published in late 1978 in V3N12P2 of 52 Notes.

My article "Improved Solution for the Quadratic" on pages 8-12 of TI PPC Notes (1991) provides a solution for the TI-59 which uses the technique in Calculus by Calculator. It will yield the correct answer to the test case.

The POLY routine in the TI-68 (1989), the QAD routine of the TI-95 (1986) and the POLY routines in the TI-85 and TI-86 all yield the correct answer to the test case but there is no discussion of the use of an improved method in the manuals. The TI-89 yields inconsistent results. The Polynomial Root Finder of the Apps Desktop returns roots of 2.E13 and 0.51. However, if the test case is entered through the home screen as either a" factor" problem or as a "solve" problem a correct answer will be returned.


hp-33s Results

The following table presents results obtained with two different solution methods on the hp-33s. The Textbook solution is from the Polynomial Root Finder program on pages 15-20 through 15-32 of the hp-33s User's Guide. The Revised solution is an implementation of the technique for avoidance of destructive cancellation defined in the HP-15C Advanced Functions Handbook. A listing appears at the end of this submission. The coefficients are from test cases 3 through 6 on page 207 of the HP-15C Advanced Functions Handbook but with the first degree coefficient multiplied by minus two to change the definition of a quadratic from ax2 - 2bx + c to ax2 + bx + c.

Coefficients             Correct              Textbook                Revised

a 1E-13 2E13 2E13 2E13
b -2 0.5 0 0.5
c 1

a 654323 1 1.000000921 1.000000000
b -1308644 0.9999969434 0.999996022 0.999996943
c 654321

a 11713 62.77179203 62.771792026 62.771792026
b -1470492 8.5375E-05 i 62.771792026 8.5375E-05 i
c 46152709

a 80841 12.2171175 12.217117552 12.217117552
b -1975288 1.374514E-03 i 1.374045E-03 i 1.374514E-03 i
c 12066163

HP-41 Results

The following table presents results obtained with three different solution methods on the HP-41. The Math-Pac solution uses a version of the textbook formula but with the coefficient of the second degree term set to one. The Textbook solution uses the HP-41 program from pages 172 through 174 of the HP-42S RPN Scientific Calculator Owner's Manual modified to complete the solution for complex roots. The Revised solution uses a program translated from the Revised hp-33s solution in this submission.

   Correct                Math-Pac                  Textbook               Revised

2E13 2E13 2E13 2E13
0.5 0 0 0.5

1 1.000000 0.999998472 0.999998472
0.9999969434 2.0000E-05 i 0.999998472 0.999998472

62.77179203 62.77352410 62.77179203 62.77179203
8.5375E-05 i 62.77006000 62.77179203 62.77179203

12.2171175 12.21711755 12.21711755 12.21711755
1.374514E-03 i 1.414214E-03 i 1.369104E-03 i 1.377460537E-03

The curious result from the Math-Pac for the second problem results from the definition of the coefficient of the second degree term as unity. That means that the coefficient of the first degee term is entered as 654321/654323 and the constant term as -1308644/654323. I obtain a similar curious result if I enter those coefficients in either the textbook solution or in the revised solution.


An hp-33s Program Which Avoids Destructive Cancellation

K0001  LBL K
K0002 CF 1
K0003 INPUT A
K0004 INPUT B
K0005 INPUT C
K0006 RCL A
K0007 x ne 0 ? not equal to
K0008 GTO L
K0009 RD roll down
K0010 x<>y
K0011 /
K0012 +/-
K0013 0
K0014 1/x
K0015 STOP
K0016 GTO K
L0001 LBL L
L0002 x
L0003 x<>y
L0004 +/-
L0005 2
L0006 /
L0007 ENTER
L0008 RD
L0009 x^2 x squared
L0010 -
L0011 x<=0?
L0012 GTO M
L0013 SF 1
L0014 square root
L0015 x<>y
L0016 /
L0017 LASTx
L0018 R^ roll up
L0019 x<>y
L0020 /
L0021 STOP
L0022 GTO K
M0001 LBL M
M0002 +/-
M0003 square root
M0004 x<>y
M0005 RD
M0006 x<>y
M0007 SGN
M0008 x
M0009 +
M0010 ENTER
M0011 R^
M0012 /
M0013 x<>yt
M0014 x=0?
M0015 STOP
M0016 x=0?
M0017 GTO K
M0018 RCL C
M0019 x<>y
M0020 /
M0021 STOP
M0022 GTO K


#9

I have been reading through Wickes's book "HP 48 insights" that another forum member sent. It is truly beautiful to see in practice what Karl mentioned a few days ago:

Quote:
...use equation-based programming when the formula is simple and straightforward, and use keystroke programming when the formula or calculation is lengthy...

Is not the quadratic equation a worthy example? I did not appreciate the simplicity of the equation until I saw it in this form:



There are really only two terms here: c/a and b/2a, the second of which comes up twice! So, a very straightforward RPN-equation would be a two-fold process:

[a b c] -> [c/a b/2a] -> [root1 root2] 



I can't find my original program from 1993, so here is a rough hack for one at 52.5 bytes. I know someone can get one below 50 bytes.
<< 
SWAP -2 / ->V2 SWAP / V-> 'sets up [c/a,-b/2a] using vector
DUP SQ ROT - (SQRROOT) 'finds SQRROOT(b^2/4a^2 -c/a)
DUP2 + 3 ROLLD - 'gets both Y+X and Y-X on the stack.
>>

Then again you could just type in:
[pre]
<<
->V3 PROOT
>>

But I think that would defeat the purpose....


Edited: 4 Mar 2007, 1:53 a.m.


#10

This is very elegant but it still has destructive cancellation and gives wrong answers for all of the tough examples Palmer gives.

Les

#11

Hi, Palmer --

Quote:
A more thorough treatment of the subject with a similar solution method appears in the appendix "Accuracy of Numerical Calculations" to the HP-15C Advanced Functions Handbook also published in 1982. The discussion on page 191 presents the textbook formula and an algebraically equvalent calculation sequence which avoids "destructive cancellation". (Readers of that reference should note that for some reason the authors chose to define the quadratic as c - 2bz + az2. The more common definition is used here.) The discussion states

"... Such a program will be listed later (page 205) and must be used in those instances, common in engineering, when the smaller root y is needed accurately despite the fact that the quadratic's other unwanted root is relatively large."


One advantage of defining the quadratic that way is that the determinant becomes b2 - ac instead of b2 - 4ac, and the magnitude of those terms is reduced by a factor of 4.

I haven't fully examined the program to determine exactly how it works. Both the b2 and a*c terms overflow a 10-digit register in cases 4, 5, and 6. "b2 - ac" nonethless is calculated and tested against zero in several ways.

-- KS

Quote:
A demonstration problem is provided. For a quadratic equation of the form ax2 + bx + c and a = 1E-13, b = -2 and c = 1 the exact solution is r1 = 2E13 and r2 = 0.5. Implementation of the textbook formula on most HP and TI machines will yield r1 = 2E13 and r2 = 0. .

I get

r1 = 19999999999999.4999999999999875     = 2E13 - 0.5 - 1.25E-14
r2 = 0.50000000000001249999999996826535 ~= 0.5 + 1.25E-14

on Windows XP Calculator, with its 80-bit extended double-precision variables. These results are quite close to exact, with about 1E-25 calculation roundoff error when re-inserted into the original equation.

x = 0.5 would be the only (and exact) solution if "a" were zero.

x = 2E13 yields ax2 + bx = zero exactly, so ax2 + bx + c = 1.


In the original problem, calculating the square root of the determinant "b2 - 4ac" as

sqrt(abs(b + 2*sqrt(a)*sqrt(c)) * sqrt(abs(b - 2*sqrt(a)*sqrt(c))

roundoff error was reduced, but 12 significant digits isn't enough without using special methods.

-- KS

Edited: 4 Mar 2007, 3:03 p.m.


#12

Quote:

I get

r1 = 19999999999999.4999999999999875     = 2E13 - 0.5 - 1.25E-14
r2 = 0.50000000000001249999999996826535 ~= 0.5 + 1.25E-14

on Windows XP Calculator, with its 80-bit extended double-precision variables. These results are quite close to exact, with about 1E-25 calculation roundoff error when re-inserted into the original equation.


You can, of course, simply do the arithmetic with pencil and paper and find that r1 will be a little less than 2E13 and r2 will be a little more than 0.5. Furthermore, for a = 1E-nn, b = -2 and c = 1 the pencil and paper answers will be r1 a little less than 2Enn and r2 a little more than 0.5 . Some tests with a number of machines show that if nn is less than the number of digits carried by the machine (but in some cases only if nn is two less than the number of digits carried by the machine) then the machine will get the 0.5 answer using the textbook solution. This means that getting the 0.5 answer for r2 from unknown routines on fifteen digit machines like the TI graphic calculators, for example, does not mean that the unknown routines use the more accurate formulation. I have tested the routines in the TI-85, TI-86 and TI-89 to verify that they will get the 0.5 answer when nn is greater than 15.

#13

Thank you, Palmer, for this thorough and enlightening post. To answer your original question, I am guessing HP's logic goes something like this:

Professional engineers etc. who are paying big $$ for the higher-end calcs and reading the advanced-function manuals are likely to run into unusual problems for which they need accurate solutions. Purchasers of simpler, less expensive machines are likely to be students for whom learning to use and program the calculator to solve textbook problems is more useful. HP may have thought that the more accurate but less obvious program would confuse and intimidate beginners and students but be worth the consideration (and extra keystrokes) of more advanced users

John


(Post edited for additional note)

Though it may be defeating some purposes, the PROOT function on the 50g returns the correct result for all 4 of Allen's examples, and does so in less than half a second. Presumably it uses the better algorithms that Palmer describes.


Edited: 4 Mar 2007, 5:48 p.m. after one or more responses were posted


#14

Quote:
I am guessing HP's logic goes something like this:

Professional engineers etc. who are paying big $$ for the higher-end calcs and reading the advanced-function manuals are likely to run into unusual problems for which they need accurate solutions. Purchasers of simpler, less expensive machines are likely to be students for whom learning to use and program the calculator to solve textbook problems is more useful. HP may have thought that the more accurate but less obvious program would confuse and intimidate beginners and students but be worth the consideration (and extra keystrokes) of more advanced users


You may be correct. I find it interesting that H-P chose to give the more accurate formulation for the HP-67 and the HP-15C, albeit not in the baseline documentation but in extended documentation. This may be consistent with the opinion expressed with some frequency in this forum that H-P isn't providing documentation as carefully as they used to.

One concern about the "more accurate" method arises from my background in aerospace and weapons systems where we always worried a little (maybe a lot) when we changed software to fix one problem for fear that we may have inadvertently caused another problem.


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP-PIRME] BUG Pretty Print, Solution: abs => |...|, ABS => ||...|| CompSystems 2 1,036 12-13-2013, 09:36 AM
Last Post: CompSystems
  Quadratic & Cubic Regression RPN progs Matt Agajanian 9 1,410 09-17-2013, 11:37 AM
Last Post: Jeff O.
  [HP-Prime xCAS] Arrays: Programming Commands (solution to the ambiguity) CompSystems 6 1,261 08-30-2013, 06:32 PM
Last Post: CompSystems
  HP-Prime: Solution to ambiguity with [ ... ] & another BUGs CompSystems 12 1,628 08-28-2013, 08:48 AM
Last Post: CompSystems
  Advanced User's Manual and solution to the ambiguity of data types [HP-Prime xCAS] CompSystems 15 2,670 08-20-2013, 03:37 PM
Last Post: Thomas Klemm
  WP34s program submission: Quadratic fit Andrew Nikitin 2 650 06-13-2013, 02:44 AM
Last Post: Paul Dale
  HHC 2012 RPN Programming Contest - 10-Step Way-After-Nashville Solution Jeff O. 32 3,874 10-12-2012, 01:41 AM
Last Post: Paul Dale
  New Quadratic Lagrangian Interpolation Method Namir 2 711 07-20-2012, 04:32 PM
Last Post: Namir
  Fast Quadratic Formula for the HP-41C Gerson W. Barbosa 21 3,079 07-18-2012, 08:53 AM
Last Post: Gerson W. Barbosa
  Quadratic formula program help Chris C 18 2,453 06-16-2012, 12:14 AM
Last Post: Matt Agajanian

Forum Jump: