HP35S - Cubic (and incidentally quadratic) equations solver

[Roberto Volpi]

Hi all,

inspired by similar programs from classic HP models, here a short program to solve cubic equations in the form of Ax^3+Bx^2+Cx+D=0

The algorithm search the 1st root x1, always real, divides the polynomial by (x-x1) with synthetic division (Ruffini rule) to obtain a quadratic equation, and solve it to have 2nd and 3rd roots x2 and x3.

In case we have a quadratic equation in the form Ax^2+Bx+C=0 to solve, it processes directly to obtain the roots, real or complex.

LINE INSTRUCTION
C001 LBL C
C002 CLVARS
C003 STO D
C004 R down
C005 STO C
C006 R down
C007 STO B
C008 R down
C009 STO/B
C010 STO/C
C011 STO/D
C012 RCL D
C013 SGN
C014 STO H
C015 LAST X
C016 ABS
C017 RCL B
C018 ABS
C019 +
C020 10
C021 STO P
C022 R down
C023 1E-6
C024 STO O
C025 RCLxP
C026 X<=Y?
C027 GTO C025
C028 STO K
C029 RCL P
C030 STO/K
C031 RCL H
C032 +/-
C033 STO H
C034 RCL K
C035 RCLxH
C036 RCL+F
C037 STO G
C038 RCL F
C039 X=Y?
C040 GTO C059
C041 RCL G
C042 RCL+B
C043 RCLxG
C044 RCL+C
C045 RCLxG
C046 RCL+D
C047 STO L
C048 ABS
C049 RCL O
C050 X>Y?
C051 GTO C059
C052 RCL G
C053 STO F
C054 RCL L
C055 RCLxH
C056 X<0?
C057 GTO C034
C058 GTO C029
C059 RCL G
C060 R/S
C061 1
C062 X<>Y
C063 RCL+B
C064 ENTER
C065 RCLxG
C066 RCL+C
C067 GTO E001
E001 LBL E
E002 -2 i 0
E003 x
E004 LAST X
E005 R up
E006 x
E007 /
E008 X<>Y
E009 LAST X
E010 /
E011 ENTER
E012 ENTER
E013 x
E014 LAST X
E015 X<>Y
E016 R up
E017 -
E018 0.5
E019 Y^X
E020 -
E021 X<>Y
E022 LAST X
E023 +
E024 RTN

Instructions to solve cubic equation:

input A
press ENTER
input B
press ENTER
input C
press ENTER
input D
press XEQ C ENTER

the calculator shows first root x1, always real, in stack X, and the accuracy in stack Y

press R/S

the calculator shows 2nd and 3rd roots, x2 and x3, in stack Y and X.

The first real root x1 can be retrieved by pressing RCL G

Instructions to solve quadratic equations

input A
press ENTER
input B
press ENTER
input C
press XEQ E ENTER

the calculator shows both roots in stacks X and Y

Enjoy!

[PedroLeiva]

Nice program Roberto, thanks for sharing. Just to see if I understood the instructions and didn't make typing mistakes, ¿could you give me two numerical examples, one for cubic and the other for quadratic?. Pedro

[Roberto Volpi]

(04-07-2023 07:55 PM)PedroLeiva Wrote:  Nice program Roberto, thanks for sharing. Just to see if I understood the instructions and didn't make typing mistakes, ¿could you give me two numerical examples, one for cubic and the other for quadratic?. Pedro

Thanks Pedro

Let's take as example the following quadratic EQN:

x^2 -x -12

Here how to calculate both roots:

Press 1
Press ENTER
Press -1
Press ENTER
Press -12
Press XEQ E ENTER

The calculator will give following result:

Stack Y: -3 i 0 (1st root)
Stack X: 4 i 0 (2nd root)

-----

Now a quick example of cubic EQN:

x^3 -13x -12

Here how to calculate all 3 roots:

Press 1
Press ENTER
Press 0
Press ENTER
Press -13
Press -12
Press XEQ C ENTER

The calculator will give as result:

Stack Y: 0.000001 (theoretical approximation)
Stack X: 4 (1st root, always real)

Press R/S

The calculator will give as result:

Stack Y: -3 i 0 (2nd root)
Stack X: -1 i 0 (3rd root)

in case you forget to take note of 1st root, just press RCL G.

[PedroLeiva]

Excellent Roberto, the program works very well. Thanks a lot. Pedro

[Thomas Klemm]

