lib/msun/src/s_cbrtf.c

   1 /* s_cbrtf.c -- float version of s_cbrt.c.
   2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   3  * Debugged and optimized by Bruce D. Evans.
   4  */
   5
   6 /*
   7  * ====================================================
   8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9  *
  10  * Developed at SunPro, a Sun Microsystems, Inc. business.
  11  * Permission to use, copy, modify, and distribute this
  12  * software is freely granted, provided that this notice
  13  * is preserved.
  14  * ====================================================
  15  */
  16
  17 #ifndef lint
  18 static char rcsid[] = "$FreeBSD$";
  19 #endif
  20
  21 #include "math.h"
  22 #include "math_private.h"
  23
  24 /* cbrtf(x)
  25  * Return cube root of x
  26  */
  27 static const unsigned
  28         B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
  29         B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
  30
  31 /* |1/cbrt(x) - p(x)| < 2**-14.5 (~[-4.37e-4, 4.366e-5]). */
  32 static const float
  33 P0 =  1.5586718321,                     /*  0x3fc7828f */
  34 P1 = -0.78271341324,                    /* -0xbf485fe8 */
  35 P2 =  0.22403796017;                    /*  0x3e656a35 */
  36
  37 float
  38 cbrtf(float x)
  39 {
  40         float r,s,t,w;
  41         int32_t hx;
  42         u_int32_t sign;
  43         u_int32_t high;
  44
  45         GET_FLOAT_WORD(hx,x);
  46         sign=hx&0x80000000;             /* sign= sign(x) */
  47         hx  ^=sign;
  48         if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */
  49         if(hx==0)
  50             return(x);                  /* cbrt(0) is itself */
  51
  52     /* rough cbrt to 5 bits */
  53         if(hx<0x00800000) {             /* subnormal number */
  54             SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
  55             t*=x;
  56             GET_FLOAT_WORD(high,t);
  57             SET_FLOAT_WORD(t,sign|((high&0x7fffffff)/3+B2));
  58         } else
  59             SET_FLOAT_WORD(t,sign|(hx/3+B1));
  60
  61     /* new cbrt to 14 bits */
  62         r=(t*t)*(t/x);
  63         t=t*((P0+r*P1)+(r*r)*P2);
  64
  65     /*
  66      * Round t away from zero to 12 bits (sloppily except for ensuring that
  67      * the result is larger in magnitude than cbrt(x) but not much more than
  68      * 1 12-bit ulp larger).  With rounding towards zero, the error bound
  69      * would be ~5/6 instead of ~4/6, and with t 2 12-bit ulps larger the
  70      * infinite-precision error in the Newton approximation would affect
  71      * the second digit instead of the third digit of 4/6 = 0.666..., etc.
  72      */
  73         GET_FLOAT_WORD(high,t);
  74         SET_FLOAT_WORD(t,(high+0x1800)&0xfffff000);
  75
  76     /* one step Newton iteration to 24 bits with error < 0.669 ulps */
  77         s=t*t;                          /* t*t is exact */
  78         r=x/s;                          /* error <= 0.5 ulps; |r| < |t| */
  79         w=t+t;                          /* t+t is exact */
  80         r=(r-t)/(w+r);                  /* r-t is exact; w+r ~= 3*t */
  81         t=t+t*r;                        /* error <= 0.5 + 0.5/3 + epsilon */
  82
  83         return(t);
  84 }