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

     1  /* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtl.c */
     2  /*-
     3   * ====================================================
     4   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     5   * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
     6   *
     7   * Developed at SunPro, a Sun Microsystems, Inc. business.
     8   * Permission to use, copy, modify, and distribute this
     9   * software is freely granted, provided that this notice
    10   * is preserved.
    11   * ====================================================
    12   *
    13   * The argument reduction and testing for exceptional cases was
    14   * written by Steven G. Kargl with input from Bruce D. Evans
    15   * and David A. Schultz.
    16   */
    17  
    18  #include "libm.h"
    19  
    20  #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
    21  long double cbrtl(long double x)
    22  {
    23  	return cbrt(x);
    24  }
    25  #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
    26  static const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
    27  
    28  long double cbrtl(long double x)
    29  {
    30  	union ldshape u = {x}, v;
    31  	union {float f; uint32_t i;} uft;
    32  	long double r, s, t, w;
    33  	double_t dr, dt, dx;
    34  	float_t ft;
    35  	int e = u.i.se & 0x7fff;
    36  	int sign = u.i.se & 0x8000;
    37  
    38  	/*
    39  	 * If x = +-Inf, then cbrt(x) = +-Inf.
    40  	 * If x = NaN, then cbrt(x) = NaN.
    41  	 */
    42  	if (e == 0x7fff)
    43  		return x + x;
    44  	if (e == 0) {
    45  		/* Adjust subnormal numbers. */
    46  		u.f *= 0x1p120;
    47  		e = u.i.se & 0x7fff;
    48  		/* If x = +-0, then cbrt(x) = +-0. */
    49  		if (e == 0)
    50  			return x;
    51  		e -= 120;
    52  	}
    53  	e -= 0x3fff;
    54  	u.i.se = 0x3fff;
    55  	x = u.f;
    56  	switch (e % 3) {
    57  	case 1:
    58  	case -2:
    59  		x *= 2;
    60  		e--;
    61  		break;
    62  	case 2:
    63  	case -1:
    64  		x *= 4;
    65  		e -= 2;
    66  		break;
    67  	}
    68  	v.f = 1.0;
    69  	v.i.se = sign | (0x3fff + e/3);
    70  
    71  	/*
    72  	 * The following is the guts of s_cbrtf, with the handling of
    73  	 * special values removed and extra care for accuracy not taken,
    74  	 * but with most of the extra accuracy not discarded.
    75  	 */
    76  
    77  	/* ~5-bit estimate: */
    78  	uft.f = x;
    79  	uft.i = (uft.i & 0x7fffffff)/3 + B1;
    80  	ft = uft.f;
    81  
    82  	/* ~16-bit estimate: */
    83  	dx = x;
    84  	dt = ft;
    85  	dr = dt * dt * dt;
    86  	dt = dt * (dx + dx + dr) / (dx + dr + dr);
    87  
    88  	/* ~47-bit estimate: */
    89  	dr = dt * dt * dt;
    90  	dt = dt * (dx + dx + dr) / (dx + dr + dr);
    91  
    92  #if LDBL_MANT_DIG == 64
    93  	/*
    94  	 * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
    95  	 * Round it away from zero to 32 bits (32 so that t*t is exact, and
    96  	 * away from zero for technical reasons).
    97  	 */
    98  	t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32;
    99  #elif LDBL_MANT_DIG == 113
   100  	/*
   101  	 * Round dt away from zero to 47 bits.  Since we don't trust the 47,
   102  	 * add 2 47-bit ulps instead of 1 to round up.  Rounding is slow and
   103  	 * might be avoidable in this case, since on most machines dt will
   104  	 * have been evaluated in 53-bit precision and the technical reasons
   105  	 * for rounding up might not apply to either case in cbrtl() since
   106  	 * dt is much more accurate than needed.
   107  	 */
   108  	t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
   109  #endif
   110  
   111  	/*
   112  	 * Final step Newton iteration to 64 or 113 bits with
   113  	 * error < 0.667 ulps
   114  	 */
   115  	s = t*t;         /* t*t is exact */
   116  	r = x/s;         /* error <= 0.5 ulps; |r| < |t| */
   117  	w = t+t;         /* t+t is exact */
   118  	r = (r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
   119  	t = t+t*r;       /* error <= 0.5 + 0.5/3 + epsilon */
   120  
   121  	t *= v.f;
   122  	return t;
   123  }
   124  #endif