As usual we divide all coefficients of the original polynomial \(ax^3 + bx^2 + cx + d\) by \(a\), so we start with the following polynomial:

\(
f(x) = x^3 + px^2 + qx + r
\)

Using Horner's method we can reuse the intermediate values \(x_0 + p \) and \(x_0(x_0 + p) + q\) while calculating the derivative:

\(
\begin{matrix}
& 1 & p & q & r \\
x_0 & 1 & x_0 + p & x_0(x_0 + p) + q & x_0(x_0(x_0 + p) + q) + r = f(x_0) \\
x_0 & 1 & 2 x_0 + p & x_0(3 x_0 + 2 p) + q = f'(x_0) & \\
\end{matrix}

\)

This allows calculating the derivative and synthetic division in one step.

Registers

R00: \(p = \frac{b}{a}\)
R01: \(q = \frac{c}{a}\)
R02: \(r = \frac{d}{a}\)
R03: \(x_0 + p\)
R04: \(x_0(x_0 + p) + q\)

HP-42S

Code:

00 { 69-Byte Prgm }
01 LBL "CBRT"
02 R↑
03 STO÷ ST T
04 STO÷ ST Z
05 ÷
06 STO 02
07 R↑
08 STO 00
09 X<> ST Z
10 STO 01
11 ÷
12 +/-
13 LBL 00
14 ENTER
15 ENTER
16 ENTER
17 RCL+ 00
18 STO 03
19 ×
20 RCL+ 01
21 STO 04
22 ×
23 RCL+ 02
24 R↑
25 ENTER
26 RCL+ 03
27 ×
28 RCL+ 04
29 ÷
30 -
31 X≠Y?
32 GTO 00
33 STOP
34 RCL 03
35 -2
36 ÷
37 ENTER
38 ENTER
39 X↑2
40 RCL- 04
41 SQRT
42 STO- ST Z
43 +
44 END

HP-15C

Code:

   001 {    42 21 13 } f LBL C
   002 {       43 33 } g R⬆
   003 {          10 } ÷
   004 {       44  2 } STO 2
   005 {          33 } R⬇
   006 {       43 36 } g LSTΧ
   007 {          10 } ÷
   008 {       44  1 } STO 1
   009 {          34 } x↔y
   010 {       43 36 } g LSTΧ
   011 {          10 } ÷
   012 {       44  0 } STO 0
   013 {          33 } R⬇
   014 {          10 } ÷
   015 {          16 } CHS
   016 {    42 21  0 } f LBL 0
   017 {          36 } ENTER
   018 {          36 } ENTER
   019 {          36 } ENTER
   020 {    45 40  0 } RCL + 0
   021 {       44  3 } STO 3
   022 {          20 } ×
   023 {    45 40  1 } RCL + 1
   024 {       44  4 } STO 4
   025 {          20 } ×
   026 {    45 40  2 } RCL + 2
   027 {       43 33 } g R⬆
   028 {    45 40  3 } RCL + 3
   029 {       43 33 } g R⬆
   030 {          20 } ×
   031 {    45 40  4 } RCL + 4
   032 {          10 } ÷
   033 {          30 } −
   034 {    43 30  6 } g TEST x≠y
   035 {       22  0 } GTO 0
   036 {          31 } R/S
   037 {       45  3 } RCL 3
   038 {           2 } 2
   039 {          16 } CHS
   040 {          10 } ÷
   041 {          36 } ENTER
   042 {          36 } ENTER
   043 {       43 11 } g x²
   044 {    45 30  4 } RCL − 4
   045 {          11 } √x̅
   046 {          30 } −
   047 {          34 } x↔y
   048 {       43 36 } g LSTΧ
   049 {          40 } +
   050 {       43 32 } g RTN

Example

\(
x^3 - x^2 - x - 1 = 0
\)

1
ENTER
-1
ENTER
ENTER
XEQ "CBRT"

X: 1.83928675521

R/S

Y: -0.41964 - i0.60629
X: -0.41964 + i0.60629

Caveat

The value \(- \frac{r}{q}\) is used as initial value for Newton's iteration.
This is the next guess when you start with \(0\).
However, this might not always be the best value.

I was a bit sloppy with the termination criterion.
It may never be reached, but I've never experienced it.

[Roberto Volpi]

(04-10-2023 02:21 PM)Thomas Klemm Wrote:  As usual we divide all coefficients of the original polynomial \(ax^3 + bx^2 + cx + d\) by \(a\), so we start with the following polynomial:

\(
f(x) = x^3 + px^2 + qx + r
\)

Using Horner's method we can reuse the intermediate values \(x_0 + p \) and \(x_0(x_0 + p) + q\) while calculating the derivative:

\(
\begin{matrix}
& 1 & p & q & r \\
x_0 & 1 & x_0 + p & x_0(x_0 + p) + q & x_0(x_0(x_0 + p) + q) + r = f(x_0) \\
x_0 & 1 & 2 x_0 + p & x_0(3 x_0 + 2 p) + q = f'(x_0) & \\
\end{matrix}

\)

This allows calculating the derivative and synthetic division in one step.

Registers

R00: \(p = \frac{b}{a}\)
R01: \(q = \frac{c}{a}\)
R02: \(r = \frac{d}{a}\)
R03: \(x_0 + p\)
R04: \(x_0(x_0 + p) + q\)

HP-42S

Code:

00 { 69-Byte Prgm }
01 LBL "CBRT"
02 R↑
03 STO÷ ST T
04 STO÷ ST Z
05 ÷
06 STO 02
07 R↑
08 STO 00
09 X<> ST Z
10 STO 01
11 ÷
12 +/-
13 LBL 00
14 ENTER
15 ENTER
16 ENTER
17 RCL+ 00
18 STO 03
19 ×
20 RCL+ 01
21 STO 04
22 ×
23 RCL+ 02
24 R↑
25 ENTER
26 RCL+ 03
27 ×
28 RCL+ 04
29 ÷
30 -
31 X≠Y?
32 GTO 00
33 STOP
34 RCL 03
35 -2
36 ÷
37 ENTER
38 ENTER
39 X↑2
40 RCL- 04
41 SQRT
42 STO- ST Z
43 +
44 END

HP-15C

Code:

   001 {    42 21 13 } f LBL C
   002 {       43 33 } g R⬆
   003 {          10 } ÷
   004 {       44  2 } STO 2
   005 {          33 } R⬇
   006 {       43 36 } g LSTΧ
   007 {          10 } ÷
   008 {       44  1 } STO 1
   009 {          34 } x↔y
   010 {       43 36 } g LSTΧ
   011 {          10 } ÷
   012 {       44  0 } STO 0
   013 {          33 } R⬇
   014 {          10 } ÷
   015 {          16 } CHS
   016 {    42 21  0 } f LBL 0
   017 {          36 } ENTER
   018 {          36 } ENTER
   019 {          36 } ENTER
   020 {    45 40  0 } RCL + 0
   021 {       44  3 } STO 3
   022 {          20 } ×
   023 {    45 40  1 } RCL + 1
   024 {       44  4 } STO 4
   025 {          20 } ×
   026 {    45 40  2 } RCL + 2
   027 {       43 33 } g R⬆
   028 {    45 40  3 } RCL + 3
   029 {       43 33 } g R⬆
   030 {          20 } ×
   031 {    45 40  4 } RCL + 4
   032 {          10 } ÷
   033 {          30 } −
   034 {    43 30  6 } g TEST x≠y
   035 {       22  0 } GTO 0
   036 {          31 } R/S
   037 {       45  3 } RCL 3
   038 {           2 } 2
   039 {          16 } CHS
   040 {          10 } ÷
   041 {          36 } ENTER
   042 {          36 } ENTER
   043 {       43 11 } g x²
   044 {    45 30  4 } RCL − 4
   045 {          11 } √x̅
   046 {          30 } −
   047 {          34 } x↔y
   048 {       43 36 } g LSTΧ
   049 {          40 } +
   050 {       43 32 } g RTN

Example

\(
x^3 - x^2 - x - 1 = 0
\)

1
ENTER
-1
ENTER
ENTER
XEQ "CBRT"

X: 1.83928675521

R/S

Y: -0.41964 - i0.60629
X: -0.41964 + i0.60629

Caveat

The value \(- \frac{r}{q}\) is used as initial value for Newton's iteration.
This is the next guess when you start with \(0\).
However, this might not always be the best value.

I was a bit sloppy with the termination criterion.
It may never be reached, but I've never experienced it.

Just tried the version for HP15C, but after obtaining the 1st root (20 seconds in an iPhone emulator, which is much much faster than a real HP15C), R/S just gives ERROR 0.

Can you confirm specific instructions also for HP15C?

