1219576Skargl/*- 2219576Skargl * ==================================================== 3219576Skargl * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4219576Skargl * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz. 5219576Skargl * 6219576Skargl * Developed at SunPro, a Sun Microsystems, Inc. business. 7219576Skargl * Permission to use, copy, modify, and distribute this 8219576Skargl * software is freely granted, provided that this notice 9219576Skargl * is preserved. 10219576Skargl * ==================================================== 11219576Skargl * 12219576Skargl * The argument reduction and testing for exceptional cases was 13219576Skargl * written by Steven G. Kargl with input from Bruce D. Evans 14219576Skargl * and David A. Schultz. 15219576Skargl */ 16219576Skargl 17219576Skargl#include <sys/cdefs.h> 18219576Skargl__FBSDID("$FreeBSD$"); 19219576Skargl 20219576Skargl#include <float.h> 21238924Skargl#ifdef __i386__ 22219576Skargl#include <ieeefp.h> 23238924Skargl#endif 24219576Skargl 25219576Skargl#include "fpmath.h" 26219576Skargl#include "math.h" 27219576Skargl#include "math_private.h" 28219576Skargl 29219576Skargl#define BIAS (LDBL_MAX_EXP - 1) 30219576Skargl 31219576Skarglstatic const unsigned 32219576Skargl B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ 33219576Skargl 34219576Skargllong double 35219576Skarglcbrtl(long double x) 36219576Skargl{ 37219576Skargl union IEEEl2bits u, v; 38219576Skargl long double r, s, t, w; 39219576Skargl double dr, dt, dx; 40219576Skargl float ft, fx; 41219576Skargl uint32_t hx; 42219576Skargl uint16_t expsign; 43219576Skargl int k; 44219576Skargl 45219576Skargl u.e = x; 46219576Skargl expsign = u.xbits.expsign; 47219576Skargl k = expsign & 0x7fff; 48219576Skargl 49219576Skargl /* 50219576Skargl * If x = +-Inf, then cbrt(x) = +-Inf. 51219576Skargl * If x = NaN, then cbrt(x) = NaN. 52219576Skargl */ 53219576Skargl if (k == BIAS + LDBL_MAX_EXP) 54219576Skargl return (x + x); 55219576Skargl 56238782Skargl ENTERI(); 57219576Skargl if (k == 0) { 58219576Skargl /* If x = +-0, then cbrt(x) = +-0. */ 59238782Skargl if ((u.bits.manh | u.bits.manl) == 0) 60238782Skargl RETURNI(x); 61219576Skargl /* Adjust subnormal numbers. */ 62219576Skargl u.e *= 0x1.0p514; 63219576Skargl k = u.bits.exp; 64219576Skargl k -= BIAS + 514; 65219576Skargl } else 66219576Skargl k -= BIAS; 67219576Skargl u.xbits.expsign = BIAS; 68219576Skargl v.e = 1; 69219576Skargl 70219576Skargl x = u.e; 71219576Skargl switch (k % 3) { 72219576Skargl case 1: 73219576Skargl case -2: 74219576Skargl x = 2*x; 75219576Skargl k--; 76219576Skargl break; 77219576Skargl case 2: 78219576Skargl case -1: 79219576Skargl x = 4*x; 80219576Skargl k -= 2; 81219576Skargl break; 82219576Skargl } 83219576Skargl v.xbits.expsign = (expsign & 0x8000) | (BIAS + k / 3); 84219576Skargl 85219576Skargl /* 86219576Skargl * The following is the guts of s_cbrtf, with the handling of 87219576Skargl * special values removed and extra care for accuracy not taken, 88219576Skargl * but with most of the extra accuracy not discarded. 89219576Skargl */ 90219576Skargl 91219576Skargl /* ~5-bit estimate: */ 92219576Skargl fx = x; 93219576Skargl GET_FLOAT_WORD(hx, fx); 94219576Skargl SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1)); 95219576Skargl 96219576Skargl /* ~16-bit estimate: */ 97219576Skargl dx = x; 98219576Skargl dt = ft; 99219576Skargl dr = dt * dt * dt; 100219576Skargl dt = dt * (dx + dx + dr) / (dx + dr + dr); 101219576Skargl 102219576Skargl /* ~47-bit estimate: */ 103219576Skargl dr = dt * dt * dt; 104219576Skargl dt = dt * (dx + dx + dr) / (dx + dr + dr); 105219576Skargl 106219576Skargl#if LDBL_MANT_DIG == 64 107219576Skargl /* 108219576Skargl * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8). 109219576Skargl * Round it away from zero to 32 bits (32 so that t*t is exact, and 110219576Skargl * away from zero for technical reasons). 111219576Skargl */ 112219576Skargl volatile double vd2 = 0x1.0p32; 113219576Skargl volatile double vd1 = 0x1.0p-31; 114219576Skargl #define vd ((long double)vd2 + vd1) 115219576Skargl 116219576Skargl t = dt + vd - 0x1.0p32; 117219576Skargl#elif LDBL_MANT_DIG == 113 118219576Skargl /* 119219576Skargl * Round dt away from zero to 47 bits. Since we don't trust the 47, 120219576Skargl * add 2 47-bit ulps instead of 1 to round up. Rounding is slow and 121219576Skargl * might be avoidable in this case, since on most machines dt will 122219576Skargl * have been evaluated in 53-bit precision and the technical reasons 123219576Skargl * for rounding up might not apply to either case in cbrtl() since 124219576Skargl * dt is much more accurate than needed. 125219576Skargl */ 126219576Skargl t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60; 127219576Skargl#else 128219576Skargl#error "Unsupported long double format" 129219576Skargl#endif 130219576Skargl 131219576Skargl /* 132219576Skargl * Final step Newton iteration to 64 or 113 bits with 133219576Skargl * error < 0.667 ulps 134219576Skargl */ 135219576Skargl s=t*t; /* t*t is exact */ 136219576Skargl r=x/s; /* error <= 0.5 ulps; |r| < |t| */ 137219576Skargl w=t+t; /* t+t is exact */ 138219576Skargl r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */ 139219576Skargl t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ 140219576Skargl 141219576Skargl t *= v.e; 142238782Skargl RETURNI(t); 143219576Skargl} 144