HP35S - Cubic (and incidentally quadratic) equations solver

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