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