Posts: 735
Threads: 34
Joined: May 2007
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.
Posts: 33
Threads: 5
Joined: Jan 2009
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]
Posts: 735
Threads: 34
Joined: May 2007
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 735
Threads: 34
Joined: May 2007
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
Posts: 735
Threads: 34
Joined: May 2007
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 735
Threads: 34
Joined: May 2007
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?
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
What is your initial guess for cube > 1 and for cube < 1?
Posts: 735
Threads: 34
Joined: May 2007
In the article it says:
Quote:
Because of the fast converging an initial approximation of 1 suffices.
That's exactly what I used.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 735
Threads: 34
Joined: May 2007
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
Thanks Thomas!!
Edited: 30 Oct 2012, 10:13 a.m.
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 73
Threads: 2
Joined: Sep 2011
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
Posts: 735
Threads: 34
Joined: May 2007
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
Posts: 735
Threads: 34
Joined: May 2007
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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 :-)
Posts: 20
Threads: 4
Joined: Apr 2012
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
|