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