s_expl.c revision 330897
1228072Sbapt/*- 2228072Sbapt * SPDX-License-Identifier: BSD-2-Clause-FreeBSD 3228072Sbapt * 4228072Sbapt * Copyright (c) 2009-2013 Steven G. Kargl 5228072Sbapt * All rights reserved. 6228072Sbapt * 7228072Sbapt * Redistribution and use in source and binary forms, with or without 8228072Sbapt * modification, are permitted provided that the following conditions 9228072Sbapt * are met: 10228072Sbapt * 1. Redistributions of source code must retain the above copyright 11228072Sbapt * notice unmodified, this list of conditions, and the following 12228072Sbapt * disclaimer. 13228072Sbapt * 2. Redistributions in binary form must reproduce the above copyright 14228072Sbapt * notice, this list of conditions and the following disclaimer in the 15228072Sbapt * documentation and/or other materials provided with the distribution. 16228072Sbapt * 17228072Sbapt * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 18228072Sbapt * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 19228072Sbapt * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 20228072Sbapt * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 21228072Sbapt * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 22228072Sbapt * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 23228072Sbapt * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 24228072Sbapt * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 25228072Sbapt * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 26228072Sbapt * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 27228072Sbapt * 28228072Sbapt * Optimized by Bruce D. Evans. 29228072Sbapt */ 30228072Sbapt 31228072Sbapt#include <sys/cdefs.h> 32228072Sbapt__FBSDID("$FreeBSD: stable/11/lib/msun/ld80/s_expl.c 330897 2018-03-14 03:19:51Z eadler $"); 33228072Sbapt 34228072Sbapt/** 35228072Sbapt * Compute the exponential of x for Intel 80-bit format. This is based on: 36228072Sbapt * 37228072Sbapt * PTP Tang, "Table-driven implementation of the exponential function 38228072Sbapt * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 39228072Sbapt * 144-157 (1989). 40228072Sbapt * 41228072Sbapt * where the 32 table entries have been expanded to INTERVALS (see below). 42228072Sbapt */ 43228072Sbapt 44228072Sbapt#include <float.h> 45228072Sbapt 46228072Sbapt#ifdef __i386__ 47228072Sbapt#include <ieeefp.h> 48228072Sbapt#endif 49228072Sbapt 50228072Sbapt#include "fpmath.h" 51228072Sbapt#include "math.h" 52228072Sbapt#include "math_private.h" 53228072Sbapt#include "k_expl.h" 54228072Sbapt 55228072Sbapt/* XXX Prevent compilers from erroneously constant folding these: */ 56228072Sbaptstatic const volatile long double 57228072Sbapthuge = 0x1p10000L, 58228072Sbapttiny = 0x1p-10000L; 59228072Sbapt 60228072Sbaptstatic const long double 61228072Sbapttwom10000 = 0x1p-10000L; 62228072Sbapt 63228072Sbaptstatic const union IEEEl2bits 64228072Sbapt/* log(2**16384 - 0.5) rounded towards zero: */ 65228072Sbapt/* log(2**16384 - 0.5 + 1) rounded towards zero for expm1l() is the same: */ 66228072Sbapto_thresholdu = LD80C(0xb17217f7d1cf79ab, 13, 11356.5234062941439488L), 67228072Sbapt#define o_threshold (o_thresholdu.e) 68228072Sbapt/* log(2**(-16381-64-1)) rounded towards zero: */ 69228072Sbaptu_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L); 70228072Sbapt#define u_threshold (u_thresholdu.e) 71228072Sbapt 72228072Sbaptlong double 73228072Sbaptexpl(long double x) 74228072Sbapt{ 75228072Sbapt union IEEEl2bits u; 76228072Sbapt long double hi, lo, t, twopk; 77228072Sbapt int k; 78228072Sbapt uint16_t hx, ix; 79228072Sbapt 80228072Sbapt DOPRINT_START(&x); 81228072Sbapt 82228072Sbapt /* Filter out exceptional cases. */ 83228072Sbapt u.e = x; 84228072Sbapt hx = u.xbits.expsign; 85228072Sbapt ix = hx & 0x7fff; 86228072Sbapt if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 87228072Sbapt if (ix == BIAS + LDBL_MAX_EXP) { 88228072Sbapt if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */ 89228072Sbapt RETURNP(-1 / x); 90228072Sbapt RETURNP(x + x); /* x is +Inf, +NaN or unsupported */ 91228072Sbapt } 92228072Sbapt if (x > o_threshold) 93228072Sbapt RETURNP(huge * huge); 94228072Sbapt if (x < u_threshold) 95228072Sbapt RETURNP(tiny * tiny); 96228072Sbapt } else if (ix < BIAS - 75) { /* |x| < 0x1p-75 (includes pseudos) */ 97228072Sbapt RETURN2P(1, x); /* 1 with inexact iff x != 0 */ 98228072Sbapt } 99228072Sbapt 100228072Sbapt ENTERI(); 101228072Sbapt 102228072Sbapt twopk = 1; 103228072Sbapt __k_expl(x, &hi, &lo, &k); 104228072Sbapt t = SUM2P(hi, lo); 105228072Sbapt 106228072Sbapt /* Scale by 2**k. */ 107228072Sbapt if (k >= LDBL_MIN_EXP) { 108228072Sbapt if (k == LDBL_MAX_EXP) 109228072Sbapt RETURNI(t * 2 * 0x1p16383L); 110228072Sbapt SET_LDBL_EXPSIGN(twopk, BIAS + k); 111228072Sbapt RETURNI(t * twopk); 112228072Sbapt } else { 113228072Sbapt SET_LDBL_EXPSIGN(twopk, BIAS + k + 10000); 114228072Sbapt RETURNI(t * twopk * twom10000); 115228072Sbapt } 116228072Sbapt} 117228072Sbapt 118228072Sbapt/** 119228072Sbapt * Compute expm1l(x) for Intel 80-bit format. This is based on: 120228072Sbapt * 121228072Sbapt * PTP Tang, "Table-driven implementation of the Expm1 function 122228072Sbapt * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 18, 123228072Sbapt * 211-222 (1992). 124228072Sbapt */ 125228072Sbapt 126228072Sbapt/* 127228072Sbapt * Our T1 and T2 are chosen to be approximately the points where method 128228072Sbapt * A and method B have the same accuracy. Tang's T1 and T2 are the 129228072Sbapt * points where method A's accuracy changes by a full bit. For Tang, 130228072Sbapt * this drop in accuracy makes method A immediately less accurate than 131228072Sbapt * method B, but our larger INTERVALS makes method A 2 bits more 132228072Sbapt * accurate so it remains the most accurate method significantly 133228072Sbapt * closer to the origin despite losing the full bit in our extended 134228072Sbapt * range for it. 135228072Sbapt */ 136228072Sbaptstatic const double 137228072SbaptT1 = -0.1659, /* ~-30.625/128 * log(2) */ 138228072SbaptT2 = 0.1659; /* ~30.625/128 * log(2) */ 139228072Sbapt 140228072Sbapt/* 141228072Sbapt * Domain [-0.1659, 0.1659], range ~[-2.6155e-22, 2.5507e-23]: 142228072Sbapt * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.6 143228072Sbapt * 144228072Sbapt * XXX the coeffs aren't very carefully rounded, and I get 2.8 more bits, 145228072Sbapt * but unlike for ld128 we can't drop any terms. 146228072Sbapt */ 147228072Sbaptstatic const union IEEEl2bits 148228072SbaptB3 = LD80C(0xaaaaaaaaaaaaaaab, -3, 1.66666666666666666671e-1L), 149228072SbaptB4 = LD80C(0xaaaaaaaaaaaaaaac, -5, 4.16666666666666666712e-2L); 150228072Sbapt 151228072Sbaptstatic const double 152228072SbaptB5 = 8.3333333333333245e-3, /* 0x1.111111111110cp-7 */ 153228072SbaptB6 = 1.3888888888888861e-3, /* 0x1.6c16c16c16c0ap-10 */ 154228072SbaptB7 = 1.9841269841532042e-4, /* 0x1.a01a01a0319f9p-13 */ 155228072SbaptB8 = 2.4801587302069236e-5, /* 0x1.a01a01a03cbbcp-16 */ 156228072SbaptB9 = 2.7557316558468562e-6, /* 0x1.71de37fd33d67p-19 */ 157228072SbaptB10 = 2.7557315829785151e-7, /* 0x1.27e4f91418144p-22 */ 158228072SbaptB11 = 2.5063168199779829e-8, /* 0x1.ae94fabdc6b27p-26 */ 159228072SbaptB12 = 2.0887164654459567e-9; /* 0x1.1f122d6413fe1p-29 */ 160228072Sbapt 161228072Sbaptlong double 162228072Sbaptexpm1l(long double x) 163228072Sbapt{ 164228072Sbapt union IEEEl2bits u, v; 165228072Sbapt long double fn, hx2_hi, hx2_lo, q, r, r1, r2, t, twomk, twopk, x_hi; 166228072Sbapt long double x_lo, x2, z; 167228072Sbapt long double x4; 168228072Sbapt int k, n, n2; 169228072Sbapt uint16_t hx, ix; 170228072Sbapt 171228072Sbapt DOPRINT_START(&x); 172228072Sbapt 173228072Sbapt /* Filter out exceptional cases. */ 174228072Sbapt u.e = x; 175228072Sbapt hx = u.xbits.expsign; 176228072Sbapt ix = hx & 0x7fff; 177228072Sbapt if (ix >= BIAS + 6) { /* |x| >= 64 or x is NaN */ 178228072Sbapt if (ix == BIAS + LDBL_MAX_EXP) { 179228072Sbapt if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */ 180228072Sbapt RETURNP(-1 / x - 1); 181228072Sbapt RETURNP(x + x); /* x is +Inf, +NaN or unsupported */ 182228072Sbapt } 183228072Sbapt if (x > o_threshold) 184228072Sbapt RETURNP(huge * huge); 185228072Sbapt /* 186228072Sbapt * expm1l() never underflows, but it must avoid 187228072Sbapt * unrepresentable large negative exponents. We used a 188228072Sbapt * much smaller threshold for large |x| above than in 189228072Sbapt * expl() so as to handle not so large negative exponents 190228072Sbapt * in the same way as large ones here. 191228072Sbapt */ 192228072Sbapt if (hx & 0x8000) /* x <= -64 */ 193228072Sbapt RETURN2P(tiny, -1); /* good for x < -65ln2 - eps */ 194228072Sbapt } 195228072Sbapt 196228072Sbapt ENTERI(); 197228072Sbapt 198250874Sjkim if (T1 < x && x < T2) { 199250875Sjkim if (ix < BIAS - 74) { /* |x| < 0x1p-74 (includes pseudos) */ 200250875Sjkim /* x (rounded) with inexact if x != 0: */ 201250875Sjkim RETURNPI(x == 0 ? x : 202250874Sjkim (0x1p100 * x + fabsl(x)) * 0x1p-100); 203250875Sjkim } 204250874Sjkim 205250874Sjkim x2 = x * x; 206250874Sjkim x4 = x2 * x2; 207250874Sjkim q = x4 * (x2 * (x4 * 208228072Sbapt /* 209228072Sbapt * XXX the number of terms is no longer good for 210228072Sbapt * pairwise grouping of all except B3, and the 211228072Sbapt * grouping is no longer from highest down. 212228072Sbapt */ 213228072Sbapt (x2 * B12 + (x * B11 + B10)) + 214228072Sbapt (x2 * (x * B9 + B8) + (x * B7 + B6))) + 215228072Sbapt (x * B5 + B4.e)) + x2 * x * B3.e; 216228072Sbapt 217228072Sbapt x_hi = (float)x; 218228072Sbapt x_lo = x - x_hi; 219228072Sbapt hx2_hi = x_hi * x_hi / 2; 220228072Sbapt hx2_lo = x_lo * (x + x_hi) / 2; 221228072Sbapt if (ix >= BIAS - 7) 222228072Sbapt RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q); 223228072Sbapt else 224228072Sbapt RETURN2PI(x, hx2_lo + q + hx2_hi); 225228072Sbapt } 226228072Sbapt 227228072Sbapt /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */ 228228072Sbapt /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 229228072Sbapt fn = x * INV_L + 0x1.8p63 - 0x1.8p63; 230228072Sbapt#if defined(HAVE_EFFICIENT_IRINTL) 231250125Sjkim n = irintl(fn); 232228072Sbapt#elif defined(HAVE_EFFICIENT_IRINT) 233228072Sbapt n = irint(fn); 234228072Sbapt#else 235228072Sbapt n = (int)fn; 236228072Sbapt#endif 237228072Sbapt n2 = (unsigned)n % INTERVALS; 238228072Sbapt k = n >> LOG2_INTERVALS; 239228072Sbapt r1 = x - fn * L1; 240228072Sbapt r2 = fn * -L2; 241228072Sbapt r = r1 + r2; 242228072Sbapt 243228072Sbapt /* Prepare scale factor. */ 244228072Sbapt v.e = 1; 245228072Sbapt v.xbits.expsign = BIAS + k; 246228072Sbapt twopk = v.e; 247228072Sbapt 248228072Sbapt /* 249228072Sbapt * Evaluate lower terms of 250228072Sbapt * expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). 251228072Sbapt */ 252228072Sbapt z = r * r; 253228072Sbapt q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6; 254228072Sbapt 255228072Sbapt t = (long double)tbl[n2].lo + tbl[n2].hi; 256228072Sbapt 257228072Sbapt if (k == 0) { 258228072Sbapt t = SUM2P(tbl[n2].hi - 1, tbl[n2].lo * (r1 + 1) + t * q + 259228072Sbapt tbl[n2].hi * r1); 260228072Sbapt RETURNI(t); 261228072Sbapt } 262228072Sbapt if (k == -1) { 263228072Sbapt t = SUM2P(tbl[n2].hi - 2, tbl[n2].lo * (r1 + 1) + t * q + 264228072Sbapt tbl[n2].hi * r1); 265228072Sbapt RETURNI(t / 2); 266228072Sbapt } 267228072Sbapt if (k < -7) { 268228072Sbapt t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1)); 269228072Sbapt RETURNI(t * twopk - 1); 270228072Sbapt } 271228072Sbapt if (k > 2 * LDBL_MANT_DIG - 1) { 272228072Sbapt t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1)); 273228072Sbapt if (k == LDBL_MAX_EXP) 274228072Sbapt RETURNI(t * 2 * 0x1p16383L - 1); 275228072Sbapt RETURNI(t * twopk - 1); 276228072Sbapt } 277228072Sbapt 278228072Sbapt v.xbits.expsign = BIAS - k; 279228072Sbapt twomk = v.e; 280228072Sbapt 281228072Sbapt if (k > LDBL_MANT_DIG - 1) 282228072Sbapt t = SUM2P(tbl[n2].hi, tbl[n2].lo - twomk + t * (q + r1)); 283228072Sbapt else 284228072Sbapt t = SUM2P(tbl[n2].hi - twomk, tbl[n2].lo + t * (q + r1)); 285228072Sbapt RETURNI(t * twopk); 286228072Sbapt} 287228072Sbapt