▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-21-2010, 10:20 AM 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 ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-21-2010, 02:49 PM 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 ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-21-2010, 03:24 PM 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 ``` ▼ Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 10-22-2010, 08:53 AM 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 ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 11:50 AM 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. ▼ Csaba Tizedes (Hungary) Member Posts: 59 Threads: 9 Joined: Apr 2008 10-24-2010, 06:08 PM 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! ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-25-2010, 03:23 PM 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. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-23-2010, 12:49 PM 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. Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 10-26-2010, 04:05 AM 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 ``` Eddie W. Shore Posting Freak Posts: 764 Threads: 118 Joined: Aug 2007 10-22-2010, 09:12 AM Wow, nice job! Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-22-2010, 03:44 AM Very nicely done! - Pauli ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 10-22-2010, 06:16 AM 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 :) ▼ Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-22-2010, 06:33 AM But when are you doing MCODE analytic solutions to higher order polynomials? :-) - Pauli ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 10-22-2010, 08:16 AM Touché :) Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 07:47 AM 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: Regards, Gerson. ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 10-22-2010, 08:23 AM 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. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 06:14 PM Thanks Paul! I hope the 34s has a built-in QUAD command to do it in only one step :-) ▼ Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-22-2010, 06:47 PM The thought had occurred to me to drop your code in as a built in function :-) - Pauli ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 07:17 PM 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. ▼ Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-22-2010, 08:00 PM 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 W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 09:06 PM 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. Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 10-22-2010, 09:12 PM 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 ``` Walter B Posting Freak Posts: 4,587 Threads: 105 Joined: Jul 2005 10-23-2010, 05:33 AM 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? ;) Thomas Chrapkiewicz Senior Member Posts: 263 Threads: 42 Joined: Aug 2007 10-22-2010, 09:59 AM 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 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 06:01 PM 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. ▼ Chuck Senior Member Posts: 320 Threads: 59 Joined: Dec 2006 10-22-2010, 08:29 PM 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 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-22-2010, 10:00 PM 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. ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 10-23-2010, 02:29 AM 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 Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-23-2010, 02:11 AM 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 ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 10-23-2010, 02:32 AM 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? ▼ Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-23-2010, 03:00 AM 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. Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-23-2010, 04:00 AM 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. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-23-2010, 06:21 AM 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. ▼ Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-23-2010, 07:41 AM 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. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-23-2010, 09:21 AM 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 ``` ▼ Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-23-2010, 11:50 AM It seems to be a recurring theme that we discuss every few years. Thanks for starting that thread (again). Best regards Thomas