HP35S - Cubic (and incidentally quadratic) equations solver

HP35S - Cubic (and incidentally quadratic) equations solver - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: HP35S - Cubic (and incidentally quadratic) equations solver (/thread-19783.html)

HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 04-07-2023 09:31 AM

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!

RE: HP35S - Cubic (and incidentally quadratic) equations solver - PedroLeiva - 04-07-2023 07:55 PM

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

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 04-08-2023 04:39 AM

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

RE: HP35S - Cubic (and incidentally quadratic) equations solver - PedroLeiva - 04-08-2023 11:46 AM

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

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 04-10-2023 02:21 PM

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.

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 04-20-2023 07:18 AM

(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!

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 04-20-2023 11:36 AM

(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

automatically, when executing f I or f Re<>Im; or

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̅

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 05-12-2023 07:49 AM

(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

automatically, when executing f I or f Re<>Im; or

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.

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 05-13-2023 10:33 AM

(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?

RE: HP35S - Cubic (and incidentally quadratic) equations solver - gor1060 - 05-13-2023 12:40 PM

[attachment=12106]

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Gerson W. Barbosa - 05-13-2023 06:00 PM

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

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 05-18-2023 06:26 PM

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.

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 05-20-2023 09:45 AM

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

Could it be that this should rather be:

C026 X<=Y?

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 05-20-2023 10:22 AM

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

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 05-20-2023 12:05 PM

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.

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 05-20-2023 12:34 PM

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.

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Roberto Volpi - 05-21-2023 10:49 AM

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

RE: HP35S - Cubic (and incidentally quadratic) equations solver - Thomas Klemm - 05-22-2023 07:53 AM

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