blob: 37c8844bb211477c6a9da42b9f621c9e7078af79 [file] [log] [blame]
caryclark@google.com639df892012-01-10 21:46:10 +00001// http://metamerist.com/cbrt/CubeRoot.cpp
2//
3
4#include <math.h>
caryclark@google.com27accef2012-01-25 18:57:23 +00005#include "CubicUtilities.h"
caryclark@google.com639df892012-01-10 21:46:10 +00006
7#define TEST_ALTERNATIVES 0
8#if TEST_ALTERNATIVES
9typedef float (*cuberootfnf) (float);
10typedef double (*cuberootfnd) (double);
11
12// estimate bits of precision (32-bit float case)
13inline int bits_of_precision(float a, float b)
14{
15 const double kd = 1.0 / log(2.0);
16
17 if (a==b)
18 return 23;
19
20 const double kdmin = pow(2.0, -23.0);
21
22 double d = fabs(a-b);
23 if (d < kdmin)
24 return 23;
25
26 return int(-log(d)*kd);
27}
28
29// estiamte bits of precision (64-bit double case)
30inline int bits_of_precision(double a, double b)
31{
32 const double kd = 1.0 / log(2.0);
33
34 if (a==b)
35 return 52;
36
37 const double kdmin = pow(2.0, -52.0);
38
39 double d = fabs(a-b);
40 if (d < kdmin)
41 return 52;
42
43 return int(-log(d)*kd);
44}
45
46// cube root via x^(1/3)
47static float pow_cbrtf(float x)
48{
49 return (float) pow(x, 1.0f/3.0f);
50}
51
52// cube root via x^(1/3)
53static double pow_cbrtd(double x)
54{
55 return pow(x, 1.0/3.0);
56}
57
58// cube root approximation using bit hack for 32-bit float
59static float cbrt_5f(float f)
60{
61 unsigned int* p = (unsigned int *) &f;
62 *p = *p/3 + 709921077;
63 return f;
64}
65#endif
66
67// cube root approximation using bit hack for 64-bit float
68// adapted from Kahan's cbrt
69static double cbrt_5d(double d)
70{
71 const unsigned int B1 = 715094163;
72 double t = 0.0;
73 unsigned int* pt = (unsigned int*) &t;
74 unsigned int* px = (unsigned int*) &d;
75 pt[1]=px[1]/3+B1;
76 return t;
77}
78
79#if TEST_ALTERNATIVES
80// cube root approximation using bit hack for 64-bit float
81// adapted from Kahan's cbrt
82#if 0
83static double quint_5d(double d)
84{
85 return sqrt(sqrt(d));
86
87 const unsigned int B1 = 71509416*5/3;
88 double t = 0.0;
89 unsigned int* pt = (unsigned int*) &t;
90 unsigned int* px = (unsigned int*) &d;
91 pt[1]=px[1]/5+B1;
92 return t;
93}
94#endif
95
96// iterative cube root approximation using Halley's method (float)
97static float cbrta_halleyf(const float a, const float R)
98{
99 const float a3 = a*a*a;
100 const float b= a * (a3 + R + R) / (a3 + a3 + R);
101 return b;
102}
103#endif
104
105// iterative cube root approximation using Halley's method (double)
106static double cbrta_halleyd(const double a, const double R)
107{
108 const double a3 = a*a*a;
109 const double b= a * (a3 + R + R) / (a3 + a3 + R);
110 return b;
111}
112
113#if TEST_ALTERNATIVES
114// iterative cube root approximation using Newton's method (float)
115static float cbrta_newtonf(const float a, const float x)
116{
117// return (1.0 / 3.0) * ((a + a) + x / (a * a));
118 return a - (1.0f / 3.0f) * (a - x / (a*a));
119}
120
121// iterative cube root approximation using Newton's method (double)
122static double cbrta_newtond(const double a, const double x)
123{
124 return (1.0/3.0) * (x / (a*a) + 2*a);
125}
126
127// cube root approximation using 1 iteration of Halley's method (double)
128static double halley_cbrt1d(double d)
129{
130 double a = cbrt_5d(d);
131 return cbrta_halleyd(a, d);
132}
133
134// cube root approximation using 1 iteration of Halley's method (float)
135static float halley_cbrt1f(float d)
136{
137 float a = cbrt_5f(d);
138 return cbrta_halleyf(a, d);
139}
140
141// cube root approximation using 2 iterations of Halley's method (double)
142static double halley_cbrt2d(double d)
143{
144 double a = cbrt_5d(d);
145 a = cbrta_halleyd(a, d);
146 return cbrta_halleyd(a, d);
147}
148#endif
149
150// cube root approximation using 3 iterations of Halley's method (double)
151static double halley_cbrt3d(double d)
152{
153 double a = cbrt_5d(d);
154 a = cbrta_halleyd(a, d);
155 a = cbrta_halleyd(a, d);
156 return cbrta_halleyd(a, d);
157}
158
159#if TEST_ALTERNATIVES
160// cube root approximation using 2 iterations of Halley's method (float)
161static float halley_cbrt2f(float d)
162{
163 float a = cbrt_5f(d);
164 a = cbrta_halleyf(a, d);
165 return cbrta_halleyf(a, d);
166}
167
168// cube root approximation using 1 iteration of Newton's method (double)
169static double newton_cbrt1d(double d)
170{
171 double a = cbrt_5d(d);
172 return cbrta_newtond(a, d);
173}
174
175// cube root approximation using 2 iterations of Newton's method (double)
176static double newton_cbrt2d(double d)
177{
178 double a = cbrt_5d(d);
179 a = cbrta_newtond(a, d);
180 return cbrta_newtond(a, d);
181}
182
183// cube root approximation using 3 iterations of Newton's method (double)
184static double newton_cbrt3d(double d)
185{
186 double a = cbrt_5d(d);
187 a = cbrta_newtond(a, d);
188 a = cbrta_newtond(a, d);
189 return cbrta_newtond(a, d);
190}
191
192// cube root approximation using 4 iterations of Newton's method (double)
193static double newton_cbrt4d(double d)
194{
195 double a = cbrt_5d(d);
196 a = cbrta_newtond(a, d);
197 a = cbrta_newtond(a, d);
198 a = cbrta_newtond(a, d);
199 return cbrta_newtond(a, d);
200}
201
202// cube root approximation using 2 iterations of Newton's method (float)
203static float newton_cbrt1f(float d)
204{
205 float a = cbrt_5f(d);
206 return cbrta_newtonf(a, d);
207}
208
209// cube root approximation using 2 iterations of Newton's method (float)
210static float newton_cbrt2f(float d)
211{
212 float a = cbrt_5f(d);
213 a = cbrta_newtonf(a, d);
214 return cbrta_newtonf(a, d);
215}
216
217// cube root approximation using 3 iterations of Newton's method (float)
218static float newton_cbrt3f(float d)
219{
220 float a = cbrt_5f(d);
221 a = cbrta_newtonf(a, d);
222 a = cbrta_newtonf(a, d);
223 return cbrta_newtonf(a, d);
224}
225
226// cube root approximation using 4 iterations of Newton's method (float)
227static float newton_cbrt4f(float d)
228{
229 float a = cbrt_5f(d);
230 a = cbrta_newtonf(a, d);
231 a = cbrta_newtonf(a, d);
232 a = cbrta_newtonf(a, d);
233 return cbrta_newtonf(a, d);
234}
235
236static double TestCubeRootf(const char* szName, cuberootfnf cbrt, double rA, double rB, int rN)
237{
238 const int N = rN;
239
240 float dd = float((rB-rA) / N);
241
242 // calculate 1M numbers
243 int i=0;
244 float d = (float) rA;
245
246 double s = 0.0;
247
248 for(d=(float) rA, i=0; i<N; i++, d += dd)
249 {
250 s += cbrt(d);
251 }
252
253 double bits = 0.0;
254 double worstx=0.0;
255 double worsty=0.0;
256 int minbits=64;
257
258 for(d=(float) rA, i=0; i<N; i++, d += dd)
259 {
260 float a = cbrt((float) d);
261 float b = (float) pow((double) d, 1.0/3.0);
262
263 int bc = bits_of_precision(a, b);
264 bits += bc;
265
266 if (b > 1.0e-6)
267 {
268 if (bc < minbits)
269 {
270 minbits = bc;
271 worstx = d;
272 worsty = a;
273 }
274 }
275 }
276
277 bits /= N;
278
279 printf(" %3d mbp %6.3f abp\n", minbits, bits);
280
281 return s;
282}
283
284
285static double TestCubeRootd(const char* szName, cuberootfnd cbrt, double rA, double rB, int rN)
286{
287 const int N = rN;
288
289 double dd = (rB-rA) / N;
290
291 int i=0;
292
293 double s = 0.0;
294 double d = 0.0;
295
296 for(d=rA, i=0; i<N; i++, d += dd)
297 {
298 s += cbrt(d);
299 }
300
301
302 double bits = 0.0;
303 double worstx = 0.0;
304 double worsty = 0.0;
305 int minbits = 64;
306 for(d=rA, i=0; i<N; i++, d += dd)
307 {
308 double a = cbrt(d);
309 double b = pow(d, 1.0/3.0);
310
311 int bc = bits_of_precision(a, b); // min(53, count_matching_bitsd(a, b) - 12);
312 bits += bc;
313
314 if (b > 1.0e-6)
315 {
316 if (bc < minbits)
317 {
318 bits_of_precision(a, b);
319 minbits = bc;
320 worstx = d;
321 worsty = a;
322 }
323 }
324 }
325
326 bits /= N;
327
328 printf(" %3d mbp %6.3f abp\n", minbits, bits);
329
330 return s;
331}
332
333static int _tmain()
334{
335 // a million uniform steps through the range from 0.0 to 1.0
336 // (doing uniform steps in the log scale would be better)
337 double a = 0.0;
338 double b = 1.0;
339 int n = 1000000;
340
341 printf("32-bit float tests\n");
342 printf("----------------------------------------\n");
343 TestCubeRootf("cbrt_5f", cbrt_5f, a, b, n);
344 TestCubeRootf("pow", pow_cbrtf, a, b, n);
345 TestCubeRootf("halley x 1", halley_cbrt1f, a, b, n);
346 TestCubeRootf("halley x 2", halley_cbrt2f, a, b, n);
347 TestCubeRootf("newton x 1", newton_cbrt1f, a, b, n);
348 TestCubeRootf("newton x 2", newton_cbrt2f, a, b, n);
349 TestCubeRootf("newton x 3", newton_cbrt3f, a, b, n);
350 TestCubeRootf("newton x 4", newton_cbrt4f, a, b, n);
351 printf("\n\n");
352
353 printf("64-bit double tests\n");
354 printf("----------------------------------------\n");
355 TestCubeRootd("cbrt_5d", cbrt_5d, a, b, n);
356 TestCubeRootd("pow", pow_cbrtd, a, b, n);
357 TestCubeRootd("halley x 1", halley_cbrt1d, a, b, n);
358 TestCubeRootd("halley x 2", halley_cbrt2d, a, b, n);
359 TestCubeRootd("halley x 3", halley_cbrt3d, a, b, n);
360 TestCubeRootd("newton x 1", newton_cbrt1d, a, b, n);
361 TestCubeRootd("newton x 2", newton_cbrt2d, a, b, n);
362 TestCubeRootd("newton x 3", newton_cbrt3d, a, b, n);
363 TestCubeRootd("newton x 4", newton_cbrt4d, a, b, n);
364 printf("\n\n");
365
366 return 0;
367}
368#endif
369
370double cube_root(double x) {
371 return halley_cbrt3d(x);
372}
373
374#if TEST_ALTERNATIVES
375// http://bytes.com/topic/c/answers/754588-tips-find-cube-root-program-using-c
376/* cube root */
377int icbrt(int n) {
378 int t=0, x=(n+2)/3; /* works for n=0 and n>=1 */
379 for(; t!=x;) {
380 int x3=x*x*x;
381 t=x;
382 x*=(2*n + x3);
383 x/=(2*x3 + n);
384 }
385 return x ; /* always(?) equal to floor(n^(1/3)) */
386}
387#endif