The Quartic Equation on the HP50g

(Note: This post is not actually finished yet.  I am posting it early so I can check it on my phone, and see if I can download all the programs from this page to my calculator)

I have a fondness for calculators, and particularly for the TI89 (/ Titanium) and HP50g (/ HP49g, HP49g+), which are capable of doing symbolic math such as algebra and even calculus.

After I got them, I tried putting them through their paces and seeing what they were and were not capable of.  Many people who take high school algebra learn to solve second-degree equations (the quadratic), and it turns out these calculators have the solution to these equations built-in.  It is also possible to solve third-degree (cubic) and fourth-degree (quartic) equations symbolically; the solutions to these equations have been known for hundreds of years.  (It has been proven that fifth-degree (quintic) equations do not have a general symbolic algebraic solution.) Many people don’t learn them, however, because the final solutions are complicated, tedious to derive or write down, and difficult to remember.  A computer algebra system will not complain about that, but solutions to those equations are not built in to those calculators.

In my opinion, the TI89 is not as well suited to dealing with those solutions, due to its tendency to try to simplify every expression. This creates problems when the expression is as complicated as the solution to the quartic. The HP50g is more passive, and (usually) only simplifies expressions when explicitly told to.

There were a few challenges.  One is that the HP50g does evaluate squareroots.  For example, if have

b^2 - 4ac

on the stack and press the squareroot key, you will not simply get the squareroot of that, but rather the calculator will re-arraign the terms.

A solution to this problem is that the calculator does not evaluate squareroots after a substitution.

The other trick was to use stack commands (eg., DUP, ROT) to prevent the calculator from evaluating pending expressions. This may cause the programs to be unreadable, but I consider that a necessary sacrifice.  RPL is not good for commenting code, anyway, and you shouldn’t really need to read the code; it can be downloaded from this page and then directly uploaded into a calculator.

Here are the raw (HP format) programs you will need.  They should all be given the same name they have here (minus the .hp extensions) (because the cubic program calls the quadratic, and the quartic calls the cubic), and uploaded to the same directory.

QUADR.hp
CUBIC.hp
QUARTIC.hp

In all cases, they assume 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.

These are the same programs in text format, so the source code is readable.  Some emulators require files to be in .hp format to be loaded, but the actual calculator should be able to load them in either format.

QUADR.txt
CUBIC.txt
QUARTIC.txt

The following files are intended to hint at the solution methods for the cubic and quartic.  They may be cryptic (and only understandable if you’ve seen the solutions before) but were the best concise solution I could come up with for this calculator.

CUBMETHD.hp
CUBMETHD.txt

QUARMETH.hp
QUARMETH.txt

Briefly, for the cubic, the substitution Y = X – B/(3A) is made to eliminate the second-order term.  Then the substitution Y = Z – M/Z is made.  M is solved to make the resulting equation a quadratic equation in Z^3.  (It turns out it is possible to eliminate both the Z^4 and Z^2 terms with the same value for M) This is solved using the normal quadratic equation solution.

For the quartic, similarly, Y = X – B/(4A) eliminates the third-order term.  The next goal is to factor the resulting quartic into the product of two quadratics.  That the quartic has no Y^3 term puts a simple condition on J and M. Now here’s a particularly clever, tricky step.  The resulting expressions for P and Q involve K-N and K+N.  They are squared, then combined to eliminate the K^2 and N^2 terms, leaving only the KN term.  But KN relates to R. This eliminates all the terms except M, which now appears in a cubic equation (in M^2).  This is solved using the previous cubic solution.  Then N is found with back substitution, giving one of the quadratics.  This is solved using the normal quadratic method.

What happens if you actually try to solve the general quartic on a physical HP50g? In about 3 minutes, you get this:

 

%%HP: T(3)A(D)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)’

 

I eventually decided I wanted a better solution.

This solution is basically just plugging into an already simplified solution. It “solves” the general case in about 10 seconds, and the result can be displayed in pretty print.

%%HP: T(3)A(R)F(.); \<< 0 0 0 0 \-> A B C D E ALPHA BETA GAMM DELTA \<< C 2 ^ 3 B D * * – 12 A E * * + ‘ALPHA’ STO 2 C 3 ^ * 9 B C D * * * – 27 A D 2 ^ * * + 27 B 2 ^ E * * + 72 A C E * * * – ‘BETA’ STO BETA ‘SQUAREROOT’ \v/ ‘SQUAREROOT’ -4 ALPHA 3 ^ * BETA 2 ^ + = SUBST + ‘GAMM’ STO ‘SQUAREROOT’ \v/ ‘SQUAREROOT’ B 2 ^ 4 A 2 ^ * / 2 C * 3 A * / – ALPHA 3 A GAMM 2 / 3 INV ^ * * / 1 3 A * / GAMM 2 / 3 INV ^ * + + = SUBST ‘DELTA’ STO B NEG 4 A * / 2 INV DELTA * – 2 INV NEG ‘SQUAREROOT’ \v/ ‘SQUAREROOT’ B 2 ^ 2 A 2 ^ * / 4 C * 3 A * / – ALPHA 3 A * GAMM 2 / 3 INV ^ * / – 3 A * INV GAMM 2 / 3 INV ^ * – B 3 ^ A 3 ^ / 4 B C * * A 2 ^ / – 8 D * A / + 4 DELTA * / + = SUBST * + \>> \>>

This is the file for downloading:

QUARTIC2.txt

The solution is:

%%HP: T(3)A(R)F(.); ‘((-(B/(4*A)))-((1/2)*(\v/((((B^2)/(4*(A^2)))-((2*C)/(3*A)))+(((((C^2)-(3*(B*D)))+(12*(A*E)))/(3*(A*XROOT(3,((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))+(\v/(((-4)*((((C^2)-(3*(B*D)))+(12*(A*E)))^3))+((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))^2))))/2))))+((1/(3*A))*XROOT(3,((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))+(\v/(((-4)*((((C^2)-(3*(B*D)))+(12*(A*E)))^3))+((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))^2))))/2)))))))-((1/2)*(\v/((((((B^2)/(2*(A^2)))-((4*C)/(3*A)))-((((C^2)-(3*(B*D)))+(12*(A*E)))/((3*A)*XROOT(3,((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))+(\v/(((-4)*((((C^2)-(3*(B*D)))+(12*(A*E)))^3))+((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))^2))))/2))))-((1/(3*A))*XROOT(3,((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))+(\v/(((-4)*((((C^2)-(3*(B*D)))+(12*(A*E)))^3))+((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))^2))))/2)))+(((((B^3)/(A^3))-((4*(B*C))/(A^2)))+((8*D)/A))/(4*(\v/((((B^2)/(4*(A^2)))-((2*C)/(3*A)))+(((((C^2)-(3*(B*D)))+(12*(A*E)))/(3*(A*XROOT(3,((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))+(\v/(((-4)*((((C^2)-(3*(B*D)))+(12*(A*E)))^3))+((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))^2))))/2))))+((1/(3*A))*XROOT(3,((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))+(\v/(((-4)*((((C^2)-(3*(B*D)))+(12*(A*E)))^3))+((((((2*(C^3))-(9*(B*(C*D))))+(27*(A*(D^2))))+(27*((B^2)*E)))-(72*(A*(C*E))))^2))))/2))))))))))’

Leave a Reply

Your email address will not be published. Required fields are marked *


four × 5 =