HP35S - Cubic (and incidentally quadratic) equations solver
04-07-2023, 09:31 AM (This post was last modified: 05-20-2023 10:17 AM by Roberto Volpi.)
Post: #1
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
HP35S - Cubic (and incidentally quadratic) equations solver
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

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!

Put a calculator into your life!
04-07-2023, 07:55 PM
Post: #2
 PedroLeiva Member Posts: 214 Joined: Jun 2014
RE: HP35S - Cubic (and incidentally quadratic) equations solver
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
04-08-2023, 04:39 AM
Post: #3
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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.

Put a calculator into your life!
04-08-2023, 11:46 AM
Post: #4
 PedroLeiva Member Posts: 214 Joined: Jun 2014
RE: HP35S - Cubic (and incidentally quadratic) equations solver
Excellent Roberto, the program works very well. Thanks a lot. Pedro
04-10-2023, 02:21 PM
Post: #5
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
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.
04-20-2023, 07:18 AM
Post: #6
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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!

Put a calculator into your life!
04-20-2023, 11:36 AM
Post: #7
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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̅
05-12-2023, 07:49 AM (This post was last modified: 05-12-2023 07:50 AM by Roberto Volpi.)
Post: #8
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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.

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.

Put a calculator into your life!
05-13-2023, 10:33 AM
Post: #9
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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?
05-13-2023, 12:40 PM
Post: #10
 gor1060 Junior Member Posts: 23 Joined: Mar 2017
RE: HP35S - Cubic (and incidentally quadratic) equations solver

05-13-2023, 06:00 PM (This post was last modified: 05-13-2023 06:06 PM by Gerson W. Barbosa.)
Post: #11
 Gerson W. Barbosa Senior Member Posts: 1,580 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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.
05-18-2023, 06:26 PM (This post was last modified: 05-18-2023 06:26 PM by Roberto Volpi.)
Post: #12
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
RE: HP35S - Cubic (and incidentally quadratic) equations solver

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.

Put a calculator into your life!
05-20-2023, 09:45 AM
Post: #13
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(04-07-2023 09:31 AM)Roberto Volpi Wrote:  C026 X>=Y?

Could it be that this should rather be:

C026 X<=Y?
05-20-2023, 10:22 AM
Post: #14
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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.

Put a calculator into your life!
05-20-2023, 12:05 PM
Post: #15
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
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.
05-20-2023, 12:34 PM (This post was last modified: 05-20-2023 12:46 PM by Thomas Klemm.)
Post: #16
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
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.
05-21-2023, 10:49 AM
Post: #17
 Roberto Volpi Member Posts: 80 Joined: Feb 2021
RE: HP35S - Cubic (and incidentally quadratic) equations solver
(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.

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.

Put a calculator into your life!
05-22-2023, 07:53 AM
Post: #18
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
RE: HP35S - Cubic (and incidentally quadratic) equations solver
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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)