cbrtf.c raw

   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  }
  67