The Museum of HP Calculators

HP Forum Archive 19

 The solution to the quartic equation on the HP50g (long)Message #1 Posted by Crawl on 21 Oct 2009, 10:56 p.m. The HP50g has a nice CAS, but it cannot (automatically) symbolically solve cubic and quartic equations, though the solutions to those equations have been known for hundreds of years. I recently got a 1GB SD card just for my calculator, and my (unrealistic) goal is to eventually fill it with calculator stuff. So, this seemed like a good project. I learned a bit about the calculator on the way. In my first attempts, the quartic solution defeated the calculator -- it was just too long, it seemed. But I figured out what the problem was. The calculator can deal with an *expression* that big -- it just has a hard time *evaluating* it. It seemed like the solution was to use stack commands. After all, if you have X on rows 1 and 2, and hit *, it displays X*X -- it does not simplify to X^2. That should make things easier on the CPU. But I learned that *sometimes* the calculator does evaluate stack commands on symbolic input. I noticed this when I was just working on the simple quadratic equation. I noticed it returned (-b-sqrt(b^2-4ca))/(2a) even though I entered it b^2-4ac, with a and c in reverse order. So, every time you take a squareroot, the CAS tries to evaluate it. This was what was slowing the program down. I tried several approaches. What finally worked was when I realized the calculator does not automatically evaluate after a substitution. This is the solution to the quadratic equation. It is the shortest and easiest to see how I take squareroots. The other programs call it and assume it has name QUART3 (the 3 is just because I tried different versions before I got the squareroot to work) %%HP: T(3)A(R)F(.); \<< ROT DUP ROT * 4 SWAP * ROT DUP 2 ^ ROT - '\v/SQUAREROOT' SWAP 'SQUAREROOT' SWAP = SUBST SWAP NEG SWAP - SWAP 2 SWAP * / \>> I didn't want to just have the solution to the cubic and quartic equation -- I also wanted to show the method. But the programs don't really clearly demonstrate what the method is. So, I put the steps in row matrix form so they can be read. I was originally planning on having notes, but matrices can't have string entries. So, they might work more as a refresher if you once kind of knew how they were solved than a teacher. This is the method for solving cubics. ``` %%HP: T(3)A(R)F(.); [[ 'A*X^3+B*X^2+C*X+D' ] [ 'X=Y+T' ] [ 'Y^3+(3*A*T+B)/A*Y^2+((3*(A*T^2)+2*(B*T)+C)/A*Y+(A*T^3+B*T^2+C*T+D)/A)' ] [ 'T=-(B/(3*A))' ] [ 'Y^3+P*Y+Q' ] [ 'P=(3*A*C-B^2)/(3*A^2)' ] [ 'Q=(27*A^2*D-9*(A*B*C)+2*B^3)/(27*A^3)' ] [ 'Y=Z-M/Z' ] [ '(Z^6-(3*M-P)*Z^4+Q*Z^3+M*(3*M-P)*Z^2-M^3)/Z^3' ] [ 'M=P/3' ] [ 'Z^6+Q*Z^3-P^3/27' ]] ``` And this program solves cubics. The quartic program assumes it has the name CUBIC. %%HP: T(3)A(R)F(.); \<< 27 5 PICK 2 ^ * SWAP * 9 5 PICK * 4 PICK * 3 PICK * - 2 4 PICK 3 ^ * + 27 5 PICK 3 ^ * / 3 5 PICK * 3 PICK * 4 PICK 2 ^ - 3 6 PICK 2 ^ * / DUP UNROT 3 ^ 27 / NEG 1 UNROT QUAD3 3 INV ^ DUP UNROT 3 SWAP * / - 4 ROLLD DROP SWAP 3 SWAP * / - \>> This is the method for solving quartics. ``` %%HP: T(3)A(R)F(.); [[ 'A*X^4+B*X^3+C*X^2+D*X+E' ] [ 'X=Y+T' ] [ 'Y^4+(4*A*T+B)/A*Y^3+(6*A*T^2+3*B*T+C)/A*Y^2+(4*A*T^3+3*B*T^2+2*C*T+D)/A*Y+(A*T^4+B*T^3+C*T^2+D*T+E)/A' ] [ 'T=-(B/(4*A))' ] [ 'Y^4+P*Y^2+Q*Y+R' ] [ 'P=(8*A*C-3*B^2)/(8*A^2)' ] [ 'Q=(8*A^2*D-4*A*B*C+B^3)/(8*A^3)' ] [ 'R=(256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)' ] [ 'Y^4+P*Y^2+Q*Y+R=(Y^2+J*Y+K)*(Y^2+M*Y+N)' ] [ 'J+M=0' ] [ 'J=-M' ] [ 'P=KM^2' ] [ 'Q=M*(K-N)' ] [ 'R=K*N' ] [ 'K^2+N^2+2*K*N=(P+M^2)^2' ] [ 'K^2+N^2-2*K*N=(Q/M)^2' ] [ '4*K*N=(4*R=(P+M^2)^2-(Q/M)^2)' ] [ 'M^6+2*P*M^4+(P^2-4*R)*M^2-Q^2=0' ] [ 'N=(P+M^2-Q/M)/2' ]] ``` And this program does it %%HP: T(3)A(R)F(.); \<< 8 6 PICK * 4 PICK * 3 6 PICK 2 ^ * - 8 7 PICK 2 ^ * / 8 7 PICK 2 ^ * 4 PICK * 4 8 PICK * 7 PICK * 6 PICK * - 6 PICK 3 ^ + 8 8 PICK 3 ^ * / 256 8 PICK 3 ^ * 4 PICK * 64 9 PICK 2 ^ * 8 PICK * 6 PICK * - 16 9 PICK * 8 PICK 2 ^ * 7 PICK * + 3 8 PICK 4 ^ * - 256 9 PICK 4 ^ * / 2 4 PICK * 4 PICK 2 ^ 4 4 PICK * - 4 PICK 2 ^ NEG 1 4 ROLLD CUBIC DUP '\v/SQUAREROOT' SWAP 'SQUAREROOT' SWAP = SUBST 5 PICK 3 PICK + 5 PICK 3 PICK / - 2 / 1 UNROT QUAD3 10 ROLLD 7 DROPN SWAP 4 SWAP * / - \>> In all cases, it assumes the coefficient of the highest x power is entered first on the stack, and so on until the constant term is in level 1 of the stack. The programs work with both numeric or symbolic inputs. If you want a numeric result, however, you should be in approximate mode at the beginning. Otherwise, you might not be able to evaluate the result. The calculator has a quirk: It acts like it can't take cube roots of imaginary numbers. If you're in numeric mode the whole time, though, it will do it. I've noticed a few other quirks. Sometimes it doesn't simplify the cube root of zero. There are some cubic equation where you'd divide zero by the cube root of zero -- really, you should take a limit, and it would work. But the CAS evaluates zero divided by an unevaluated cube root of zero AS zero, even if a proper limit would give a different answer. In fact, this is specifically why I choose the sign of the squareroot in the quadratic solution the way I did. Otherwise, it leads to x^3+1=0 having solution ? (in exact mode) or 0.^(1/3) (in approx. mode). (The resolvent quadratic is x^2+x=0, and you want the -1 solution, not the 0 solution). There may be cases for which even the current solution might be numerically unstable, but I haven't checked everything. Anyway, what happens if you try to solve the general quartic symbolically? You get this %%HP: T(3)A(R)F(.); '(-\v/(XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)-(3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))-2*((8*A*C-3*B^2)/(8*A^2))/3)-\v/(\v/(XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)-(3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))-2*((8*A*C-3*B^2)/(8*A^2))/3)^2-4*(((8*A*C-3*B^2)/(8*A^2)+(XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)-(3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))-2*((8*A*C-3*B^2)/(8*A^2))/3)-(8*A^2*D-4*A*B*C+B^3)/(8*A^3)/\v/(XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)-(3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,(-((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)-\v/(((-(27*((8*A^2*D-4*A*B*C+B^3)/(8*A^3))^2)-9*(2*((8*A*C-3*B^2)/(8*A^2)))*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))+2*(2*((8*A*C-3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C-3*B^2)/(8*A^2))^2-4*((256*A^3*E-64*A^2*B*D+16*A*B^2*C-3*B^4)/(256*A^4)))-(2*((8*A*C-3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))-2*((8*A*C-3*B^2)/(8*A^2))/3))/2)))/2-B/(4*A)' It takes a couple of minutes to get that. I haven't timed it, but I believe it takes much less than 15 minutes (and on a real calculator, not an emulator). It seems the expression is too big to view in "pretty print". It displays fine as text, but when I select "graph", it just displays garbled nonsense. It's still an impressively large expression for the calculator to handle.

 Re: The solution to the quartic equation on the HP50g (long)Message #2 Posted by Crawl on 21 Oct 2009, 11:08 p.m.,in response to message #1 by Crawl I see there is at least one mistake there The quartic method matrix shows P=KM^2 when it should be P=K+N-M^2 It is correct on my calculator; an HTML formatting bug must have messed it up on the post. I'll have to check later if there are other errors.

 Re: The solution to the quartic equation on the HP50g (long)Message #3 Posted by Bruce Bergman on 22 Oct 2009, 12:28 a.m.,in response to message #1 by Crawl Goodness, how EVER did you type that in? Who cares about the equation! I just want to know how long it took you to type that. :-)

 Re: The solution to the quartic equation on the HP50g (long)Message #4 Posted by Crawl on 22 Oct 2009, 6:08 a.m.,in response to message #3 by Bruce Bergman That's why I needed the calculator to do the solution -- I didn't want to type it in. After the calculator got the solution, I saved it as a variable, then transferred it to PC with the USB cable, then just copy pasted. So, it actually was fairly quick and easy.

 Re: The solution to the quartic equation on the HP50g (long)Message #5 Posted by C.Ret on 23 Oct 2009, 3:58 a.m.,in response to message #3 by Bruce Bergman Quote:Who cares about the equation! It is always a good idea to give a minimum of documentation. Any reader of this forum may be interesting in such documentation and any educated people may have recognized the Cardan's method of resolution. Fortunately, this is also the method I used to solve 3rd degree equation. Despite Crawl's code, which is using a lot of obscure PICK command, my solution use local variables, which make it a bit more readable. But still, in my production the exhaustive list of intermediate equations is missing. As Crawl's code, the EQU3 function takes equation the coefficients from the stack (the coefficient of the highest x power is entered first on the stack) and uses the EQU2 to solve the intermediate quadratic equation. This code is originally for HP28C/S, assuming that LAST command is enable. But, it may be easily ported to any HP48/49 or 50, as it only simples RPL commands: ```« -> a b c @ a.x2+b.x+c=0 where a<>0 « b NEG 2 / a / @ term1: -b/2a b SQ 4 a c * * - SQRT 2 / a / @ term2: SQRT(4ac-b^2)/2a - @ x2 : term1-term2 LAST @ recover term1 and term2 + @ x1 : term1+term2 » ‘EQUA2’ STO Usage 3: a ---- EQUA2 ----> 2: x2 2: b 1: x1 1: c ``` a,b,c coefficients of the quadratic equation (real, complex or symbolic) x1,x2 two roots of the equation (real, complex or symbolic). Example: 4 1 CHS 1 CHS EQUA2 returns 2: –0.3904 and 1: 0.6404 (assuming 4 FIX display mode) 1 –2 1 EQUA2 returns 2: 1.0000 1: 1.0000 (double roots are simply repeated) ```« -> a b c d @ a.x3 + b.x2 + c.x +d = 0 « c a / 3 / b SQ a SQ / 9 / - @ p : B/3–A^2/27 where B=c/a and A=b/a d a / b c * a / 3 / - 2 b 3 ^ * a 3 ^ / 27 / + @ q : C-A.B/3+2A^3/27 where C=d/a -> p q « @ quadratic equation t^2+q.t-p^3=0 1 @ first coefficient (t^2) q @ second coefficient (t^1) p 3 ^ NEG @ third coefficient (t^0) EQUA2 @ compute the two roots u3 and v3 3 INV ^ @ extract v from v3 SWAP 3 INV ^ @ extract u from u3 b a / 3 @ dx : B/3 for future usage (translation) -> v u dx « u v + @ DUP dx - @ x1 : u+v – dx SWAP –2 / @ term1: -(u+v)/2 -3 SQRT 2 /u v - * @ term2: i.SQRT(3).(u-v)/2 + @ x2+dx: term1+term2 LAST - @ x3+dx: term1-term2 dx – SWAP dx - » » » » ‘EQUA3’ STO Usage 4: a 3: b ---- EQUA3 ----> 3: x1 2: c 2: x2 1: d 1: x3 ``` Exemple 1 –3 3 –1 EQUA3 returns 1.0000 (1.0000,0.0000) (1.0000,0.0000) Most of the time x1 is expressed as a real,and, x2 & x3 are express as complex. [/code]