Cube root [HP-35]

[Gerson W. Barbosa]

27
√ √
ENTER
√ √ √ √
×
ENTER
√ √ √ √ √ √ √ √
×
ENTER
√ √
×

-> 2.999949713

https://youtu.be/gU_gty53GfU

PS: On the 35 we can use logs, but we can still use the algorithm on the HP-16C, which lacks them.

[EdS2]

I've a feeling we're using a truncated binary approximation to a third. Clever!

[Gerson W. Barbosa]

(03-06-2020 09:50 AM)EdS2 Wrote:  I've a feeling we're using a truncated binary approximation to a third. Clever!

I only converted the method in the YouTube video to RPN. It took me a while to figure out what he meant by "ten cross", though :-)

Algorithms using square roots to approximate non-available functions might have been interesting when scientific calculators were not affordable. There was one for ln(x) which I managed to improve a bit. See here.

[Gene]

Clever yes.

Is there a need to use it however? :-)

3 1/x 27 X^y ?

[Albert Chan]

This is a rare case which CHAIN take less keystrokes than RPN

I would swap the order of the keys though, so successive iterations improve the cube root

Code:

; CHAIN calculator, N^(1/3) = N^(5/15)
27                   ; N
√ √ ×                ; N^0x0.4    = 2.279507057
√ √ ×                ; N^0x0.5    = 2.800923042
√ √ √ √ ×            ; N^0x0.55   = 2.987153223
√ √ √ √ √ √ √ √ ×    ; N^0x0.5555 = 2.999949710

[Gerson W. Barbosa]

(03-06-2020 11:57 AM)Gene Wrote:  Is there a need to use it however? :-)

3 1/x 27 X^y ?

No, there is not, hence the PS in post #1. I chose the HP-35 because it’s my oldest RPN calculator with square root as a primary function key.

Is there an HP calculator other than the HP-16C this might be useful to?

[Gene]

ah, I saw that but it didn't hit me early this morning about the implementation of X^y. :-)

[Gerson W. Barbosa]

(03-06-2020 01:41 PM)Albert Chan Wrote:  I would swap the order of the keys though, so successive iterations improve the cube root

Code:

; CHAIN calculator, N^(1/3) = N^(5/15)
27                   ; N
√ √ ×                ; N^0x0.4    = 2.279507057
√ √ ×                ; N^0x0.5    = 2.800923042
√ √ √ √ ×            ; N^0x0.55   = 2.987153223
√ √ √ √ √ √ √ √ ×    ; N^0x0.5555 = 2.999949710

Nice change!

Another method, but more keystrokes, at least on the 8-digit Canon LC-37:

27
√ √ √ √ √ √ √ √ √ √ √
-
1 / 3
+ 1 =
× = × = × = × = × =
× = × = × = × = × =
× =
-> 3.000989

Replace 3 with 5 and you’ll get the 5th root: 1.9338459

From https://www.quora.com/Using-a-simple-cal...root-n-odd

PS: Again, it looks like they have trouble with /th/ – “...den press √ (12 times)” :-)

[Juan14]

Fixed point iteration:

From x^3 – A = 0 we have: x^2 = A/x and x = √(A/x), we can use the last equation for the fixed point iteration method.
For example for A=27 and (with x0 = 27) we load the stack with 27.

27 ENTER ENTER ENTER

for the first steep we press

√ √

then for each iteration we press

/ √

after 10 iterations we have: x = 2.99919546023

[Albert Chan]

(03-08-2020 03:23 PM)Juan14 Wrote:  27 ENTER ENTER ENTER

for the first steep we press

√ √

then for each iteration we press

/ √

after 10 iterations we have: x = 2.99919546023

Slight improvement to above cube roots fixed point iteration, by replacing "/ √" as "× √ √"
This cut down iteration counts in half.

[Gerson W. Barbosa]

(03-08-2020 03:23 PM)Juan14 Wrote:  27 ENTER ENTER ENTER

for the first steep we press

√ √

then for each iteration we press

/ √

after 10 iterations we have: x = 2.99919546023

Great!

After 9 iterations we have

3.001609727

Now we store that result and do a final iteration

STO / √

then

