The solution to the quartic equation on the HP50g (long)

[Crawl]

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.

[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.

[Bruce Bergman]

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.

:-)

[Crawl]

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.

[C.Ret]

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]

[Crawl]

Thanks for the interest in my program.

After I noticed that error in the quartic method, I decided to upload the files to an online directory. So those should all be error free. I also uploaded C. Ret's Equa2, Equa3, as well as Equa4, which is my quartic code, but which calls his quadratic and cubic routines rather than my own.

You should be able to upload either the text files or binary (.hp) files to a calculator. The advantage of the text files is you can read the source code; the advantage of the .hp files is they can be loaded directly into EMU48.

Yes, the PICK commands are difficult to read, but there was a good reason for it. The goal was to be able to even solve the general quartic, and in reasonable time. You can be lackadaisical with the quadratic, and even the cubic to some degree. But the quartic is MUCH harder.

Here are some timed tests (using a stopwatch).

For my programs...
The general quartic: 3 minutes, 21 seconds. In my humble opinion, this is very fast for a pocket calculator.

The depressed quartic. That is, x^4+p.x^2+q.x+r=0.: 11 seconds. That's really fast! Plus, it can be viewed in pretty print through the equation editor.

The general cubic: Less than 2 seconds.

For C. Ret's programs.

The general cubic: About 23 seconds.

The depressed quartic with Equa4: Don't know. I let it run for over 15 minutes on a real calculator, in which time it didn't finish. I verified that it does work on EMU48, and pretty quickly (probably less than a minute).

The general quartic: Tried it on EMU48. Eventually got the message "INSUFFICIENT MEMORY"

I'm skeptical that the TI89's CAS can be used to solve the general quartic at all. This page has a quartic solver for the TI89, but it doesn't work sometimes even with numeric input (I think x^4-4.x^3+6.x^2-4.x+1=0 led to undefined). I ran it on VTI, which seems to be at least 10 times faster than a real TI89. For the general quartic, it ran for over an hour with no solution, and I believe it eventually caused the PC to crash.

Oh, but the TI89 has one thing that saved its honor during this challenge. I discovered another quirk with the HP50: When I tried making the substitution x=y+t into x^4+(b/a).x^3+...+(e/a), and then taking the taylor series, for some reason the HP50 just hung. It didn't finish in even 10 minutes. The TI89 does it almost instantly.

HOWEVER, if you have a=1, the HP50 also does it almost instantly. Maybe it gets caught up by worrying that a=0 or something?