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.