compiler/stdc/math/s_cbrtl.c

   1 /* s_cbrtf.c -- float version of s_cbrt.c.
   2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   3  * Debugged and optimized by Bruce D. Evans.
   4  */
   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  * The argument reduction and testing for exceptional cases was
  17  * written by Steven G. Kargl with input from Bruce D. Evans
  18  * and David A. Schultz.
  19  */
  20
  21 #ifndef lint
  22 static char rcsid[] = "$FreeBSD: src/lib/msun/src/s_cbrtl.c,v 1.1 2011/03/12 19:37:35 kargl Exp $";
  23 #endif
  24
  25 #include <float.h>
  26 #include "math.h"
  27
  28 #include "fpmath.h"
  29 #include "math_private.h"
  30 #if defined(__x86_64__)
  31 #include "x86_64/ieeefp.h"
  32 #elif defined(__i386__)
  33 #include "i386/ieeefp.h"
  34 #endif
  35
  36
  37 #define BIAS    (LDBL_MAX_EXP - 1)
  38
  39 static const unsigned
  40     B1 = 709958130;     /* B1 = (127-127.0/3-0.03306235651)*2**23 */
  41
  42 long double
  43 cbrtl(long double x)
  44 {
  45         union IEEEl2bits u, v;
  46         long double r, s, t, w;
  47         double dr, dt, dx;
  48         float ft, fx;
  49         uint32_t hx;
  50         uint16_t expsign;
  51         int k;
  52
  53         u.e = x;
  54         expsign = u.xbits.expsign;
  55         k = expsign & 0x7fff;
  56
  57         /*
  58          * If x = +-Inf, then cbrt(x) = +-Inf.
  59          * If x = NaN, then cbrt(x) = NaN.
  60          */
  61         if (k == BIAS + LDBL_MAX_EXP)
  62                 return (x + x);
  63
  64         ENTERI();
  65
  66         if (k == 0) {
  67                 /* If x = +-0, then cbrt(x) = +-0. */
  68                 if ((u.bits.manh | u.bits.manl) == 0)
  69                         RETURNI(x);
  70                 /* Adjust subnormal numbers. */
  71                 u.e *= 0x1.0p514;
  72                 k = u.bits.exp;
  73                 k -= BIAS + 514;
  74         } else
  75                 k -= BIAS;
  76         u.xbits.expsign = BIAS;
  77         v.e = 1;
  78
  79         x = u.e;
  80         switch (k % 3) {
  81         case 1:
  82         case -2:
  83                 x = 2*x;
  84                 k--;
  85                 break;
  86         case 2:
  87         case -1:
  88                 x = 4*x;
  89                 k -= 2;
  90                 break;
  91         }
  92         v.xbits.expsign = (expsign & 0x8000) | (BIAS + k / 3);
  93
  94         /*
  95          * The following is the guts of s_cbrtf, with the handling of
  96          * special values removed and extra care for accuracy not taken,
  97          * but with most of the extra accuracy not discarded.
  98          */
  99
 100         /* ~5-bit estimate: */
 101         fx = x;
 102         GET_FLOAT_WORD(hx, fx);
 103         SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
 104
 105         /* ~16-bit estimate: */
 106         dx = x;
 107         dt = ft;
 108         dr = dt * dt * dt;
 109         dt = dt * (dx + dx + dr) / (dx + dr + dr);
 110
 111         /* ~47-bit estimate: */
 112         dr = dt * dt * dt;
 113         dt = dt * (dx + dx + dr) / (dx + dr + dr);
 114
 115 #if LDBL_MANT_DIG == 64
 116         /*
 117          * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
 118          * Round it away from zero to 32 bits (32 so that t*t is exact, and
 119          * away from zero for technical reasons).
 120          */
 121         volatile double vd2 = 0x1.0p32;
 122         volatile double vd1 = 0x1.0p-31;
 123         #define vd ((long double)vd2 + vd1)
 124
 125         t = dt + vd - 0x1.0p32;
 126 #elif LDBL_MANT_DIG == 113
 127         /*
 128          * Round dt away from zero to 47 bits.  Since we don't trust the 47,
 129          * add 2 47-bit ulps instead of 1 to round up.  Rounding is slow and
 130          * might be avoidable in this case, since on most machines dt will
 131          * have been evaluated in 53-bit precision and the technical reasons
 132          * for rounding up might not apply to either case in cbrtl() since
 133          * dt is much more accurate than needed.
 134          */
 135         t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
 136 #else
 137 #error "Unsupported long double format"
 138 #endif
 139
 140         /*
 141          * Final step Newton iteration to 64 or 113 bits with
 142          * error < 0.667 ulps
 143          */
 144         s=t*t;                          /* t*t is exact */
 145         r=x/s;                          /* error <= 0.5 ulps; |r| < |t| */
 146         w=t+t;                          /* t+t is exact */
 147         r=(r-t)/(w+r);                  /* r-t is exact; w+r ~= 3*t */
 148         t=t+t*r;                        /* error <= 0.5 + 0.5/3 + epsilon */
 149
 150         t *= v.e;
 151         RETURNI(t);
 152 }