github.com/afumu/libc@v0.0.6/musl/src/math/cbrtf.c (about)

     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  }