Heartfelt thanks!

[Thomas Klemm]

(04-20-2023 07:18 AM)Roberto Volpi Wrote:  R/S just gives ERROR 0.

From HP-15C Owner's Handbook p. 121:

Quote:Complex mode is activated

1. automatically, when executing f I or f Re<>Im; or
2. by setting flag 8, the Complex mode flag.
   

Either set flag 8 before running the program or set it only if \(X < 0\) before calculating the square root in the following line:

Code:

   045 {          11 } √x̅

[Roberto Volpi]

(04-20-2023 11:36 AM)Thomas Klemm Wrote:  

(04-20-2023 07:18 AM)Roberto Volpi Wrote:  R/S just gives ERROR 0.

From HP-15C Owner's Handbook p. 121:

Quote:Complex mode is activated

1. automatically, when executing f I or f Re<>Im; or
2. by setting flag 8, the Complex mode flag.
   

Thanks for your prompt answer.

Either set flag 8 before running the program or set it only if \(X < 0\) before calculating the square root in the following line:

Code:

   045 {          11 } √x̅

Tried again with complex mode activated, and it is ok, except when a cubic polynomial with multiple roots is entered.

Ex. x^3+3x^2+3x+1, instead of giving -1 gives ERROR 0

The same for x^3-3x^3+3x-1.

[Thomas Klemm]

(05-12-2023 07:49 AM)Roberto Volpi Wrote:  … it is ok, except when a cubic polynomial with multiple roots is entered.

Thank you for pointing that out.
In these cases the derivative at this root is zero leading to the \(y \div 0\) error in line 32.
You can wiggle the coefficients a bit to get roots that are close to the exact roots.

Example

1.000000001
ENTER
3
ENTER
ENTER
1

f C
-0.999031396

R/S
-1.000484301

f (i)
0.000839643

x<>y
-1.000484301

f (i)
0.000839643

As you can see they might turn out to be complex though.

---

How do you handle this special case in your program?

[gor1060]

   

[Gerson W. Barbosa]

(04-10-2023 02:21 PM)Thomas Klemm Wrote:  Caveat

I was a bit sloppy with the termination criterion.
It may never be reached, but I've never experienced it.

I have, on Free42 with the first example I tried: 1 ENTER 2 ENTER 3 ENTER 4 XEQ “CBRT” (No problem on the real 42S, at least with this example).
The insertion of RND between steps 30 and 31 might be a workaround, but accuracy becomes limited to only twelve significant digits on Free42. Also, this will cause your example to run forever with FIX 1.

Gerson.

[Roberto Volpi]

---

How do you handle this special case in your program?
[/quote]

No special care is needed.

Just input the cubic equation, and in case it has repeated roots they will be shown on the stacks.

EX. X^3+3X^2+3X+1

1st root always real: -1

2nd and 3rd root:

stack y: -1
stack x: -1

Just it.

[Thomas Klemm]

(04-07-2023 09:31 AM)Roberto Volpi Wrote:  C026 X>=Y?

Could it be that this should rather be:

C026 X<=Y?

[Roberto Volpi]

(05-20-2023 09:45 AM)Thomas Klemm Wrote:  

(04-07-2023 09:31 AM)Roberto Volpi Wrote:  C026 X>=Y?

Could it be that this should rather be:

C026 X<=Y?

Thanks Thomas,

Just chequed the instruction C026 and you were right! It was a bug.

Correction has been made in the original text.

Thanks once more for your valuable contribution.

[Thomas Klemm]

It took me a while to understand how your program works.
The following Python program helped me to do so:

Code:

from math import copysign

def cbeq(a, b, c, d):
    B, C, D = b / a, c / a, d / a  # C002 - C011
    F = 0
    H = copysign(1, D)  # C012 - C014
    X = abs(D) + abs(B)  # C015 - C019
    P = 10  # C020 - C022
    O = 1e-6  # C023 - C024
    Y, X = X, O
    X *= P  # C025
    print(f"C025: {X=:15.12f} {Y=:15.12f} {X <= Y}")
    while X <= Y:  # C026 - C027
        X *= P  # C025
        print(f"C025: {X=:15.12f} {Y=:15.12f} {X <= Y}")
    K = X  # C028
    L = 0
    while True:
        if L * H >= 0:  # C054 - C058
            K /= P  # C029 - C030
            H = -H  # C031 - C033
        G = K * H + F  # C034 - C037
        if G == F:  # C038 - C039
            return G  # C040
        L = ((G + B) * G + C) * G + D  # C041 - C047
        print(f"C048: {H=:2.0f} {K=:15.12f} {G=:15.12f} {L=:15.12f}")
        if O > abs(L):  # C048 - C050
            return G  # C051
        F = G  # C052 - C053

