Cube root function. More...

#include "prim.h"

Include dependency graph for cbrt.c:

Functions
double	cbrt (double x)
	Cube root function. More...

Variables
static const unsigned	B1 = 715094163

static const unsigned	B2 = 696219795

static const double	C = 5.42857142857142815906e-01

static const double	D = -7.05306122448979611050e-01

static const double	E = 1.41428571428571436819e+00

static const double	F = 1.60714285714285720630e+00

static const double	G = 3.57142857142857150787e-01

Detailed Description

Cube root function.

Definition in file cbrt.c.

Function Documentation

◆ cbrt()

double cbrt ( double x )

Cube root function.

Definition at line 62 of file cbrt.c.

References B1, B2, C, D, E, F, and G.

Referenced by RealNumber::CubeRoot().

 {
     int32_t hx, lx, ht;
     double r, s, t = 0.0, w;
     uint32_t sign;
 
     GET_HIGH_WORD(hx, x);
 
     // sign = sign(x)
     sign = hx & 0x80000000;
     hx ^= sign;
 
     // cbrt(NaN,INF) is itself
     if (hx >= 0x7FF00000)
     {
         return x + x;
     }
 
     GET_LOW_WORD(lx, x);
 
     // cbrt(0) is itself
     if ((hx | lx) == 0)
     {
         return x;
     }
 
     // x <- |x|
     SET_HIGH_WORD(x, hx);
 
     // rough cbrt to 5 bits
     if (hx < 0x00100000) // subnormal number
     {
         SET_HIGH_WORD(t, 0x43500000); // set t= 2**54
         t *= x;
         GET_HIGH_WORD(ht, t);
         SET_HIGH_WORD(t, ht / 3 + B2);
     }
     else
     {
         SET_HIGH_WORD(t, hx / 3 + B1);
     }
 
     // new cbrt to 23 bits, may be implemented in single precision
     r = t * t / x;
     s = C + r * t;
     t *= G + F / (s + E + D / s);
 
     // chopped to 20 bits and make it larger than cbrt(x)
     SET_LOW_WORD(t, 0);
     GET_HIGH_WORD(ht, t);
     SET_HIGH_WORD(t, ht + 0x00000001);
 
     // one step newton iteration to 53 bits with error less than 0.667 ulps
     s = t * t; // t*t is exact
     r = x / s;
     w = t + t;
     r = (r - t) / (w + r); // r-s is exact
     t = t + t * r;
 
     // retore the sign bit
     GET_HIGH_WORD(ht, t);
     SET_HIGH_WORD(t, ht | sign);
 
     return (t);
 }