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 "libm.h"
  21
  22 static const unsigned
  23 B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
  24 B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
  25
  26 float cbrtf(float x)
  27 {
  28         double r,T;
  29         float t;
  30         int32_t hx;
  31         uint32_t sign;
  32         uint32_t high;
  33
  34         GET_FLOAT_WORD(hx, x);
  35         sign = hx & 0x80000000;
  36         hx ^= sign;
  37         if (hx >= 0x7f800000)  /* cbrt(NaN,INF) is itself */
  38                 return x + x;
  39
  40         /* rough cbrt to 5 bits */
  41         if (hx < 0x00800000) {  /* zero or subnormal? */
  42                 if (hx == 0)
  43                         return x;  /* cbrt(+-0) is itself */
  44                 SET_FLOAT_WORD(t, 0x4b800000);  /* set t = 2**24 */
  45                 t *= x;
  46                 GET_FLOAT_WORD(high, t);
  47                 SET_FLOAT_WORD(t, sign|((high&0x7fffffff)/3+B2));
  48         } else
  49                 SET_FLOAT_WORD(t, sign|(hx/3+B1));
  50
  51         /*
  52          * First step Newton iteration (solving t*t-x/t == 0) to 16 bits.  In
  53          * double precision so that its terms can be arranged for efficiency
  54          * without causing overflow or underflow.
  55          */
  56         T = t;
  57         r = T*T*T;
  58         T = T*((double)x+x+r)/(x+r+r);
  59
  60         /*
  61          * Second step Newton iteration to 47 bits.  In double precision for
  62          * efficiency and accuracy.
  63          */
  64         r = T*T*T;
  65         T = T*((double)x+x+r)/(x+r+r);
  66
  67         /* rounding to 24 bits is perfect in round-to-nearest mode */
  68         return T;
  69 }