HP35S - Cubic (and incidentally quadratic) equations solver

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