Cube root on standard calculator



#2

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.


#3

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]

#4

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


#5

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

#6

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.


#7

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.


#8

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.


#9

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

#10

Using with n = 210 = 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


#11

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 :-)

#12

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.


#13

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


#14

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


#15

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?


#16

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.


#17

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

#18

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.

#19

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


#20

In the article it says:

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

That's exactly what I used.


#21

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.


#22

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.


#23

Thanks Thomas!!

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

#24

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP Prime] Using n-root symbol and exponent problem uklo 7 409 11-11-2013, 01:39 AM
Last Post: Alberto Candel
  Cubic root (-8) = 2 ? Gilles Carpentier 37 1,423 08-12-2013, 10:26 PM
Last Post: jep2276
  Square Root Simplifier for HP39gII Mic 4 269 03-11-2013, 08:25 AM
Last Post: Eddie W. Shore
  HP-71B - thanks to Marcus von Cube for MATH ROM article Michael Lopez 2 206 03-03-2013, 07:19 AM
Last Post: Paul Berger (Canada)
  ROOT bug? HP 48S/48G Eddie W. Shore 8 445 07-13-2012, 07:05 PM
Last Post: Eddie W. Shore
  Test HP-41C PPC ROM and STANDARD modules Robert (Simi Valley) 2 216 06-14-2012, 08:05 AM
Last Post: Frido Bohn
  HP-67 standard pac Wytnucls 3 220 04-13-2012, 11:23 AM
Last Post: Jim Horn
  x root y on hp42s David Griffith 14 596 04-08-2012, 12:43 PM
Last Post: Walter B
  HP 32sII Integration Error of Standard Normal Curve Anthony (USA) 4 254 03-14-2012, 03:18 AM
Last Post: Nick_S
  Calculators on Standard Tests Norman Dziedzic 18 641 02-14-2012, 10:26 PM
Last Post: bill platt

Forum Jump: