HP-45 Application Book (Not getting Correct Answer))

[Namir]

On Page 75 of the HP-45 Application Book, there are theory and instructions to solve for the roots of a cubic equation. I have tried both example. I get matching results with the first example, however I get different results with the second exmple. I also tried the calculation using MATLAB and they match the ones I got with the Free42. I get the same results which are:

d = 1.7778 (determinant)
x1 = 3.0000 + 0.5774i
x2 = 1.0000 + 0.5774i
x3 = 0.0000 - 1.1547i

While the book gives the following results for the secon example (x^3-4*x^2+8*x-8=0):

x1 = 2.00
x2 = 1 + 1.73i
x3 = 1 - 1.73i
(d = 1.78)

I am applying the equations in the HP-45 Application Book still getting different answers. Here is the Free-42 code:

Code:

01 LBL "CUBIC"
02 STO "C"
03 R↓
04 STO "B"
05 R↓
06 STO "A"
07 STO "A0"
08 X↑2
09 -3
10 ÷
11 RCL "B"
12 +
13 STO "P"
14 RCL "A"
15 3
16 ÷
17 3
18 Y↑X
19 2
20 ×
21 RCL "A"
22 RCL "B"
23 ×
24 3
25 ÷
26 -
27 RCL "C"
28 +
29 STO "Q"
30 X↑2
31 4
32 ÷
33 RCL "P"
34 3
35 Y↑X
36 27
37 ÷
38 +
39 STO "D"
40 SQRT
41 RCL "Q"
42 2
43 ÷
44 -
45 3
46 1/X
47 Y↑X
48 STO "A"
49 RCL "D"
50 SQRT
51 RCL "Q"
52 2
53 ÷
54 +
55 +/-
56 3
57 1/X
58 Y↑X
59 STO "B"
60 RCL "A"
61 +
62 STO "Y1"
63 RCL "A"
64 RCL "B"
65 -
66 2
67 ÷
68 3
69 SQRT
70 ×
71 STO "I1"
72 RCL "A"
73 RCL "B"
74 +
75 -2
76 ÷
77 STO "R1"
78 RCL "I1"
79 -1
80 SQRT
81 ×
82 +
83 STO "Y2"
84 RCL "R1"
85 RCL "I1"
86 -1
87 SQRT
88 ×
89 -
90 STO "Y3"
91 RCL "A0"
92 3
93 ÷
94 STO "A0"
95 -
96 RCL "Y2"
97 RCL "A0"
98 -
99 RCL "Y1"
100 RCL "A0"
101 -
102 RTN

The MATLAB command roots([1 -4 8 -8]) gives answers that match the HP-45 Application book. So, what am I doing wrong? You need to download the application book from Eric's website and go to page 74 and 75 to read the equations and the examples.

Namir

[Albert Chan]

Hi, Namir

I have not read the book, nor your code, but this is easy to explain.

³√(-8) = -2, but (-8)^(1/3) = -2*ω, where ω = cis(2*pi/3) = (-1 + i√3) / 2

Because of this ambiguity, below is a *bad* cubic formula.
x might not be a root al all! One term might require ω factor correction.



Instead of doing 2 cube roots, do only 1, and derive the other, guaranteed correct one.
Bonus: we can get the one with least cancellation errors, and have both accurate terms.

x^3 - 4*x^2 + 8*x - 8 = 0            // let x = y + 4/3
y^3 + 8/3*y - 56/27 = 0               // c = 8/3, d = 56/27, to match above formula

--> (y - (α + β)) * (y - (α*ω + β/ω)) * (y - (α/ω + β*ω)) = 0

α = ³√(d/2 + √((d/2)^2 + (c/3)^3)) = ³√(28/27 + 4/3) = ³√(64/27) = 4/3

β = (-c/3) / α = (-8/9) / (4/3) = -2/3

β is correctly paired with α --> all generated roots are correct.

x = 4/3 + y

x1 = 4/3 + (4/3)     + (-2/3)     = 2
x2 = 4/3 + (4/3)*ω + (-2/3)/ω = 1 + i√3
x3 = 4/3 + (4/3)/ω + (-2/3)*ω = -1/3 − i√3

Had we use cube root for β, we may get (28/27 - 4/3)^(1/3) = -2/3/ω

XCas> ω := exp(i*2*pi/3)
XCas> float(normal([[1,1],[ω,1/ω],[1/ω,ω]] * [4/3, -2/3/ω] .+ 4/3))

[3.0+0.57735026919*i, 1.0+0.57735026919*i, -1.15470053838*i]

[Thomas Klemm]

(06-28-2024 10:18 PM)Namir Wrote:  

Code:

49 RCL "D"
50 SQRT
51 RCL "Q"
52 2
53 ÷
54 +
55 +/-
56 3
57 1/X
58 Y↑X
59 STO "B"

So, what am I doing wrong?

While calculating B in line 56, the x-value is negative.
Calculating the cubic root returns the principal value.
This is a complex number.

What you want to calculate instead is the negative of the cubic root of the negative value.

This is handled in the HP-45 Application Book.

[Namir]

Thanks Albert and Thomas!

Namir

[Namir]

Good news! I was able to adapt the HP-45 listing for the HP-41C and posted the adaptation in the HP41C Software Library section. The solution uses flags and labels to correctly direct the flow of calculations.

Namir