▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Here is a 42Sspecific quadratic solver. It is 28byte long, uses only the stack and preserves the original stack register X. I have also a 27byte version, but it requires one extra step.
QUADRATIC SOLVER  HP42S
X Y Z T L
00 { 28Byte 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 x^{2}  5x + 6 = 0, compute 123456*x_{1} and 123456*x_{2}
123456 ENTER 1 ENTER 5 +/ ENTER 6 XEQ Q > Y: 3 ; x_{2} = 3
X: 2 ; x_{1} = 2
Rv * > X: 370,368 ; 123456*x_{2}
Rv * > X: 246,912 ; 123456*x_{1}
2) 3x^{2}  12x + 87 = 0
3 ENTER 12 +/ ENTER 87 XEQ Q > Y: 2 i5 ; x_{2} = 2 + 5i
X: 2 i5 ; x_{1} = 2  5i
2) x^{2}  (4 + 6i)x + (5 + 10i) = 0
1 ENTER 4 ENTER 6 COMPLEX +/ 5 +/ ENTER 10 COMPLEX XEQ Q > Y: 3 i4 ; x_{2} = 3 + 4i
X: 1 i2 ; x_{1} = 1 + 2i
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Slightly different, but still 28byte and 16step long:
QUADRATIC SOLVER  HP42S
00 { 28Byte 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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
One step shorter:
00 { 28Byte 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
▼
Posts: 163
Threads: 7
Joined: Jul 2007
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 { 27Byte 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^2y 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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 HP41C) is 16step 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 42Sspecific program might break the 15step length barrier.
Regards,
Gerson.
▼
Posts: 59
Threads: 9
Joined: Apr 2008
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 (lowend) 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!
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 HP42S 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 Tstates 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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Werner, that's the 27byte version I had mentioned. It is one step shorter and will solve all my three examples as stated. Neither HP41 compatible nor accurate enough for the case you handle though. Doing it using only HP41 instructions would be quite challenging.
QUADRATIC SOLVER  HP42S
00 { 27Byte 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.
Posts: 163
Threads: 7
Joined: Jul 2007
Strange how you can always improve upon an old program when you look at it again:
00 { 25Byte 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
Posts: 764
Threads: 118
Joined: Aug 2007
Posts: 3,229
Threads: 42
Joined: Jul 2006
Very nicely done!
 Pauli
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
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 :)
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
But when are you doing MCODE analytic solutions to higher order polynomials? :)
 Pauli
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 HP15C handle complex numbers also, as you know:
http://www.hpmuseum.org/cgisys/cgiwrap/hpmuseum/archv017.cgi?read=114345
Regards,
Gerson.
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
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 Zsteps, 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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Thanks Paul!
I hope the 34s has a builtin QUAD command to do it in only one step :)
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
The thought had occurred to me to drop your code in as a built in function :)
 Pauli
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Accuracywise, Werner's code is better though. For instance, for a = 1E13, b = 2 and c = 1 it returns the correct 12digit 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.
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 1,253
Threads: 117
Joined: Nov 2005
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=C1 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=C1 MS b
0AE A<>C ALL
261 ?NC XQ
060 >1898 [DV2_10]
0E8 WRIT 3(X) register reuse…
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)^2c/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)^2c/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)^2c/a)
2BE C=C1 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
Posts: 4,587
Threads: 105
Joined: Jul 2005
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? ;)
Posts: 263
Threads: 42
Joined: Aug 2007
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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))^2c/a)'
'x2=b/(2*a)sqrt((b/(2*a))^2c/a)'
which is the quadratic formula for
'x^2b/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)^2c'
'x2=b/2sqrt((b/2)^2c'
Regards,
Gerson.
▼
Posts: 320
Threads: 59
Joined: Dec 2006
Gerson, I had done this many many years ago for, I think, my HP11c, 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+5x8=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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Chuck,
This appears quite suitable for hand calculations also, as the formula is easier to use and remember. My 14year 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.
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
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
Posts: 735
Threads: 34
Joined: May 2007
This program returns only one of the possible solutions of the equation x^{2} + px + q = 0. The quadratic solver hidden in ACOSH is used:
00 { 17Byte 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 HP15C as well. However it won't work with the HP35s due to its lack of proper support of complex functions.
Usage:
1) x^{2}  5x + 6 = 0
5 +/ ENTER 6 XEQ Q > X: 3
2) x^{2}  4x + 29 = 0
4 +/ ENTER 29 XEQ Q > X: 2 i5
3) x^{2}  (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
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
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 PR instead of ACOSH for the 41?
▼
Posts: 735
Threads: 34
Joined: May 2007
This could only be used when q <= 0. But even then you'd have to calculate the square root of q previously to use RP. So I'm afraid we won't gain anything. Only now I see you mentioned PR, not RP. What did you have in mind?
Edited: 23 Oct 2010, 11:02 a.m.
Posts: 735
Threads: 34
Joined: May 2007
Okay the following works only for certain values of p and q:
0 < 4q < p^{2}
00 { 20Byte 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^{2}
00 { 21Byte 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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hi Thomas,
Quote:
PS: Finding the other solution is left as an exercise.
Still unoptimized, as I'm in a hurry:
00 { 25Byte 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.
▼
Posts: 735
Threads: 34
Joined: May 2007
A little bit shorter:
00 { 21Byte 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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
You're right! x_{2} = q/x_{1}, for x_{1} <> 0, is shorter than x_{2} = (x_{1}+p).
This can be made shorter, I think.
Thanks for presenting your interesting insight.
Regards,
Gerson.
00 { 31Byte 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 x^{2}  5x + 6 = 0, compute 123456*x_{1} and 123456*x_{2}
123456 ENTER 1 ENTER 5 +/ ENTER 6 XEQ Q > Y: 2.99999999999 ; x_{2} = 3
X: 2.00000000001 ; x_{1} = 2
Rv * > X: 370,367.999999 ; 123456*x_{2}
Rv * > X: 246,912.000001 ; 123456*x_{1}
2) 3x^{2}  12x + 87 = 0
3 ENTER 12 +/ ENTER 87 XEQ Q > Y: 2.00E0 i5.00E0 ; x_{2} = 2 + 5i
X: 2.00000000001 i5 ; x_{1} = 2  5i
3) x^{2}  (4 + 6i)x + (5 + 10i) = 0
1 ENTER 4 ENTER 6 COMPLEX +/ 5 +/ ENTER 10 COMPLEX XEQ Q > Y: 3 i3.99999999999 ; x_{2} = 3 + 4i
X: 1.00E0 i2.00E0 ; x_{1} = 1 + 2i
It will work for a=1E13, b=2 and c=1 also:
1.99999999996E13
0.50000000001
▼
Posts: 735
Threads: 34
Joined: May 2007
It seems to be a recurring theme that we discuss every few years. Thanks for starting that thread (again).
Best regards
Thomas
