Short Quadratic Solver (HP-42S)



#37

Here is a 42S-specific quadratic solver. It is 28-byte long, uses only the stack and preserves the original stack register X. I have also a 27-byte version, but it requires one extra step.

QUADRATIC SOLVER - HP-42S

X Y Z T L
00 { 28-Byte Prgm }
01 LBL "Q" c b a T L
02 RCL/ ST Z c/a b a T c
03 X<> ST Z a b c/a T c
04 ST0+ ST X 2a b c/a T c
05 / b/(2a) c/a T T 2a
06 +/- -b/(2a) c/a T T 2a
07 STO ST Z -b/(2a) c/a -b/(2a) T 2a
08 x^2 b^2/(4a^2) c/a -b/(2a) T -b/(2a)
09 X<>Y c/a b^2/(4a^2) -b/(2a) T -b/(2a)
10 - b^2/(4a^2)-c/a -b/(2a) T T c/a
11 SQRT sqrt(b^2/(4a^2)-c/a) -b/(2a) T T b^2/(4a^2)-c/a
12 STO ST L sqrt(b^2/(4a^2)-c/a) -b/(2a) T T sqrt(b^2/(4a^2)-c/a)
13 X<>Y -b/(2a) sqrt(b^2/(4a^2)-c/a) T T sqrt(b^2/(4a^2)-c/a)
14 STO+ ST Y -b/(2a) -b/(2a)+sqrt(b^2/(4a^2)-c/a) T T sqrt(b^2/(4a^2)-c/a)
15 RCL- ST L -b/(2a)-sqrt(b^2/(4a^2)-c/a) -b/(2a)+sqrt(b^2/(4a^2)-c/a) T T -b/(2a)
16 END

Examples:

1) Given x2 - 5x + 6 = 0, compute 123456*x1 and 123456*x2

123456 ENTER 1 ENTER 5 +/- ENTER 6 XEQ Q -> Y: 3 ; x2 = 3
X: 2 ; x1 = 2
Rv * -> X: 370,368 ; 123456*x2
Rv * -> X: 246,912 ; 123456*x1

2) 3x2 - 12x + 87 = 0

3 ENTER 12 +/- ENTER 87 XEQ Q -> Y: 2 i5 ; x2 = 2 + 5i
X: 2 -i5 ; x1 = 2 - 5i

2) x2 - (4 + 6i)x + (-5 + 10i) = 0

1 ENTER 4 ENTER 6 COMPLEX +/- 5 +/- ENTER 10 COMPLEX XEQ Q -> Y: 3 i4 ; x2 = 3 + 4i
X: 1 i2 ; x1 = 1 + 2i


#38

Slightly different, but still 28-byte and 16-step long:

QUADRATIC SOLVER - HP-42S

00 { 28-Byte Prgm }
01 LBL "Q"
02 RCL/ ST Z
03 X<> ST Z
04 /
05 2
06 /
07 STO ST Z
08 x^2
09 X<>Y
10 -
11 SQRT
12 RCL- ST Y
13 X<>Y
14 RCL+ ST L
15 +/-
16 END


#39

One step shorter:

00 { 28-Byte Prgm }                        
01 LBL "Q"
02 RCL/ ST Z
03 X<> ST Z
04 /
05 -2
06 /
07 STO ST Z
08 x^2
09 X<>Y
10 -
11 SQRT
12 RCL+ ST Y
13 X<>Y
14 RCL- ST L
15 END

#40

From my archives ;-)
One byte shorter, more accurate for cases where b^2 is large compared to 4*a*c, and compatible with the 41C:

 x := -b/(2*a);
y := c/a ;
d := SQRT(x^2 - y);

r1 := x + sign(x)*d;
r2 := y/r1;

00 { 27-Byte Prgm } X Y Z T
01*LBL "Q" c b a
02 X<> ST Z a b c
03 STO/ ST Z a b c/a
04 STO+ ST X 2*a b c/a
05 +/- -2*a b c/a
06 / x y
07 ENTER x x y
08 X^2 x^2 x y
09 RCL ST Z y x^2 x y
10 - x^2-y x y y
11 SQRT d x y y
12 RCL ST Y x d x y
13 SIGN +/-1 d x y
14 * +/-d x y
15 + r1 y y y
16 STO ST Z
17 /
18 END


#41

Hello Werner,

Quote:
One byte shorter, more accurate for cases where b^2 is large compared to 4*a*c, and compatible with the 41C

These are valuable features. Michael Harwood's program in the Software Library (Fast Quadratic Formula for the HP-41C) is 16-step long (when the prompt is removed) and would require 25 bytes on the 42S, but it may not be as accurate as yours. I have been trying to preserve the stack register X. Without this constraint, I suspect a 42S-specific program might break the 15-step length barrier.

Regards,

Gerson.


#42

Maybe the HP's SOLVE algorithm is accurate?

Here is my solution (from my memories...) - not shortest, uses variables, but not uses the classical quadratic formula. The 32SII (low-end) version:

--- The firs root:

I01 LBL I
I02 CF 0
I03 FN=Q
I04 SOLVE X
I05 GTO X
I06 SF 0
I07 RTN

7 steps / CK=368D / 10.5 byte


--- The second:

X01 LBL X
X02 STOP
X03 RCL B
X04 RCL div A
X05 +
X06 +/-
X07 RTN

7 steps / CK=0339 / 10.5 byte


--- The quadratic equation:

Q01 LBL Q
Q02 RCL X
Q03 RCL mul A
Q04 RCL add B
Q05 RCL mul X
Q06 RCL add C
Q07 RTN

7 steps / CK=8AAE / 10.5 byte

'add', 'sub', 'mul' and 'div' means 'add', 'substact', 'multiple' and 'divide' in this order.

Solve for real roots of an 'a*x^2+b*x+c=0' equation:
a STO A, b STO B, c STO C then XEQ I.
If FLAG 0 set this equation has no real roots (location of extremum on display).
If FLAG 0 not set, the first real root appear on LCD, then press R/S, then the second root will be on the display.

Enjoy!


#43

Hello Csaba,

Quote:
Maybe the HP's SOLVE algorithm is accurate?

It's quite accurate, as we know.

The idea is getting the shortest possible HP-42S quadratic solver. By short I mean less number of steps (I know on the 42S less bytes would be more appropriate, but I guess less steps means faster execution - to be sure, I would have to check out a table of T-states per instruction, but I don't have any).

Anyway, you approach using the SOLVER is interesting, thanks for posting it. When I tested it on the 33s I added the lines below to make it even easier to use, but then I noticed you had done that already here :-)

...
I0002 INPUT A
I0002 INPUT B
I0003 INPUT C
...

I expected the equivalent 42S version to handle complex cases, but it didn't. Perhaps I just didn't implement it right...

Best regards,

Gerson.

#44

Werner, that's the 27-byte version I had mentioned. It is one step shorter and will solve all my three examples as stated. Neither HP-41 compatible nor accurate enough for the case you handle though. Doing it using only HP-41 instructions would be quite challenging.

QUADRATIC SOLVER - HP-42S

00 { 27-Byte Prgm }
01 LBL "Q"
02 RCL/ ST Z
03 Rv
04 X<>Y
05 STO+ ST X
06 /
07 +/-
08 ENTER
09 X^2
10 RCL- ST T
11 SQRT
12 X<>Y
13 ENTER
14 X<> ST Z
15 STO+ ST Z
16 -
17 END

Regards,

Gerson.

#45

Strange how you can always improve upon an old program when you look at it again:

00 { 25-Byte Prgm }	X	Y	Z	T
01*LBL"Q" c b a
02 X<> Z a b c
03 ST/ Z a b y
04 ENTER a a b y
05 + 2*a b c/a c/a
06 CHS
07 / x y y y
08 ENTER
09 X^2 x^2 x y y
10 RUP
11 -
12 SQRT d x y y
13 RCL Y
14 SIGN
15 *
16 +
17 ST/ Y
18 END

To be sure, it exhibits the same pro's and con's as the previous:

+ 41C compatible
+ avoids cancellation
- no complex support (doesn't work for complex roots, nor for complex arguments)
- doesn't work for a=0 or b,c=0
#46

Wow, nice job!

#47

Very nicely done!

- Pauli


#48

Does it also work with complex numbers? I guess so, since the 42 can hold them in the normal stack?

Reminds me those contests for the shortest implementation of a given program, and I'd have to say the following takes the prize:

01 QROOT

(see SandMath module - of course that's cheating, as the MCODE listing is longer than that :)


#49

But when are you doing MCODE analytic solutions to higher order polynomials? :-)


- Pauli


#50

Touché :)

#51

Hello Ángel,

Quote:
Does it also work with complex numbers? I guess so, since the 42 can hold them in the normal stack?

Yes, as in examples #2 and #3 above.