ENTER + RCL + 3 /

and we get

3.000000216

(on the HP-35)

[Gerson W. Barbosa]

just a proof of concept on the hp 33s to reduce the number of iterations in this variation of Juan’s algorithm. Line D0013 should be replaced by a short routine to compute ln(x) on the HP-16C, where this might make some sense.

Code:


C0001 LBL C
C0002 12
C0003 STO i
C0004 x<>y
C0005 ENTER 
C0006 ENTER
C0007 ENTER
C0008 √x
D0001 LBL D
D0002 ÷
D0003 √x
D0004 DSE i
D0005 GTO D
D0006 LASTx
D0007 R↑
D0008 x<>y
D0009 ÷
D0010 x<>y
D0011 8192
D0012 R↑
D0013 LN
D0014 ÷
D0015 5..18 ; 5/18 = 2.77777777778E-1
D0016 +
D0017 1/x
D0018 +/-
D0019 2
D0020 +
D0021 ×
D0022 LASTx
D0023 R↓
D0024 +
D0025 R↑
D0026 1
D0027 +
D0028 ÷
D0029 RTN

27 XEQ C -> 3

1E50 XEQ C -> 4.64158883316E16

1E50 ³√x -> 4.64158883316E16

Algorithm:

r₀ = √x
r₁ = √(x/r₀)
r₂ = √(x/r₁)
...
rn = √(x/rn-1)

a = rn
b = rn-1

³√x ~ {a[2 - 1/(2ⁿ⁺¹/lnx + 5/18)] + b}/[3 - 1/(2ⁿ⁺¹/lnx + 5/18)]

[Albert Chan]

We can use Pade approximation to speedup convergence.
Example, simple Pade[1,1] gives cubic convergence.

Pade[(1+z)^(1/n), {z,0,1,1}] = \(\large{2n + (n+1)z \over 2n + (n-1)z}\)

\(\large\sqrt[n]k = x \left(\sqrt[n] {1 + ({k \over x^n} - 1)}\right)
≈ x \left({k(n+1)\;+\; (n-1)x^n \over k(n-1)\;+\;(n+1)x^n } \right) \)

For cube roots, we have: \(\large\sqrt[3]k ≈ x \left({2k\;+\; x^3 \over k\;+\;2x^3 } \right) \)

Example, cube root of 0.1 (guess x=0.5):

0.5
0.46428 57142 85714 28571
0.46415 88833 67588 52147
0.46415 88833 61277 88924 10076 35091 94543
0.46415 88833 61277 88924 10076 35091 94465 76551 34912 50112 ... (98 correct digits)

Above formula is *exactly* the same as Halley's method, H_f(x), with f(x) = x^n - k
see Method of obtaining Third Order Iterative Formulas

[Gerson W. Barbosa]

(03-16-2020 02:42 PM)Albert Chan Wrote:  We can use Pade approximation to speedup convergence.

...

For cube roots, we have: \(\large\sqrt[3]k ≈ x \left({2k\;+\; x^3 \over k\;+\;2x^3 } \right) \)

Very nice! For k = 27 and x₀ = √√27 I get 3.00000000001 after 3 iterations. Thanks!

[Albert Chan]

(03-16-2020 07:49 PM)Gerson W. Barbosa Wrote:  

(03-16-2020 02:42 PM)Albert Chan Wrote:  We can use Pade approximation to speedup convergence.

...

For cube roots, we have: \(\large\sqrt[3]k ≈ x \left({2k\;+\; x^3 \over k\;+\;2x^3 } \right) \)

Very nice! For k = 27 and x₀ = √√27 I get 3.00000000001 after 3 iterations. Thanks!

One issue is optimal range of k, to maximize convergence rate.
In other words, with guess = k^y, y≠1/3, should we convert \(\sqrt[3]{27} = 10 \sqrt[3]{0.027}\) ?

At the break-even point, errors should match, but with opposite sign:

\(\sqrt[3]k - k^y = 10(k/1000)^y - \sqrt[3]k\)
\(2 \sqrt[3]k = (1 + 10/1000^y) k^y \)
At break-even, \(\large k = (1/2 + 5/1000^y)^{1 \over 1/3\;-\;y}\)

With guess, x₀ = √√k = k^(1/4), break-even k ≈ 51.636 → optimal k in range [0.051636, 51.636)

Example, try \(\;\sqrt[3]{100}\) = 4.64158 88336 12728 89241 00763 ...

x₀ = √√(100) → x₃ = 4.64158 88336 12724 68166 58655
x₀ = 10√√0.1 → x₃ = 4.64158 88336 12728 89241 08943

Keeping k inside optimal range, \(\sqrt[3]{100} = 10\;\sqrt[3]{0.1}\), x₃ gained 7 digits accuracy.

---

OP, we estimate k^(1/3) as k^0x0.5555, break-even k ≈ 31.624 → optimal k in range [0.031624, 31.624)

       100^0x0.5555 = 4.641480114, error ≈ +0.000109
10 * 0.1^0x0.5555 = 4.641643194, error ≈ −0.000054, better, as expected

[Albert Chan]

(03-16-2020 10:54 PM)Albert Chan Wrote:  One issue is optimal range of k, to maximize convergence rate.
In other words, with guess = k^y, y≠1/3, should we convert \(\sqrt[3]{27} = 10 \sqrt[3]{0.027}\) ?

At the break-even point, errors should match, but with opposite sign:

\(\sqrt[3]k - k^y = 10(k/1000)^y - \sqrt[3]k\)
\(2 \sqrt[3]k = (1 + 10/1000^y) k^y \)
At break-even, \(\large k = (1/2 + 5/1000^y)^{1 \over 1/3\;-\;y}\)

I made an error assuming better guess speedup Pade convergence.
Plotting the errors suggest my pade setup prefer over-estimated guess.

For guess x₀ = √√k = k^(1/4), optimal k is before break-even point, at k ≈ √1000
Thus, optimal k should be in range [0.031623, 31.623)

\(\large\sqrt[3]k ≈ x \left({2k\;+\; x^3 \over k\;+\;2x^3 } \right) \)

Example, try \(\quad\sqrt[3]{27} =\; 3\)
x₀ = √√(27)   → x₃ = 2.9999 99999 99999 99934     // error =  +66 ULP
x₀ = 10√√.027 → x₃ = 3.0000 00000 00000 00709     // error = -709 ULP

Example, try \(\quad\sqrt[3]{32}\) = 3.1748 02103 93639 89495 ...
x₀ = √√(32)   → x₃ = 3.1748 02103 93639 89237     // error = +258 ULP
x₀ = 10√√.032 → x₃ = 3.1748 02103 93639 89710     // error = -215 ULP

Edit: rough estimate, with k in range [0.031623, 31.623), x₀ = √√k :

\(\large\sqrt[3]k ≈ x_1 = {½\;+\; x_0 \over ½\;+\; 1/x_0}\)     // relative error < 1.5%

[Gerson W. Barbosa]

(03-17-2020 04:17 PM)Albert Chan Wrote:  Edit: rough estimate, with k in range [0.031623, 31.623), x₀ = √√k :

\(\large\sqrt[3]k ≈ x_1 = {½\;+\; x_0 \over ½\;+\; 1/x_0}\)     // relative error < 1.5%

My particular HP-16C application requires the evaluation of cubic roots in a very narrow range around 17/4, so I use

\(1+\sqrt{\frac{{x}-2}{6}}\), exact at k = 8

(Taken from http://ajmonline.org/2008/5.pdf)

Code:



067- LBL 3
068- 2
069- STO I
070- x⇆y
071- STO 0
072- ENTER 
073- +
074- STO 2
075- LASTx
076- 2
077- -
078- 6
079- /
080- √x
081- 1
082- +
083- LBL 4
084- ENTER 
085- ENTER 
086- ×
087- LASTx
088- ×
089- RCL 2
090- x⇆y
091- +
092- LASTx 
093- ENTER 
094- +
095- RCL 0
096- +
097- /
098- ×
099- DSZ
100- GTO 4
101- RTN

4.1 GSB 3 -> 1.600520664

[Albert Chan]

(03-20-2020 02:10 PM)Gerson W. Barbosa Wrote:  (Taken from http://ajmonline.org/2008/5.pdf)

Thanks for the link.

Setup as interation formula, this also have cubic convergence, slightly better than Pade[1,1]

\(\large \sqrt[3]k ≈\left(x + \sqrt{4k\;-\; x^3 \over 3x}\right) ÷ 2\)

Note: guess x can be at most 4^(1/3)-1 ≈ 58% above true value of \(\sqrt[3]k\)

Code:

STO 0    ; HP-12C code for cube root
Enter
Enter
×
×
CHS
X<>Y
4
×
+
RCL 0
3
×
/
SQRT
RCL 0
+
2
/

4.1 Enter 2
R/S     → 1.591607979
R/S     → 1.600520756
R/S     → 1.600520664

100 Enter 5
R/S     → 4.640872096
R/S     → 4.641588834

[Gerson W. Barbosa]

(03-20-2020 05:36 PM)Albert Chan Wrote:  Setup as interation formula, this also have cubic convergence, slightly better than Pade[1,1]

\(\large \sqrt[3]k ≈\left(x + \sqrt{4k\;-\; x^3 \over 3x}\right) ÷ 2\)

It’s just perfect for my purpose. Thanks!

Now I can get logs to 7 places on the 16C, like I used to do on my old logarithm tables book :-)

2 GSB A -> 0.6931471(360)
GSB B -> 2.0000(13165)

12345 GSB A -> 9.421006(848)
GSB B -> 12345.0(1957)

6.789 EEX 79 GSB A -> 183.8195(343) GSB B -> 6.(686887E79)

0.12345 GSB A -> -2.091919(104)
GSB B -> 0.123450(119)

230 GSB B -> 7.496895E99
GSB A -> 229.97041(15)

1 GSB A -> 0.000000000
GSB B -> 1.000000000
GSB B -> 2.7182(92170)
GSB B -> 15.1544(4181)
GSB A -> 2.71829(4016)
GSB A -> 1.000004(736)
GSB A -> 0.000004736

0.693147181 CHS GSB B -> 0.500000(16)

10 GSB A -> 2.302585(344) STO 2

2 GSB A RCL 2 / -> 0.3010299(435)

0 GSB A -> Error 0

2 CHS GSB A -> Error 0

Code:


001- LBL A; LN
002- 1
003- EEX
004- CHS
005- 2
006- CF 1
007- 1
008- +
009- STO I
010- CLx
011- LASTx
012- x⇆y
013- x≤y
014- GTO 2
015- LBL 0
016- √x
017- RCL I
018- x>y
019- GSB 1
020- R↓
021- x⇆y
022- ENTER 
023- +
024- x⇆y
025- GTO 0
026- RTN
027- LBL 1
028- R↓
029- STO 0
030- R↓
031- 2
032- ×
033- STO 1
034- RCL 0
035- ENTER 
036- ×
037- 9
038- ×
039- 6
040- RCL 0
041- ×
042- -
043- 2
044- +
045- √x
046- 3
047- RCL 0
048- ×
049- +
050- 1
051- -
052- GSB 3
053- ENTER 
054- 1/x
055- -
056- 1
057- -
058- RCL 1
059- ×
060- F?1
061- CHS
062- RTN
063- LBL 2
064- SF1
065- 1/x
066- RTN 
067- LBL 3
068- 2
069- STO I
070- ENTER 
071- +
072- x⇆y
073- ×
074- STO 0
075- LASTx
076- 2
077- -
078- 6
079- /
080- √x
081- 1
082- +
083- LBL 4
084- ENTER 
085- ENTER 
086- RCL 0
087- x⇆y
088- /
089- LASTx
090- ENTER 
091- ×
092- -
093- 3
094- /
095- √x
096- +
097- 2
098- /
099- DSZ
100- GTO 4
101- RTN
102- LBL B; eᵡ
103- 1
104- 3
105- STO I
106- R↓
107- 8
108- 1
109- 9
110- 2
111- +
112- LASTx
113- /
114- ENTER
115- ×
116- 1
117- +
118- 2
119- /
120- LBL 5
121- ENTER 
122- ×
123- DSZ
124- GTO 5
125- RTN

Edited to fix errors in line 019 and 045