CubeRoot.cpp revision 639df891487e40925a4f8d9a34fd3dc0c18b40a7
1// http://metamerist.com/cbrt/CubeRoot.cpp 2// 3 4#include <math.h> 5#include "CubicIntersection.h" 6 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 388