I've been using this formula:

(Presented by Allen in this thread (Message #2).

The HP-15C handle complex numbers also, as you know:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.cgi?read=114345

Regards,

Gerson.


#52

Since the 41 lacks the "native" ability to handle complex numbers on its stack, I added the following short routine to the 41Z to solve the quadratic roots in such case:

1   LBL "ZQRT" 
2 ZENTER^
3 ZR^
4 Z/
5 LASTZ
6 ZR^
7 Z<>W
8 Z/
9 ZHALF
10 ZNEG
11 ZENTER^
12 ZENTER^
13 Z^2
14 ZR^
15 Z-
16 ZSQRT
17 ZENTER^
18 ZNEG
19 ZR^
20 Z+
21 ZRDN
22 Z+
23 ZRUP
24 END
just 24 steps (or Z-steps, rather). It takes the three coefs in complex stack leves 1,2,3, leaving both solutions on complex stack levels 1 and 2.

Cheers,
ÁM


Edited: 22 Oct 2010, 8:25 a.m.

#53

Thanks Paul!

I hope the 34s has a built-in QUAD command to do it in only one step :-)


#54

The thought had occurred to me to drop your code in as a built in function :-)

- Pauli


#55

Accuracy-wise, Werner's code is better though. For instance, for a = 1E-13, b = -2 and c = 1 it returns the correct 12-digit results, that is 2E13 and 0.5. Mine returns 2E13 and 0. Unless, of course, you want to preserve at least the stack register X.

Gerson.


#56

Werner's code fails for b=c=0 BTW.

Easy enough to detect and return the correct answer for but that means more steps.


- Pauli


#57

It also fails for example #3, I noticed that soon after I wrote. I had checked only the second example. But this can be fixed. By the way, I never needed to solve a quadratic equation with complex coefficients before. I simply expanded (x-(1+2i))*(x-(3+4i)) and there it was.

Regards,

Gerson.

typo correction


Edited: 22 Oct 2010, 10:10 p.m.

#58

Or alternatively, drop the MCODE (oh if only CPU's were compatible :)

094	"T"	
00F "O"
00F "O"
012 "R"
011 "Q"
078 READ 1(Z)
361 ?NC XQ (includes SETDEC)
050 ->14D8 [CHK_NO_S]
0B8 READ 2(Y)
10E A=C ALL
0F8 READ 3(X)
351 ?NC GO (includes SETDEC)
052 ->14D4 [CHK_NO_S1]
2BE C=-C-1 MS -c
10E A=C ALL
078 READ 1(Z) a
261 ?NC XQ
060 ->1898 [DV2_10]
128 WRIT 4(L) -c/a
078 READ 1(Z) a
10E A=C ALL
01D ?NC XQ Adds normalized A and C
060 ->1807 [AD2_10]
10E A=C ALL 2a
0B8 READ 2(Y) b
2BE C=-C-1 MS -b
0AE A<>C ALL
261 ?NC XQ
060 ->1898 [DV2_10]
0E8 WRIT 3(X) register re-use…
10E A=C ALL -b/2a
135 ?NC XQ Multiplies A times C
060 ->184D [MP2_10]
10E A=C ALL (b/2a)^2
138 READ 4(L) -c/a
01D ?NC XQ Adds normalized A and C
060 ->1807 [AD2_10]
244 CLRF 9 (b/2a)^2-c/a
2FE ?C#0 MS is it negative?
01B JNC +03
05E C=0 MS Sign change (ABS)
248 SETF 9 flag as complex
2F9 ?NC XQ
060 ->18BE [SQR10]
24C ?FSET 9 Complex roots?
02B JNC +05
0A8 WRIT 2(Y) Write Im(Z) in Y
04E C=0 ALL Zero Z so we know they're
068 WRIT 1(Z) Conjugated Complex Roots
073 JNC +14d
128 WRIT 4(L) sqr((b/2a)^2-c/a)
10E A=C ALL
0F8 READ 3(X) -b/2a
01D ?NC XQ Adds normalized A and C
060 ->1807 [AD2_10]
0A8 WRIT 2(Y) x1 = -b/2a + sqrt()
138 READ 4(L) sqr((b/2a)^2-c/a)
2BE C=-C-1 MS
10E A=C ALL
0F8 READ 3(X) -b/2a
01D ?NC XQ Adds normalized A and C
060 ->1807 [AD2_10]
0E8 WRIT 3(X) x2 = -b/2a - sqrt()
3E0 RTN

#59

Input of a, b, c as z, y, x and output of possibly complex x1 and x2 in x, y, z, and t ? Do we need more reasons for an 8 level stack? ;)

