19e49fb63d355446b91d20ff78ad78b297e89a50dcaryclark@google.com/* 29e49fb63d355446b91d20ff78ad78b297e89a50dcaryclark@google.com * Copyright 2012 Google Inc. 39e49fb63d355446b91d20ff78ad78b297e89a50dcaryclark@google.com * 49e49fb63d355446b91d20ff78ad78b297e89a50dcaryclark@google.com * Use of this source code is governed by a BSD-style license that can be 59e49fb63d355446b91d20ff78ad78b297e89a50dcaryclark@google.com * found in the LICENSE file. 69e49fb63d355446b91d20ff78ad78b297e89a50dcaryclark@google.com */ 7d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com// http://metamerist.com/cbrt/CubeRoot.cpp 8639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// 9639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 10639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#include <math.h> 1127accef223a27fba437f5e825d99edbae20a045bcaryclark@google.com#include "CubicUtilities.h" 12639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 13639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#define TEST_ALTERNATIVES 0 14639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 15639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comtypedef float (*cuberootfnf) (float); 16639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comtypedef double (*cuberootfnd) (double); 17639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 18639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// estimate bits of precision (32-bit float case) 19639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.cominline int bits_of_precision(float a, float b) 20639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 21d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const double kd = 1.0 / log(2.0); 22639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 23d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (a==b) 24d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return 23; 25639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 26d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const double kdmin = pow(2.0, -23.0); 27639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 28d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double d = fabs(a-b); 29d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (d < kdmin) 30d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return 23; 31639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 32d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return int(-log(d)*kd); 33639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 34639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 35639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// estiamte bits of precision (64-bit double case) 36639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.cominline int bits_of_precision(double a, double b) 37639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 38d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const double kd = 1.0 / log(2.0); 39639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 40d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (a==b) 41d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return 52; 42639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 43d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const double kdmin = pow(2.0, -52.0); 44639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 45d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double d = fabs(a-b); 46d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (d < kdmin) 47d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return 52; 48639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 49d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return int(-log(d)*kd); 50639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 51639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 52639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root via x^(1/3) 53639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float pow_cbrtf(float x) 54639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 55d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return (float) pow(x, 1.0f/3.0f); 56639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 57639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 58639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root via x^(1/3) 59639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double pow_cbrtd(double x) 60639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 61d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return pow(x, 1.0/3.0); 62639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 63639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 64639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using bit hack for 32-bit float 65639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float cbrt_5f(float f) 66639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 67d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com unsigned int* p = (unsigned int *) &f; 68d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com *p = *p/3 + 709921077; 69d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return f; 70639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 71639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 72639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 73d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com// cube root approximation using bit hack for 64-bit float 74639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// adapted from Kahan's cbrt 75639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double cbrt_5d(double d) 76639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 77d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const unsigned int B1 = 715094163; 78d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double t = 0.0; 79d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com unsigned int* pt = (unsigned int*) &t; 80d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com unsigned int* px = (unsigned int*) &d; 81d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com pt[1]=px[1]/3+B1; 82d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return t; 83639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 84639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 85639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 86d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com// cube root approximation using bit hack for 64-bit float 87639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// adapted from Kahan's cbrt 88639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if 0 89639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double quint_5d(double d) 90639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 91d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return sqrt(sqrt(d)); 92d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 93d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const unsigned int B1 = 71509416*5/3; 94d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double t = 0.0; 95d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com unsigned int* pt = (unsigned int*) &t; 96d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com unsigned int* px = (unsigned int*) &d; 97d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com pt[1]=px[1]/5+B1; 98d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return t; 99639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 100639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 101639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 102639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Halley's method (float) 103639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float cbrta_halleyf(const float a, const float R) 104639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 105d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const float a3 = a*a*a; 106639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const float b= a * (a3 + R + R) / (a3 + a3 + R); 107d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return b; 108639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 109639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 110639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 111639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Halley's method (double) 112639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double cbrta_halleyd(const double a, const double R) 113639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 114d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const double a3 = a*a*a; 115639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com const double b= a * (a3 + R + R) / (a3 + a3 + R); 116d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return b; 117639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 118639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 119639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 120639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Newton's method (float) 121639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float cbrta_newtonf(const float a, const float x) 122639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 123639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// return (1.0 / 3.0) * ((a + a) + x / (a * a)); 124d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return a - (1.0f / 3.0f) * (a - x / (a*a)); 125639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 126639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 127639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// iterative cube root approximation using Newton's method (double) 128639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double cbrta_newtond(const double a, const double x) 129639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 130d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return (1.0/3.0) * (x / (a*a) + 2*a); 131639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 132639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 133639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 1 iteration of Halley's method (double) 134639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double halley_cbrt1d(double d) 135639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 136d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 137d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_halleyd(a, d); 138639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 139639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 140639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 1 iteration of Halley's method (float) 141639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float halley_cbrt1f(float d) 142639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 143d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt_5f(d); 144d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_halleyf(a, d); 145639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 146639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 147639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Halley's method (double) 148639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double halley_cbrt2d(double d) 149639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 150d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 151d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_halleyd(a, d); 152d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_halleyd(a, d); 153639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 154639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 155639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 156639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 3 iterations of Halley's method (double) 157639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double halley_cbrt3d(double d) 158639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 159d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 160d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_halleyd(a, d); 161d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_halleyd(a, d); 162d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_halleyd(a, d); 163639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 164639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 165639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 166639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Halley's method (float) 167639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float halley_cbrt2f(float d) 168639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 169d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt_5f(d); 170d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_halleyf(a, d); 171d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_halleyf(a, d); 172639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 173639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 174639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 1 iteration of Newton's method (double) 175639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt1d(double d) 176639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 177d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 178d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtond(a, d); 179639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 180639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 181639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Newton's method (double) 182639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt2d(double d) 183639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 184d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 185d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtond(a, d); 186d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtond(a, d); 187639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 188639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 189639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 3 iterations of Newton's method (double) 190639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt3d(double d) 191639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 192d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 193d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtond(a, d); 194d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtond(a, d); 195d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtond(a, d); 196639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 197639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 198639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 4 iterations of Newton's method (double) 199639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double newton_cbrt4d(double d) 200639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 201d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt_5d(d); 202d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtond(a, d); 203d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtond(a, d); 204d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtond(a, d); 205d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtond(a, d); 206639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 207639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 208639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Newton's method (float) 209639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt1f(float d) 210639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 211d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt_5f(d); 212d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtonf(a, d); 213639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 214639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 215639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 2 iterations of Newton's method (float) 216639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt2f(float d) 217639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 218d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt_5f(d); 219d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtonf(a, d); 220d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtonf(a, d); 221639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 222639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 223639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 3 iterations of Newton's method (float) 224639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt3f(float d) 225639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 226d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt_5f(d); 227d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtonf(a, d); 228d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtonf(a, d); 229d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtonf(a, d); 230639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 231639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 232639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// cube root approximation using 4 iterations of Newton's method (float) 233639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic float newton_cbrt4f(float d) 234639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 235d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt_5f(d); 236d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtonf(a, d); 237d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtonf(a, d); 238d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com a = cbrta_newtonf(a, d); 239d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return cbrta_newtonf(a, d); 240639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 241639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 242639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double TestCubeRootf(const char* szName, cuberootfnf cbrt, double rA, double rB, int rN) 243639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 244d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const int N = rN; 245d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 246d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float dd = float((rB-rA) / N); 247d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 248d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com // calculate 1M numbers 249d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int i=0; 250d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float d = (float) rA; 251d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 252d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double s = 0.0; 253d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 254d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com for(d=(float) rA, i=0; i<N; i++, d += dd) 255d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 256d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com s += cbrt(d); 257d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 258d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 259d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double bits = 0.0; 260d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double worstx=0.0; 261d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double worsty=0.0; 262d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int minbits=64; 263d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 264d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com for(d=(float) rA, i=0; i<N; i++, d += dd) 265d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 266d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float a = cbrt((float) d); 267d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com float b = (float) pow((double) d, 1.0/3.0); 268d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 269d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int bc = bits_of_precision(a, b); 270d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com bits += bc; 271d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 272d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (b > 1.0e-6) 273d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 274d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (bc < minbits) 275d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 276d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com minbits = bc; 277d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com worstx = d; 278d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com worsty = a; 279d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 280d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 281d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 282d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 283d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com bits /= N; 284639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 285639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf(" %3d mbp %6.3f abp\n", minbits, bits); 286639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 287d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return s; 288639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 289639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 290639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 291639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic double TestCubeRootd(const char* szName, cuberootfnd cbrt, double rA, double rB, int rN) 292639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 293d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com const int N = rN; 294d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 295d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double dd = (rB-rA) / N; 296d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 297d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int i=0; 298d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 299d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double s = 0.0; 300d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double d = 0.0; 301d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 302d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com for(d=rA, i=0; i<N; i++, d += dd) 303d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 304d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com s += cbrt(d); 305d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 306d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 307d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 308d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double bits = 0.0; 309d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double worstx = 0.0; 310d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double worsty = 0.0; 311d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int minbits = 64; 312d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com for(d=rA, i=0; i<N; i++, d += dd) 313d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 314d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = cbrt(d); 315d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double b = pow(d, 1.0/3.0); 316d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 317d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int bc = bits_of_precision(a, b); // min(53, count_matching_bitsd(a, b) - 12); 318d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com bits += bc; 319d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 320d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (b > 1.0e-6) 321d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 322d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com if (bc < minbits) 323d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com { 324d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com bits_of_precision(a, b); 325d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com minbits = bc; 326d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com worstx = d; 327d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com worsty = a; 328d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 329d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 330d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com } 331d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 332d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com bits /= N; 333639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 334639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com printf(" %3d mbp %6.3f abp\n", minbits, bits); 335639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 336d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return s; 337639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 338639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 339639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comstatic int _tmain() 340639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com{ 341d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com // a million uniform steps through the range from 0.0 to 1.0 342d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com // (doing uniform steps in the log scale would be better) 343d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double a = 0.0; 344d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com double b = 1.0; 345d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com int n = 1000000; 346d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 347d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com printf("32-bit float tests\n"); 348d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com printf("----------------------------------------\n"); 349d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("cbrt_5f", cbrt_5f, a, b, n); 350d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("pow", pow_cbrtf, a, b, n); 351d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("halley x 1", halley_cbrt1f, a, b, n); 352d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("halley x 2", halley_cbrt2f, a, b, n); 353d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("newton x 1", newton_cbrt1f, a, b, n); 354d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("newton x 2", newton_cbrt2f, a, b, n); 355d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("newton x 3", newton_cbrt3f, a, b, n); 356d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootf("newton x 4", newton_cbrt4f, a, b, n); 357d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com printf("\n\n"); 358d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 359d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com printf("64-bit double tests\n"); 360d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com printf("----------------------------------------\n"); 361d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("cbrt_5d", cbrt_5d, a, b, n); 362d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("pow", pow_cbrtd, a, b, n); 363d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("halley x 1", halley_cbrt1d, a, b, n); 364d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("halley x 2", halley_cbrt2d, a, b, n); 365d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("halley x 3", halley_cbrt3d, a, b, n); 366d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("newton x 1", newton_cbrt1d, a, b, n); 367d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("newton x 2", newton_cbrt2d, a, b, n); 368d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("newton x 3", newton_cbrt3d, a, b, n); 369d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com TestCubeRootd("newton x 4", newton_cbrt4d, a, b, n); 370d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com printf("\n\n"); 371d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com 372d6176b0dcacb124539e0cfd051e6d93a9782f020rmistry@google.com return 0; 373639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 374639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 375639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 376639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comdouble cube_root(double x) { 3779f60291c5375457f8adf228dbe6e8ff1186b13e1caryclark@google.com if (approximately_zero_cubed(x)) { 378235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com return 0; 379235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com } 380235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com double result = halley_cbrt3d(fabs(x)); 381235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com if (x < 0) { 382235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com result = -result; 383235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com } 384235f56a92f6eb6accbb243e11b3c45e3798f38f2caryclark@google.com return result; 385639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 386639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com 387639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#if TEST_ALTERNATIVES 388639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com// http://bytes.com/topic/c/answers/754588-tips-find-cube-root-program-using-c 389639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com/* cube root */ 390639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.comint icbrt(int n) { 391639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int t=0, x=(n+2)/3; /* works for n=0 and n>=1 */ 392639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com for(; t!=x;) { 393639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com int x3=x*x*x; 394639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com t=x; 395639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com x*=(2*n + x3); 396639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com x/=(2*x3 + n); 397639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com } 398639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com return x ; /* always(?) equal to floor(n^(1/3)) */ 399639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com} 400639df891487e40925a4f8d9a34fd3dc0c18b40a7caryclark@google.com#endif 401