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
(bsqrt(b^24ca))/(2a)
even though I entered it b^24ac, 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*CB^2)/(3*A^2)' ]
[ 'Q=(27*A^2*D9*(A*B*C)+2*B^3)/(27*A^3)' ]
[ 'Y=ZM/Z' ]
[ '(Z^6(3*MP)*Z^4+Q*Z^3+M*(3*MP)*Z^2M^3)/Z^3' ]
[ 'M=P/3' ]
[ 'Z^6+Q*Z^3P^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*C3*B^2)/(8*A^2)' ]
[ 'Q=(8*A^2*D4*A*B*C+B^3)/(8*A^3)' ]
[ 'R=(256*A^3*E64*A^2*B*D+16*A*B^2*C3*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*(KN)' ]
[ 'R=K*N' ]
[ 'K^2+N^2+2*K*N=(P+M^2)^2' ]
[ 'K^2+N^22*K*N=(Q/M)^2' ]
[ '4*K*N=(4*R=(P+M^2)^2(Q/M)^2)' ]
[ 'M^6+2*P*M^4+(P^24*R)*M^2Q^2=0' ]
[ 'N=(P+M^2Q/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*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)(3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))2*((8*A*C3*B^2)/(8*A^2))/3)\v/(\v/(XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)(3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))2*((8*A*C3*B^2)/(8*A^2))/3)^24*(((8*A*C3*B^2)/(8*A^2)+(XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)(3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))2*((8*A*C3*B^2)/(8*A^2))/3)(8*A^2*D4*A*B*C+B^3)/(8*A^3)/\v/(XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2)(3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3/(3*XROOT(3,((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)\v/((((27*((8*A^2*D4*A*B*C+B^3)/(8*A^3))^2)9*(2*((8*A*C3*B^2)/(8*A^2)))*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))+2*(2*((8*A*C3*B^2)/(8*A^2)))^3)/27)^2+4*(((3*(((8*A*C3*B^2)/(8*A^2))^24*((256*A^3*E64*A^2*B*D+16*A*B^2*C3*B^4)/(256*A^4)))(2*((8*A*C3*B^2)/(8*A^2)))^2)/3)^3/27)))/2))2*((8*A*C3*B^2)/(8*A^2))/3))/2)))/2B/(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.