#60

Very nice!

You are to be highly commended for the high level of documentation of the process!! This is all too often negelected, both in calculator and computer algorithms.

TomC


#61

Thanks! However the documentation is still incomplete. The formulas used in the first version are, as we have seen,

'x1=-b/(2*a)+sqrt((b/(2*a))^2-c/a)'

'x2=-b/(2*a)-sqrt((b/(2*a))^2-c/a)'

which is the quadratic formula for

'x^2-b/a*x+c/a=0'
In the second and third versions, an attempt has been made to carry this one step further, by making the inputs coefficients a, b and c equal to 1, b/a and c/a, respectively. This requires additional steps, but the formulas become somewhat simpler:
'x1=b/2+sqrt((b/2)^2-c'

'x2=b/2-sqrt((b/2)^2-c'

Regards,

Gerson.


#62

Gerson, I had done this many many years ago for, I think, my HP-11c, but I believe mine was slightly different. (Actually, looking again, it's the same.)

Start with the usual ax^2 + bx + c = 0. Then define:

m = -b/(2a) and n = c/a

The solutions become

x1 = m + sqrt( m^2 - n)

x2 = m - sqrt(m^2 - n)

I can't remember if I used the stack or storage registers. Usually the upfront calculations of m and n are trivial, and allows for a very compact program, and no fractions in the final calculations.

For something like: 3x^2+5x-8=0 I'd key -5/6 run LBL A, -8/3, press R/S and get the solutions. I forget the actual program length, but it was quite small. I'll have to dig out the 11C and try it again.

CHUCK


#63

Chuck,

This appears quite suitable for hand calculations also, as the formula is easier to use and remember. My 14-year old daughter has been studying quadratic equations and application problems in school, but so far she's been taught the traditional method only.

Regards,

Gerson.


#64

The "classic" formula is still taught in high schools (also on my kids' case), perhaps on the grounds of being easier to memorize (?) - yet the other one is much more suitable for programming. That's the one I used in the MCODE listing, it saves time and storage space.

Cheers

#65

This program returns only one of the possible solutions of the equation x2 + px + q = 0. The quadratic solver hidden in ACOSH is used:

00 { 17-Byte Prgm }
01>LBL "Q"
02 SQRT
03 ÷
04 LASTX
05 X<>Y
06 -2
07 ÷
08 ACOSH
09 E^X
10 ×
11 END

With a minor change (-2 -> 2 CHS) it may be used with the HP-15C as well. However it won't work with the HP-35s due to its lack of proper support of complex functions.

Usage:

1) x2 - 5x + 6 = 0
5 +/- ENTER 6 XEQ Q -> X: 3

2) x2 - 4x + 29 = 0
4 +/- ENTER 29 XEQ Q -> X: 2 i5

3) x2 - (4 + 6i)x + (-5 + 10i) = 0
4 ENTER 6 COMPLEX +/- 5 +/- ENTER 10 COMPLEX XEQ Q -> X: 3 i4

Since neither addition nor subtraction is used there's no loss
of significance due to cancellation.

Best regards

Thomas

PS: Finding the other solution is left as an exercise.


Edited: 23 Oct 2010, 2:41 a.m. after one or more responses were posted


#66

Very nice idea - now we're shortening the program listing although at the expense of using other functions - a valid rule in the "brevity contest" :-)

What about using P-R instead of ACOSH for the 41?


#67

This could only be used when q <= 0. But even then you'd have to calculate the square root of -q previously to use R-P. So I'm afraid we won't gain anything. Only now I see you mentioned P-R, not R-P. What did you have in mind?

Edited: 23 Oct 2010, 11:02 a.m.

#68

Okay the following works only for certain values of p and q:

0 < 4q < p2

00 { 20-Byte Prgm }
01>LBL "Q"
02 SQRT
03 ÷
04 LASTX
05 X<>Y
06 -2
07 X<>Y
08 ÷
09 ASIN
10 2
11 ÷
12 TAN
13 ÷
14 END

0 < -4q < p2

00 { 21-Byte Prgm }
01>LBL "Q"
02 +/-
03 SQRT
04 ÷
05 LASTX
06 X<>Y
07 -2
08 X<>Y
09 ÷
10 ATAN
11 2
12 ÷
13 TAN
14 ÷
15 END

