So I stumbled upon this section of the Wikipedia article concerning cube root. This formula is used to calculate it on a simple four-banger:

Of course I tried to figure out how to do it with RPN. Here's the program I came up with:

{ 60-Byte Prgm }

LBL "CURT"

XEQ 01

ENTER

XEQ 02

*

ENTER

XEQ 03

*

ENTER

XEQ 04

*

ENTER

XEQ 05

*

ENTER

XEQ 06

*

GTO 01

LBL 06

XEQ 05

LBL 05

XEQ 04

LBL 04

XEQ 03

LBL 03

XEQ 02

LBL 02

XEQ 01

LBL 01

SQRT

END

Is this algorithm well known? I'm just curious because I've never heard of it before.

Cheers

Thomas

Edited: 28 Oct 2012, 4:20 p.m.

This looks like an algorithm I came across many years ago in some manual for a simple square root calculator, maybe a TI:

               Cube root

Key in number

Store in memory

Take square root

then repeat until no change:

  X MR = [sqrt] [sqrt]

It appears to use the following formula instead (another variant to write the geometric series):


Or looking at it from another point of view, it's a fixed point iteration for

However the advantage of the previous method is:

Quote:

---

No memory is required.

---

Kind regards

Thomas

Thomas,

I like your fixed point iteration algorithm. It is very simple, it works, but is slower than Newton's method (as shown the cited Wikipedia article).

Namir

Transforming your fixed iteration point into:

f(x) = (ax)^k - x

Where k = 1/4. We then have:

f'(x) = k a (ax)^(k-1)

and

f''(x) = k (k-1) a^2 (ax)^(k-2)

Using Halley's method, you get very fast convergence. This is faster than Newton and Halley methods applied to the original cubic form (as show in the reference Wikipedia article).

Edited: 29 Oct 2012, 8:51 a.m.

Thomas,

The Wikipedia article you reference has an Alternative Method that uses the following equation:

x = 4/3*(ax)^1/4 - x/3 = 4/3*sqrt(sqrt(a*x)) - x/3

The accompanying text says that this method works very fast. I tried it and it bombed out since the first iteration gave a negative x!! Can someone else try it????????????

Namir

Edited: 29 Oct 2012, 4:43 p.m.

Hi Namir

Just tested it with 10 and 0.5:

>>> from math import sqrt

>>> def curt(a):

...   x = 1

...   for i in range(10):

...     x = (4*sqrt(sqrt(a*x))-x)/3

...     print x

... 

>>> curt(10)

2.03770588005

2.15361817766

2.15443465134

2.15443469003

2.15443469003

2.15443469003

2.15443469003

2.15443469003

2.15443469003

2.15443469003

>>> curt(0.5)

0.787861887005

0.793695134037

0.79370052598

0.793700525984

0.793700525984

0.793700525984

0.793700525984

0.793700525984

0.793700525984

0.793700525984

It seems to converge quite fast. What was your example?

Kind regards

Thomas

I was not so much interested in the speed of convergence of the algorithm but more amused to find something similar to calculate the Dottie Number: just punch some keys repeatedly and wonder what's going to happen and why.

Yes the method works fine AND fast for the cube root of 10, 0.5, and even 12^3. Try say 123^3. I get a negative number after the first iteration! So it seems that there is a limit.

After comparing various methods, using Excel, i found that the best iterative method is Halley's method applied to:

f(x) = (ax)^1/4 - x

Namir

Edited: 29 Oct 2012, 6:39 p.m. after one or more responses were posted

I get the expected result:

>>> 123**3    

1860867

>>> curt(1860867)

48.9122822015

113.929207914

122.91259267

122.999992233

123.0

123.0

123.0

123.0

123.0

123.0

What's your starting guess for x?

I use sqrt(cube_number) as the initial guess for all the methods I am comparing (except the calculator method). What is your initial guess?

Edited: 29 Oct 2012, 6:46 p.m.

When I try a low initial guess like 1, the method finds the cube root for large cubes like 123^3. The initial guess of 1 gives me the same sequence of numbers you got!!!! Cool!!!!!!!!

