1 /*-
2  * SPDX-License-Identifier: BSD-2-Clause
4  * Copyright (c) 2009-2013 Steven G. Kargl
32  * Compute the exponential of x for Intel 80-bit format.  This is based on:
34  *   PTP Tang, "Table-driven implementation of the exponential function
35  *   in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
36  *   144-157 (1989).
55 tiny = 0x1p-10000L;
58 twom10000 = 0x1p-10000L;
61 /* log(2**16384 - 0.5) rounded towards zero: */
62 /* log(2**16384 - 0.5 + 1) rounded towards zero for expm1l() is the same: */
65 /* log(2**(-16381-64-1)) rounded towards zero: */
66 u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
73 	long double hi, lo, t, twopk;  in expl()  local
83 			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */  in expl()
84 				RETURNF(-1 / x);  in expl()
91 	} else if (ix < BIAS - 75) {	/* |x| < 0x1p-75 (includes pseudos) */  in expl()
98 	__k_expl(x, &hi, &lo, &k);  in expl()
99 	t = SUM2P(hi, lo);  in expl()
114  * Compute expm1l(x) for Intel 80-bit format.  This is based on:
116  *   PTP Tang, "Table-driven implementation of the Expm1 function
117  *   in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 18,
118  *   211-222 (1992).
132 T1 = -0.1659,				/* ~-30.625/128 * log(2) */
136  * Domain [-0.1659, 0.1659], range ~[-2.6155e-22, 2.5507e-23]:
137  * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.6
143 B3 = LD80C(0xaaaaaaaaaaaaaaab, -3,  1.66666666666666666671e-1L),
144 B4 = LD80C(0xaaaaaaaaaaaaaaac, -5,  4.16666666666666666712e-2L);
147 B5  =  8.3333333333333245e-3,		/*  0x1.111111111110cp-7 */
148 B6  =  1.3888888888888861e-3,		/*  0x1.6c16c16c16c0ap-10 */
149 B7  =  1.9841269841532042e-4,		/*  0x1.a01a01a0319f9p-13 */
150 B8  =  2.4801587302069236e-5,		/*  0x1.a01a01a03cbbcp-16 */
151 B9  =  2.7557316558468562e-6,		/*  0x1.71de37fd33d67p-19 */
152 B10 =  2.7557315829785151e-7,		/*  0x1.27e4f91418144p-22 */
153 B11 =  2.5063168199779829e-8,		/*  0x1.ae94fabdc6b27p-26 */
154 B12 =  2.0887164654459567e-9;		/*  0x1.1f122d6413fe1p-29 */
161 	long double x_lo, x2, z;  in expm1l()  local
172 			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */  in expm1l()
173 				RETURNF(-1 / x - 1);  in expm1l()
185 		if (hx & 0x8000)	/* x <= -64 */  in expm1l()
186 			RETURNF(tiny - 1);	/* good for x < -65ln2 - eps */  in expm1l()
192 		if (ix < BIAS - 74) {	/* |x| < 0x1p-74 (includes pseudos) */  in expm1l()
195 			    (0x1p100 * x + fabsl(x)) * 0x1p-100);  in expm1l()
211 		x_lo = x - x_hi;  in expm1l()
214 		if (ix >= BIAS - 7)  in expm1l()
225 	r1 = x - fn * L1;  in expm1l()
226 	r2 = fn * -L2;  in expm1l()
238 	z = r * r;  in expm1l()
239 	q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;  in expm1l()
241 	t = (long double)tbl[n2].lo + tbl[n2].hi;  in expm1l()
244 		t = SUM2P(tbl[n2].hi - 1, tbl[n2].lo * (r1 + 1) + t * q +  in expm1l()
245 		    tbl[n2].hi * r1);  in expm1l()
248 	if (k == -1) {  in expm1l()
249 		t = SUM2P(tbl[n2].hi - 2, tbl[n2].lo * (r1 + 1) + t * q +  in expm1l()
250 		    tbl[n2].hi * r1);  in expm1l()
253 	if (k < -7) {  in expm1l()
254 		t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));  in expm1l()
255 		RETURNI(t * twopk - 1);  in expm1l()
257 	if (k > 2 * LDBL_MANT_DIG - 1) {  in expm1l()
258 		t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));  in expm1l()
260 			RETURNI(t * 2 * 0x1p16383L - 1);  in expm1l()
261 		RETURNI(t * twopk - 1);  in expm1l()
264 	v.xbits.expsign = BIAS - k;  in expm1l()
267 	if (k > LDBL_MANT_DIG - 1)  in expm1l()
268 		t = SUM2P(tbl[n2].hi, tbl[n2].lo - twomk + t * (q + r1));  in expm1l()
270 		t = SUM2P(tbl[n2].hi - twomk, tbl[n2].lo + t * (q + r1));  in expm1l()