So all we have to do is combine these two programs into one.

BTW: I used here the following formulas:

                 2
sin(2x) = ---------------
cot(x) + tan(x)


2
tan(2x) = ---------------
cot(x) - tan(x)

Cheers

Thomas


Edited: 23 Oct 2010, 7:15 a.m.

#69

Hi Thomas,

Quote:
PS: Finding the other solution is left as an exercise.

Still unoptimized, as I'm in a hurry:

00 { 25-Byte Prgm }
01>LBL "Q"
02 X<>Y
03 STO Z
04 X<>Y
05 SQRT
06 ÷
07 LASTX
08 X<>Y
09 -2
10 ÷
11 ACOSH
12 E^X
13 ×
14 ENTER
15 R^
16 +
17 +/-
18 END

Regards,

Gerson.


#70

A little bit shorter:

00 { 21-Byte Prgm }
01>LBL "Q"
02 STO ST Z
03 SQRT
04 ÷
05 LASTX
06 X<>Y
07 -2
08 ÷
09 ACOSH
10 E^X
11 ×
12 ÷
13 LASTX
14 END

However I had something else in my mind.


Edited: 23 Oct 2010, 7:51 a.m.


#71

You're right! x2 = q/x1, for x1 <> 0, is shorter than x2 = -(x1+p).
This can be made shorter, I think.

Thanks for presenting your interesting insight.

Regards,

Gerson.

00 { 31-Byte Prgm }
01>LBL "Q"
02 RCL/ ST Z
03 X<> ST Z
04 /
05 X<>Y
06 STO ST Z
07 SQRT
08 ÷
09 -2
10 ÷
11 ACOSH
12 E^X
13 RCL ST Y
14 SQRT
15 ×
16 X<>Y
17 RCL/ ST Y
18 END

It does work for my three examples:

1) Given x2 - 5x + 6 = 0, compute 123456*x1 and 123456*x2 

123456 ENTER 1 ENTER 5 +/- ENTER 6 XEQ Q -> Y: 2.99999999999 ; x2 = 3
X: 2.00000000001 ; x1 = 2
Rv * -> X: 370,367.999999 ; 123456*x2
Rv * -> X: 246,912.000001 ; 123456*x1

2) 3x2 - 12x + 87 = 0

3 ENTER 12 +/- ENTER 87 XEQ Q -> Y: 2.00E0 i5.00E0 ; x2 = 2 + 5i
X: 2.00000000001 -i5 ; x1 = 2 - 5i

3) x2 - (4 + 6i)x + (-5 + 10i) = 0

1 ENTER 4 ENTER 6 COMPLEX +/- 5 +/- ENTER 10 COMPLEX XEQ Q -> Y: 3 i3.99999999999 ; x2 = 3 + 4i
X: 1.00E0 i2.00E0 ; x1 = 1 + 2i

It will work for a=1E-13, b=-2 and c=1 also:

1.99999999996E13
0.50000000001

#72

It seems to be a recurring theme that we discuss every few years. Thanks for starting that thread (again).

Best regards

Thomas


Possibly Related Threads...
Thread Author Replies Views Last Post
  hp-prime solver and variable name fabrice48 22 1,532 12-10-2013, 03:25 AM
Last Post: fabrice48
  Solver issue with HP 17BII - different from 19BII Jeff Kearns 13 756 11-28-2013, 02:36 AM
Last Post: Don Shepherd
  HP Prime Triangle solver BruceH 29 1,742 11-28-2013, 12:03 AM
Last Post: Dale Reed
  HP Prime - Short "learning" modules CR Haeger 1 254 11-27-2013, 02:13 PM
Last Post: Jonathan Cameron
  HP prime: linear solver app Alberto Candel 1 307 11-21-2013, 01:57 AM
Last Post: Michael Carey
  I have written a short introduction to the HP Prime Michael Carey 7 592 11-18-2013, 08:04 PM
Last Post: Michael Carey
  HP Prime: Linear Solver app bug BruceH 0 204 11-15-2013, 06:36 PM
Last Post: BruceH
  HP-65 short circuit Ignacio Sánchez 2 314 10-22-2013, 08:27 AM
Last Post: Ignacio Sánchez Reig
  HP Prime Solver Variables Issue Anibal Morones Ruelas 8 658 10-19-2013, 09:45 AM
Last Post: Harold A Climer
  HP Prime triangle solver oddity BruceH 0 194 10-13-2013, 09:08 PM
Last Post: BruceH

Forum Jump: