HP35S - Cubic (and incidentally quadratic) equations solver

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