Edited: 29 Oct 2012, 6:47 p.m.

What is your initial guess for cube > 1 and for cube < 1?

In the article it says:

Quote:

---

Because of the fast converging an initial approximation of 1 suffices.

---

That's exactly what I used.

The initial guess of 1 does not work for cubes that are very small ... try the cube root of 5.37384E-07.

I know ... I know ... the devil is in the details.

Namir

Edited: 29 Oct 2012, 8:36 p.m.

Good catch! But that's easy to fix: just set x = a if a < 1.

>>> from math import sqrt

>>> def curt(a):

...   x = 1

...   if a < 1:

...     x = a

...   for i in range(10):

...     x = (4*sqrt(sqrt(a*x))-x)/3

...     print x

... 

>>> curt(5.37384E-07)

0.000977240604425

0.00605703930305

0.00805204305187

0.00812998755053

0.00813008171241

0.00813008171254

0.00813008171254

0.00813008171254

0.00813008171254

0.00813008171254

If I did my math correctly the limit in the condition is not 1 but 1/256.

Thanks Thomas!!

Edited: 30 Oct 2012, 10:13 a.m.

You mean like repeating:

X^2

X^2

LN

With an initial value of 8 or 9.

You can execute the X^2 function more than twice and get different converging numbers.

Namir

Edited: 30 Oct 2012, 10:17 a.m.

I tried this on a machine that uses single-precision IEEE floating point (32-bit), which has a range of about 1E+/-38 (and a bit smaller for non-normalized numbers). It works for input values from about 7.0E-25 to 4.0E22 (if I remember correctly -- did the test at work and I'm home now). Any smaller and the value diverges toward zero, and any bigger and the first iteration gives "infinity" (+1.$), and, so, the second gives "Not a Number".

I use the square root as the starting guess. Doing this, I can find some "goofy" extreme values that take maybe 14 or 15 loops to converge, but most "reasonable" numbers take 4 or 5. And the number it converges to is more accurate (in the cases I tested) than doing input ** (1.0 / 3.0) on the same machine. (I checked by taking result * result * result; some were easily checked by inspection: cuberoot(1000.0) converges to exactly 10.0 vs 9.9999...)

Unfortunately, I didn't have time to compare execution time.

Thanks for everybody's input on this thread. I learned a lot, and I'm saving this one in my "toolkit".

Dale

x = log(x4) = 4 log x

- x/4 = - log x

e- x/4 = 1/x

- x/4 e- x/4 = - 1/4

- x/4 = W(- 1/4)

x = - 4 W(- 1/4)

x = 8.61316945644

Nice! Let's have another thread about the
Lambert W function.

Cheers

Thomas

Using with n = 2¹⁰ = 1024 we can do that even on a 4-banger without logarithm:

00 { 17-Byte Prgm }

01 SQRT

02 SQRT

03 SQRT

04 SQRT

05 SQRT

06 SQRT

07 SQRT

08 SQRT

09 1

10 -

11 1024

12 *

13 END

It converges to 8.68125 so it's not very accurate. But still not bad.

Cheers

Thomas

I was able to figure it out that algorithm in 1977. I also was able to figure our the reverse algorithm that calculates ln(x) using the sqrt function.

Ahh .... old memories :-)

Hi

I try the "alternative method" with my emulator.

two versions (AOS et RPN)

It's working well... and fast.

Thank for this exercise.

RPN Version

LBL A
STO 02 1 STO 01 1 0 STO 00
LBL SQR
RCL 02 RCL 01 * SQR SQR 4 * RCL 01 - 3 / STO 01 PRT
DSZ 00 SQR
R/S

AOS Version

LBL A
STO 02 1 STO 01 1 0 STO 00
LBL SQR
( ( ( RCL 02 * RCL 01 ) SQR SQR * 4 ) - RCL 01 ) / 3 =
STO 01 PRT
DSZ 00 SQR
R/S

Example with 9993948264

421.2398826

1769.463631

2144.398571

2153.994636

2154

2154

2154

2154

2154

2154

Best regards
Pierre