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

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

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

Very nicely done!

- Pauli

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 :)

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

- Pauli

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.

Touché :)

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

Cheers,
ÁM

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

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

Wow, nice job!

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

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.

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'

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

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

Regards,

Gerson.

Thanks Paul!

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

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

- Pauli

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.

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

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

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.

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	

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.

This program returns only one of the possible solutions of the equation x² + 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

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

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?

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.

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

0 < 4q < p²

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 < p²

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.

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? ;)

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.

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.

You're right! x₂ = q/x₁, for x₁ <> 0, is shorter than x₂ = -(x₁+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

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

Best regards

Thomas

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.

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!

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.

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