First it tries to find an integral power of 10 that is bigger or equal to \(|\frac{B}{A}| + |\frac{D}{A}|\).
Then a digit by digit algorithm is used to approach the solution.
Whenever a change of sign is found the next decimal place to the right is considered.

Here's an example for the equation \(x^3 - x^2 - x - 1 = 0\):

cbeq(1, -1, -1, -1)

Code:

C025: X= 0.000010000000 Y= 2.000000000000 True
C025: X= 0.000100000000 Y= 2.000000000000 True
C025: X= 0.001000000000 Y= 2.000000000000 True
C025: X= 0.010000000000 Y= 2.000000000000 True
C025: X= 0.100000000000 Y= 2.000000000000 True
C025: X= 1.000000000000 Y= 2.000000000000 True
C025: X=10.000000000000 Y= 2.000000000000 False
C048: H= 1 K= 1.000000000000 G= 1.000000000000 L=-2.000000000000
C048: H= 1 K= 1.000000000000 G= 2.000000000000 L= 1.000000000000
C048: H=-1 K= 0.100000000000 G= 1.900000000000 L= 0.349000000000
C048: H=-1 K= 0.100000000000 G= 1.800000000000 L=-0.208000000000
C048: H= 1 K= 0.010000000000 G= 1.810000000000 L=-0.156359000000
C048: H= 1 K= 0.010000000000 G= 1.820000000000 L=-0.103832000000
C048: H= 1 K= 0.010000000000 G= 1.830000000000 L=-0.050413000000
C048: H= 1 K= 0.010000000000 G= 1.840000000000 L= 0.003904000000
C048: H=-1 K= 0.001000000000 G= 1.839000000000 L=-0.001568281000
C048: H= 1 K= 0.000100000000 G= 1.839100000000 L=-0.001021459529
C048: H= 1 K= 0.000100000000 G= 1.839200000000 L=-0.000474547712
C048: H= 1 K= 0.000100000000 G= 1.839300000000 L= 0.000072454457
C048: H=-1 K= 0.000010000000 G= 1.839290000000 L= 0.000017750174
C048: H=-1 K= 0.000010000000 G= 1.839280000000 L=-0.000036953205
C048: H= 1 K= 0.000001000000 G= 1.839281000000 L=-0.000031482908
C048: H= 1 K= 0.000001000000 G= 1.839282000000 L=-0.000026012602
C048: H= 1 K= 0.000001000000 G= 1.839283000000 L=-0.000020542286
C048: H= 1 K= 0.000001000000 G= 1.839284000000 L=-0.000015071962
C048: H= 1 K= 0.000001000000 G= 1.839285000000 L=-0.000009601629
C048: H= 1 K= 0.000001000000 G= 1.839286000000 L=-0.000004131286
C048: H= 1 K= 0.000001000000 G= 1.839287000000 L= 0.000001339065
C048: H=-1 K= 0.000000100000 G= 1.839286900000 L= 0.000000792030

And we end up with:

G= 1.839286900000


But this made me think, what happens if the solution gets closer and closer to 1e-6.

cbeq(1e3, -1e2, -1e1, -1)

Code:

C025: X= 0.000010000000 Y= 0.101000000000 True
C025: X= 0.000100000000 Y= 0.101000000000 True
C025: X= 0.001000000000 Y= 0.101000000000 True
C025: X= 0.010000000000 Y= 0.101000000000 True
C025: X= 0.100000000000 Y= 0.101000000000 True
C025: X= 1.000000000000 Y= 0.101000000000 False
C048: H= 1 K= 0.100000000000 G= 0.100000000000 L=-0.002000000000
C048: H= 1 K= 0.100000000000 G= 0.200000000000 L= 0.001000000000
C048: H=-1 K= 0.010000000000 G= 0.190000000000 L= 0.000349000000
C048: H=-1 K= 0.010000000000 G= 0.180000000000 L=-0.000208000000
C048: H= 1 K= 0.001000000000 G= 0.181000000000 L=-0.000156359000
C048: H= 1 K= 0.001000000000 G= 0.182000000000 L=-0.000103832000
C048: H= 1 K= 0.001000000000 G= 0.183000000000 L=-0.000050413000
C048: H= 1 K= 0.001000000000 G= 0.184000000000 L= 0.000003904000
C048: H=-1 K= 0.000100000000 G= 0.183900000000 L=-0.000001568281
C048: H= 1 K= 0.000010000000 G= 0.183910000000 L=-0.000001021460
C048: H= 1 K= 0.000010000000 G= 0.183920000000 L=-0.000000474548

G= 0.183920000000

---

cbeq(1e6, -1e4, -1e2, -1)

Code:

C025: X= 0.000010000000 Y= 0.010001000000 True
C025: X= 0.000100000000 Y= 0.010001000000 True
C025: X= 0.001000000000 Y= 0.010001000000 True
C025: X= 0.010000000000 Y= 0.010001000000 True
C025: X= 0.100000000000 Y= 0.010001000000 False
C048: H= 1 K= 0.010000000000 G= 0.010000000000 L=-0.000002000000
C048: H= 1 K= 0.010000000000 G= 0.020000000000 L= 0.000001000000
C048: H=-1 K= 0.001000000000 G= 0.019000000000 L= 0.000000349000

G= 0.019000000000

---

cbeq(1e9, -1e6, -1e3, -1)

Code:

C025: X= 0.000010000000 Y= 0.001000001000 True
C025: X= 0.000100000000 Y= 0.001000001000 True
C025: X= 0.001000000000 Y= 0.001000001000 True
C025: X= 0.010000000000 Y= 0.001000001000 False
C048: H= 1 K= 0.001000000000 G= 0.001000000000 L=-0.000000002000

G= 0.001000000000


Unfortunately it looks like fewer and fewer digits are calculated.

[Thomas Klemm]

We can use Horner's method to calculate the Taylor polynomial at \(x=1\) of the initial cubic polynomial:

\(
ax^3 + bx^2 + cx + d = a (x - 1)^3 + (3 a + b) (x - 1)^2 + (3 a + 2 b + c) (x - 1) + (a + b + c + d)
\)

This leads to the following outline:

\(
\begin{matrix}
a & b & c & d \\
a & a + b & a + b + c & a + b + c + d \\
a & 2a + b & 3a + 2b + c & \\
a & 3a + b & & \\
a & & & \\
\end{matrix}
\)

The following Python program uses the shift function to do this:

Code:

def shift(x, p):
    [a, b, c, d] = p
    b += a
    c += b
    d += c
    b += a
    c += b
    b += a
    return x+1, [a, b, c, d]

def zoom(x, p):
    [a, b, c, d] = p
    return 10*x, [a, 10*b, 100*c, 1000*d]

def solve(p, n):
    x = 0
    for _ in range(n):
        while p[3] <= 0:
            y, q = x, p
            x, p = shift(x, p)
        x, p = zoom(y, q)
    return x

Once a change of sign is found the zoom function is used with the previous values.
From then on the next decimal place to the right is considered.
This allows calculating the root using only addition and multiplication by ten which is a simple shift in BCD.

Examples:

solve([1, -1, -1, -1], 75)

1839286755214161132551852564653286600424178746097592246778758639404203222080

\(x \approx 1.83928675521416113255185256465328660042417874609759224677875863940420322208\)

solve([1, -2, 3, -4], 75)

1650629191439388218880800967426197435895495231706172591999594185884157996140

\(x \approx 1.65062919143938821888080096742619743589549523170617259199959418588415799614\)

solve([1, -3, 3, -1], 10)

10000000000

\(x=1\)

---

Please consider this just a sketch of an algorithm.
But those familiar with computing square roots in HP calculators may recognize the similarities.

[Roberto Volpi]

(05-20-2023 12:05 PM)Thomas Klemm Wrote:  It took me a while to understand how your program works.
The following Python program helped me to do so:

Code:

from math import copysign

def cbeq(a, b, c, d):
    B, C, D = b / a, c / a, d / a  # C002 - C011
    F = 0
    H = copysign(1, D)  # C012 - C014
    X = abs(D) + abs(B)  # C015 - C019
    P = 10  # C020 - C022
    O = 1e-6  # C023 - C024
    Y, X = X, O
    X *= P  # C025
    print(f"C025: {X=:15.12f} {Y=:15.12f} {X <= Y}")
    while X <= Y:  # C026 - C027
        X *= P  # C025
        print(f"C025: {X=:15.12f} {Y=:15.12f} {X <= Y}")
    K = X  # C028
    L = 0
    while True:
        if L * H >= 0:  # C054 - C058
            K /= P  # C029 - C030
            H = -H  # C031 - C033
        G = K * H + F  # C034 - C037
        if G == F:  # C038 - C039
            return G  # C040
        L = ((G + B) * G + C) * G + D  # C041 - C047
        print(f"C048: {H=:2.0f} {K=:15.12f} {G=:15.12f} {L=:15.12f}")
        if O > abs(L):  # C048 - C050
            return G  # C051
        F = G  # C052 - C053

First it tries to find an integral power of 10 that is bigger or equal to \(|\frac{B}{A}| + |\frac{D}{A}|\).
Then a digit by digit algorithm is used to approach the solution.
Whenever a change of sign is found the next decimal place to the right is considered.

Here's an example for the equation \(x^3 - x^2 - x - 1 = 0\):

cbeq(1, -1, -1, -1)

Code:

C025: X= 0.000010000000 Y= 2.000000000000 True
C025: X= 0.000100000000 Y= 2.000000000000 True
C025: X= 0.001000000000 Y= 2.000000000000 True
C025: X= 0.010000000000 Y= 2.000000000000 True
C025: X= 0.100000000000 Y= 2.000000000000 True
C025: X= 1.000000000000 Y= 2.000000000000 True
C025: X=10.000000000000 Y= 2.000000000000 False
C048: H= 1 K= 1.000000000000 G= 1.000000000000 L=-2.000000000000
C048: H= 1 K= 1.000000000000 G= 2.000000000000 L= 1.000000000000
C048: H=-1 K= 0.100000000000 G= 1.900000000000 L= 0.349000000000
C048: H=-1 K= 0.100000000000 G= 1.800000000000 L=-0.208000000000
C048: H= 1 K= 0.010000000000 G= 1.810000000000 L=-0.156359000000
C048: H= 1 K= 0.010000000000 G= 1.820000000000 L=-0.103832000000
C048: H= 1 K= 0.010000000000 G= 1.830000000000 L=-0.050413000000
C048: H= 1 K= 0.010000000000 G= 1.840000000000 L= 0.003904000000
C048: H=-1 K= 0.001000000000 G= 1.839000000000 L=-0.001568281000
C048: H= 1 K= 0.000100000000 G= 1.839100000000 L=-0.001021459529
C048: H= 1 K= 0.000100000000 G= 1.839200000000 L=-0.000474547712
C048: H= 1 K= 0.000100000000 G= 1.839300000000 L= 0.000072454457
C048: H=-1 K= 0.000010000000 G= 1.839290000000 L= 0.000017750174
C048: H=-1 K= 0.000010000000 G= 1.839280000000 L=-0.000036953205
C048: H= 1 K= 0.000001000000 G= 1.839281000000 L=-0.000031482908
C048: H= 1 K= 0.000001000000 G= 1.839282000000 L=-0.000026012602
C048: H= 1 K= 0.000001000000 G= 1.839283000000 L=-0.000020542286
C048: H= 1 K= 0.000001000000 G= 1.839284000000 L=-0.000015071962
C048: H= 1 K= 0.000001000000 G= 1.839285000000 L=-0.000009601629
C048: H= 1 K= 0.000001000000 G= 1.839286000000 L=-0.000004131286
C048: H= 1 K= 0.000001000000 G= 1.839287000000 L= 0.000001339065
C048: H=-1 K= 0.000000100000 G= 1.839286900000 L= 0.000000792030

And we end up with:

G= 1.839286900000


But this made me think, what happens if the solution gets closer and closer to 1e-6.

cbeq(1e3, -1e2, -1e1, -1)

Code:

C025: X= 0.000010000000 Y= 0.101000000000 True
C025: X= 0.000100000000 Y= 0.101000000000 True
C025: X= 0.001000000000 Y= 0.101000000000 True
C025: X= 0.010000000000 Y= 0.101000000000 True
C025: X= 0.100000000000 Y= 0.101000000000 True
C025: X= 1.000000000000 Y= 0.101000000000 False
C048: H= 1 K= 0.100000000000 G= 0.100000000000 L=-0.002000000000
C048: H= 1 K= 0.100000000000 G= 0.200000000000 L= 0.001000000000
C048: H=-1 K= 0.010000000000 G= 0.190000000000 L= 0.000349000000
C048: H=-1 K= 0.010000000000 G= 0.180000000000 L=-0.000208000000
C048: H= 1 K= 0.001000000000 G= 0.181000000000 L=-0.000156359000
C048: H= 1 K= 0.001000000000 G= 0.182000000000 L=-0.000103832000
C048: H= 1 K= 0.001000000000 G= 0.183000000000 L=-0.000050413000
C048: H= 1 K= 0.001000000000 G= 0.184000000000 L= 0.000003904000
C048: H=-1 K= 0.000100000000 G= 0.183900000000 L=-0.000001568281
C048: H= 1 K= 0.000010000000 G= 0.183910000000 L=-0.000001021460
C048: H= 1 K= 0.000010000000 G= 0.183920000000 L=-0.000000474548

G= 0.183920000000

---

cbeq(1e6, -1e4, -1e2, -1)

Code:

C025: X= 0.000010000000 Y= 0.010001000000 True
C025: X= 0.000100000000 Y= 0.010001000000 True
C025: X= 0.001000000000 Y= 0.010001000000 True
C025: X= 0.010000000000 Y= 0.010001000000 True
C025: X= 0.100000000000 Y= 0.010001000000 False
C048: H= 1 K= 0.010000000000 G= 0.010000000000 L=-0.000002000000
C048: H= 1 K= 0.010000000000 G= 0.020000000000 L= 0.000001000000
C048: H=-1 K= 0.001000000000 G= 0.019000000000 L= 0.000000349000

G= 0.019000000000

---

cbeq(1e9, -1e6, -1e3, -1)

Code:

C025: X= 0.000010000000 Y= 0.001000001000 True
C025: X= 0.000100000000 Y= 0.001000001000 True
C025: X= 0.001000000000 Y= 0.001000001000 True
C025: X= 0.010000000000 Y= 0.001000001000 False
C048: H= 1 K= 0.001000000000 G= 0.001000000000 L=-0.000000002000

G= 0.001000000000


Unfortunately it looks like fewer and fewer digits are calculated.

Long story short, the algorithm to calculate the first, always real root, is an open numerical method.

We start with a guess, and decimals tlll 10E-6 are calculated.

I have also tried to use native solver of HP35S. The program was shorter, but sometimes the accuracy was poor, sometime execution time was too long, or both.

I have used the constant value 1E-6 in order to balance speed and accuracy. I consider extremely important to obtain the first root in a decent time, but those who prefer a different value may change instruction C023 as they like.

[Thomas Klemm]

With little extra effort, we can also use Halley's method:

\(
\begin{align}
x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})-{\frac {f(x_{n})}{f'(x_{n})}}{\frac {f''(x_{n})}{2}}}}
\end{align}
\)

The advantage is that the rate of convergence to the root is cubic.

Code:

00 { 80-Byte Prgm }
01▸LBL "CBEQ"
02 R↑
03 STO÷ ST T
04 STO÷ ST Z
05 ÷
06 STO 02
07 R↑
08 STO 00
09 X<> ST Z
10 STO 01
11 ÷
12 +/-
13▸LBL 00
14 ENTER
15 ENTER
16 ENTER
17 RCL+ 00
18 STO 03
19 ×
20 RCL+ 01
21 STO 04
22 ×
23 RCL+ 02
24 X<> ST Z
25 RCL+ 03
26 STO+ ST Y
27 RCL× ST T
28 RCL+ 04
29 STO÷ ST Y
30 X<>Y
31 RCL× ST Z
32 -
33 ÷
34 -
35 X≠Y?
36 GTO 00
37 STOP
38 RCL 03
39 -2
40 ÷
41 ENTER
42 ENTER
43 X↑2
44 RCL- 04
45 SQRT
46 STO- ST Z
47 +
48 END

Example

But unfortunately the previous example \(x^3-x^2-x-1=0\) doesn't work anymore:

1
ENTER
-1
ENTER
ENTER
XEQ "CBEQ"

Divide by 0
x: 0

In such a case, you can just use a different initial value and start the inner loop with:

2
XEQ 00

1.83928675521


Feel free to include a RND command before this line:

Code:

35 X≠Y?

This avoids the endless loop that Gerson experienced.
Also it allows to limit the execution time when full accuracy is not needed.

These are the intermediate results of the previous example when rounded with FIX 12:

2
1.84090909091
1.83928675734
1.83928675521
1.83928675521