CubeRoot.cpp revision 639df891487e40925a4f8d9a34fd3dc0c18b40a7
1639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// http://metamerist.com/cbrt/CubeRoot.cpp 2639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// 3639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 4639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#include <math.h> 5639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#include "CubicIntersection.h" 6639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 7639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#define TEST_ALTERNATIVES 0 8639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 9639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comtypedef float (*cuberootfnf) (float); 10639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comtypedef double (*cuberootfnd) (double); 11639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 12639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// estimate bits of precision (32-bit float case) 13639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.cominline int bits_of_precision(float a, float b) 14639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 15639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double kd = 1.0 / log(2.0); 16639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 17639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (a==b) 18639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return 23; 19639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 20639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double kdmin = pow(2.0, -23.0); 21639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 22639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double d = fabs(a-b); 23639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (d < kdmin) 24639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return 23; 25639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 26639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return int(-log(d)*kd); 27639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 28639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 29639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// estiamte bits of precision (64-bit double case) 30639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.cominline int bits_of_precision(double a, double b) 31639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 32639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double kd = 1.0 / log(2.0); 33639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 34639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (a==b) 35639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return 52; 36639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 37639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double kdmin = pow(2.0, -52.0); 38639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 39639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double d = fabs(a-b); 40639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (d < kdmin) 41639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return 52; 42639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 43639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return int(-log(d)*kd); 44639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 45639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 46639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root via x^(1/3) 47639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float pow_cbrtf(float x) 48639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 49639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return (float) pow(x, 1.0f/3.0f); 50639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 51639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 52639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root via x^(1/3) 53639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double pow_cbrtd(double x) 54639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 55639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return pow(x, 1.0/3.0); 56639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 57639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 58639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using bit hack for 32-bit float 59639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float cbrt_5f(float f) 60639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 61639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com unsigned int* p = (unsigned int *) &f; 62639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com *p = *p/3 + 709921077; 63639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return f; 64639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 65639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 66639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 67639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using bit hack for 64-bit float 68639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// adapted from Kahan's cbrt 69639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double cbrt_5d(double d) 70639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 71639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const unsigned int B1 = 715094163; 72639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double t = 0.0; 73639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com unsigned int* pt = (unsigned int*) &t; 74639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com unsigned int* px = (unsigned int*) &d; 75639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com pt[1]=px[1]/3+B1; 76639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return t; 77639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 78639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 79639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 80639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using bit hack for 64-bit float 81639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// adapted from Kahan's cbrt 82639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if 0 83639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double quint_5d(double d) 84639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 85639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return sqrt(sqrt(d)); 86639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 87639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const unsigned int B1 = 71509416*5/3; 88639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double t = 0.0; 89639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com unsigned int* pt = (unsigned int*) &t; 90639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com unsigned int* px = (unsigned int*) &d; 91639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com pt[1]=px[1]/5+B1; 92639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return t; 93639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 94639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 95639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 96639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Halley's method (float) 97639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float cbrta_halleyf(const float a, const float R) 98639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 99639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const float a3 = a*a*a; 100639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const float b= a * (a3 + R + R) / (a3 + a3 + R); 101639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return b; 102639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 103639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 104639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 105639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Halley's method (double) 106639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double cbrta_halleyd(const double a, const double R) 107639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 108639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double a3 = a*a*a; 109639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double b= a * (a3 + R + R) / (a3 + a3 + R); 110639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return b; 111639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 112639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 113639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 114639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Newton's method (float) 115639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float cbrta_newtonf(const float a, const float x) 116639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 117639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// return (1.0 / 3.0) * ((a + a) + x / (a * a)); 118639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return a - (1.0f / 3.0f) * (a - x / (a*a)); 119639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 120639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 121639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Newton's method (double) 122639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double cbrta_newtond(const double a, const double x) 123639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 124639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return (1.0/3.0) * (x / (a*a) + 2*a); 125639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 126639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 127639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 1 iteration of Halley's method (double) 128639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double halley_cbrt1d(double d) 129639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 130639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 131639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_halleyd(a, d); 132639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 133639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 134639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 1 iteration of Halley's method (float) 135639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float halley_cbrt1f(float d) 136639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 137639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt_5f(d); 138639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_halleyf(a, d); 139639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 140639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 141639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Halley's method (double) 142639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double halley_cbrt2d(double d) 143639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 144639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 145639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_halleyd(a, d); 146639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_halleyd(a, d); 147639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 148639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 149639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 150639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 3 iterations of Halley's method (double) 151639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double halley_cbrt3d(double d) 152639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 153639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 154639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_halleyd(a, d); 155639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_halleyd(a, d); 156639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_halleyd(a, d); 157639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 158639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 159639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 160639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Halley's method (float) 161639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float halley_cbrt2f(float d) 162639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 163639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt_5f(d); 164639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_halleyf(a, d); 165639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_halleyf(a, d); 166639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 167639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 168639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 1 iteration of Newton's method (double) 169639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt1d(double d) 170639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 171639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 172639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtond(a, d); 173639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 174639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 175639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Newton's method (double) 176639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt2d(double d) 177639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 178639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 179639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtond(a, d); 180639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtond(a, d); 181639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 182639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 183639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 3 iterations of Newton's method (double) 184639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt3d(double d) 185639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 186639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 187639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtond(a, d); 188639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtond(a, d); 189639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtond(a, d); 190639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 191639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 192639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 4 iterations of Newton's method (double) 193639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt4d(double d) 194639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 195639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt_5d(d); 196639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtond(a, d); 197639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtond(a, d); 198639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtond(a, d); 199639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtond(a, d); 200639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 201639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 202639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Newton's method (float) 203639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt1f(float d) 204639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 205639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt_5f(d); 206639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtonf(a, d); 207639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 208639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 209639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Newton's method (float) 210639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt2f(float d) 211639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 212639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt_5f(d); 213639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtonf(a, d); 214639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtonf(a, d); 215639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 216639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 217639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 3 iterations of Newton's method (float) 218639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt3f(float d) 219639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 220639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt_5f(d); 221639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtonf(a, d); 222639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtonf(a, d); 223639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtonf(a, d); 224639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 225639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 226639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 4 iterations of Newton's method (float) 227639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt4f(float d) 228639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 229639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt_5f(d); 230639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtonf(a, d); 231639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtonf(a, d); 232639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com a = cbrta_newtonf(a, d); 233639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return cbrta_newtonf(a, d); 234639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 235639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 236639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double TestCubeRootf(const char* szName, cuberootfnf cbrt, double rA, double rB, int rN) 237639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 238639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const int N = rN; 239639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 240639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float dd = float((rB-rA) / N); 241639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 242639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com // calculate 1M numbers 243639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int i=0; 244639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float d = (float) rA; 245639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 246639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double s = 0.0; 247639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 248639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com for(d=(float) rA, i=0; i<N; i++, d += dd) 249639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 250639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com s += cbrt(d); 251639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 252639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 253639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double bits = 0.0; 254639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double worstx=0.0; 255639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double worsty=0.0; 256639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int minbits=64; 257639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 258639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com for(d=(float) rA, i=0; i<N; i++, d += dd) 259639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 260639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float a = cbrt((float) d); 261639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com float b = (float) pow((double) d, 1.0/3.0); 262639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 263639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int bc = bits_of_precision(a, b); 264639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com bits += bc; 265639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 266639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (b > 1.0e-6) 267639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 268639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (bc < minbits) 269639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 270639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com minbits = bc; 271639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com worstx = d; 272639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com worsty = a; 273639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 274639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 275639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 276639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 277639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com bits /= N; 278639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 279639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf(" %3d mbp %6.3f abp\n", minbits, bits); 280639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 281639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return s; 282639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 283639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 284639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 285639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double TestCubeRootd(const char* szName, cuberootfnd cbrt, double rA, double rB, int rN) 286639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 287639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const int N = rN; 288639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 289639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double dd = (rB-rA) / N; 290639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 291639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int i=0; 292639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 293639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double s = 0.0; 294639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double d = 0.0; 295639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 296639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com for(d=rA, i=0; i<N; i++, d += dd) 297639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 298639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com s += cbrt(d); 299639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 300639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 301639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 302639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double bits = 0.0; 303639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double worstx = 0.0; 304639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double worsty = 0.0; 305639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int minbits = 64; 306639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com for(d=rA, i=0; i<N; i++, d += dd) 307639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 308639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = cbrt(d); 309639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double b = pow(d, 1.0/3.0); 310639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 311639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int bc = bits_of_precision(a, b); // min(53, count_matching_bitsd(a, b) - 12); 312639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com bits += bc; 313639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 314639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (b > 1.0e-6) 315639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 316639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com if (bc < minbits) 317639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com { 318639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com bits_of_precision(a, b); 319639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com minbits = bc; 320639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com worstx = d; 321639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com worsty = a; 322639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 323639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 324639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 325639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 326639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com bits /= N; 327639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 328639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf(" %3d mbp %6.3f abp\n", minbits, bits); 329639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 330639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return s; 331639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 332639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 333639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic int _tmain() 334639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 335639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com // a million uniform steps through the range from 0.0 to 1.0 336639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com // (doing uniform steps in the log scale would be better) 337639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double a = 0.0; 338639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com double b = 1.0; 339639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int n = 1000000; 340639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 341639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf("32-bit float tests\n"); 342639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf("----------------------------------------\n"); 343639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("cbrt_5f", cbrt_5f, a, b, n); 344639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("pow", pow_cbrtf, a, b, n); 345639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("halley x 1", halley_cbrt1f, a, b, n); 346639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("halley x 2", halley_cbrt2f, a, b, n); 347639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("newton x 1", newton_cbrt1f, a, b, n); 348639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("newton x 2", newton_cbrt2f, a, b, n); 349639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("newton x 3", newton_cbrt3f, a, b, n); 350639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootf("newton x 4", newton_cbrt4f, a, b, n); 351639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf("\n\n"); 352639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 353639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf("64-bit double tests\n"); 354639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf("----------------------------------------\n"); 355639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("cbrt_5d", cbrt_5d, a, b, n); 356639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("pow", pow_cbrtd, a, b, n); 357639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("halley x 1", halley_cbrt1d, a, b, n); 358639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("halley x 2", halley_cbrt2d, a, b, n); 359639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("halley x 3", halley_cbrt3d, a, b, n); 360639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("newton x 1", newton_cbrt1d, a, b, n); 361639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("newton x 2", newton_cbrt2d, a, b, n); 362639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("newton x 3", newton_cbrt3d, a, b, n); 363639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com TestCubeRootd("newton x 4", newton_cbrt4d, a, b, n); 364639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf("\n\n"); 365639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 366639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return 0; 367639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 368639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 369639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 370639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comdouble cube_root(double x) { 371639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return halley_cbrt3d(x); 372639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 373639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 374639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 375639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// http://bytes.com/topic/c/answers/754588-tips-find-cube-root-program-using-c 376639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com/* cube root */ 377639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comint icbrt(int n) { 378639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int t=0, x=(n+2)/3; /* works for n=0 and n>=1 */ 379639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com for(; t!=x;) { 380639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int x3=x*x*x; 381639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com t=x; 382639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com x*=(2*n + x3); 383639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com x/=(2*x3 + n); 384639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 385639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return x ; /* always(?) equal to floor(n^(1/3)) */ 386639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 387639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 388