Cube root on standard calculator Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-28-2012, 04:16 PM 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. Bob Patton Junior Member Posts: 33 Threads: 5 Joined: Jan 2009 10-28-2012, 05:41 PM ```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] ``` Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-28-2012, 10:00 PM 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 Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 03:47 AM 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 Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 08:50 AM 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. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 04:42 PM 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. Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-29-2012, 05:56 PM 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 Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-29-2012, 06:31 PM 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. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 06:31 PM 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 Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-29-2012, 06:36 PM 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? Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 06:40 PM 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. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 06:44 PM 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. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 06:46 PM What is your initial guess for cube > 1 and for cube < 1? Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-29-2012, 07:08 PM In the article it says: Quote: Because of the fast converging an initial approximation of 1 suffices. That's exactly what I used. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2012, 07:35 PM 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. Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-29-2012, 08:47 PM 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. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-30-2012, 04:29 AM Thanks Thomas!! Edited: 30 Oct 2012, 10:13 a.m. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-30-2012, 10:16 AM 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. Dale Reed Member Posts: 73 Threads: 2 Joined: Sep 2011 10-30-2012, 08:53 PM 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 Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-30-2012, 10:55 PM ```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 Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-31-2012, 02:06 AM 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 Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-31-2012, 04:21 PM 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 :-) Pierre Junior Member Posts: 20 Threads: 4 Joined: Apr 2012 11-09-2012, 06:11 AM 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 « Next Oldest | Next Newest »

 Possibly Related Threads... Thread Author Replies Views Last Post [HP Prime] Using n-root symbol and exponent problem uklo 7 2,143 11-11-2013, 01:39 AM Last Post: Alberto Candel Cubic root (-8) = 2 ? Gilles Carpentier 37 7,353 08-12-2013, 10:26 PM Last Post: jep2276 Square Root Simplifier for HP39gII Mic 4 1,561 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 1,079 03-03-2013, 07:19 AM Last Post: Paul Berger (Canada) ROOT bug? HP 48S/48G Eddie W. Shore 8 3,174 07-13-2012, 07:05 PM Last Post: Eddie W. Shore Test HP-41C PPC ROM and STANDARD modules Robert (Simi Valley) 2 1,112 06-14-2012, 08:05 AM Last Post: Frido Bohn HP-67 standard pac Wytnucls 3 940 04-13-2012, 11:23 AM Last Post: Jim Horn x root y on hp42s David Griffith 14 3,362 04-08-2012, 12:43 PM Last Post: Walter B HP 32sII Integration Error of Standard Normal Curve Anthony (USA) 4 1,180 03-14-2012, 03:18 AM Last Post: Nick_S Calculators on Standard Tests Norman Dziedzic 18 3,310 02-14-2012, 10:26 PM Last Post: bill platt

Forum Jump: