| OLD | NEW |
|---|
| (Empty) |
| 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, dou
ble 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, dou
ble 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 | |
| OLD | NEW |
|---|
Chromium Code Reviews
Issue 867213004:
remove prototype pathops code (Closed)
Base URL: https://skia.googlesource.com/skia.git@master