src/math/cbrtf.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */
   2 /*
   3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   4  * Debugged and optimized by Bruce D. Evans.
   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 /* cbrtf(x)
  17  * Return cube root of x
  18  */
  19
  20 #include <math.h>
  21 #include <stdint.h>
  22
  23 static const unsigned
  24 B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
  25 B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
  26
  27 float cbrtf(float x)
  28 {
  29         double_t r,T;
  30         union {float f; uint32_t i;} u = {x};
  31         uint32_t hx = u.i & 0x7fffffff;
  32
  33         if (hx >= 0x7f800000)  /* cbrt(NaN,INF) is itself */
  34                 return x + x;
  35
  36         /* rough cbrt to 5 bits */
  37         if (hx < 0x00800000) {  /* zero or subnormal? */
  38                 if (hx == 0)
  39                         return x;  /* cbrt(+-0) is itself */
  40                 u.f = x*0x1p24f;
  41                 hx = u.i & 0x7fffffff;
  42                 hx = hx/3 + B2;
  43         } else
  44                 hx = hx/3 + B1;
  45         u.i &= 0x80000000;
  46         u.i |= hx;
  47
  48         /*
  49          * First step Newton iteration (solving t*t-x/t == 0) to 16 bits.  In
  50          * double precision so that its terms can be arranged for efficiency
  51          * without causing overflow or underflow.
  52          */
  53         T = u.f;
  54         r = T*T*T;
  55         T = T*((double_t)x+x+r)/(x+r+r);
  56
  57         /*
  58          * Second step Newton iteration to 47 bits.  In double precision for
  59          * efficiency and accuracy.
  60          */
  61         r = T*T*T;
  62         T = T*((double_t)x+x+r)/(x+r+r);
  63
  64         /* rounding to 24 bits is perfect in round-to-nearest mode */
  65         return T;
  66 }