Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	lib/libm/src/ld80/s_exp2l.c	Lines:	0	42	0.0 %
Date:	2017-11-07	Branches:	0	22	0.0 %


/*	$OpenBSD: s_exp2l.c,v 1.3 2014/07/21 01:51:11 guenther Exp $	*/
/*-
 * Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG>
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#include <sys/types.h>
#include <machine/ieee.h>
#include <float.h>
#include <math.h>
#include <stdint.h>

#define	TBLBITS	7
#define	TBLSIZE	(1 << TBLBITS)

#define	BIAS	(LDBL_MAX_EXP - 1)
#define	EXPMASK	(BIAS + LDBL_MAX_EXP)

static const long double huge = 0x1p10000L;
#if 0 /* XXX Prevent gcc from erroneously constant folding this. */
static const long double twom10000 = 0x1p-10000L;
#else
static volatile long double twom10000 = 0x1p-10000L;
#endif

static const double
    redux     = 0x1.8p63 / TBLSIZE,
    P1        = 0x1.62e42fefa39efp-1,
    P2        = 0x1.ebfbdff82c58fp-3,
    P3        = 0x1.c6b08d7049fap-5,
    P4        = 0x1.3b2ab6fba4da5p-7,
    P5        = 0x1.5d8804780a736p-10,
    P6        = 0x1.430918835e33dp-13;

static const double tbl[TBLSIZE * 2] = {
	0x1.6a09e667f3bcdp-1,	-0x1.bdd3413b2648p-55,
	0x1.6c012750bdabfp-1,	-0x1.2895667ff0cp-57,
	0x1.6dfb23c651a2fp-1,	-0x1.bbe3a683c88p-58,
	0x1.6ff7df9519484p-1,	-0x1.83c0f25860fp-56,
	0x1.71f75e8ec5f74p-1,	-0x1.16e4786887bp-56,
	0x1.73f9a48a58174p-1,	-0x1.0a8d96c65d5p-55,
	0x1.75feb564267c9p-1,	-0x1.0245957316ep-55,
	0x1.780694fde5d3fp-1,	 0x1.866b80a0216p-55,
	0x1.7a11473eb0187p-1,	-0x1.41577ee0499p-56,
	0x1.7c1ed0130c132p-1,	 0x1.f124cd1164ep-55,
	0x1.7e2f336cf4e62p-1,	 0x1.05d02ba157ap-57,
	0x1.80427543e1a12p-1,	-0x1.27c86626d97p-55,
	0x1.82589994cce13p-1,	-0x1.d4c1dd41533p-55,
	0x1.8471a4623c7adp-1,	-0x1.8d684a341cep-56,
	0x1.868d99b4492edp-1,	-0x1.fc6f89bd4f68p-55,
	0x1.88ac7d98a6699p-1,	 0x1.994c2f37cb5p-55,
	0x1.8ace5422aa0dbp-1,	 0x1.6e9f156864bp-55,
	0x1.8cf3216b5448cp-1,	-0x1.0d55e32e9e4p-57,
	0x1.8f1ae99157736p-1,	 0x1.5cc13a2e397p-56,
	0x1.9145b0b91ffc6p-1,	-0x1.dd6792e5825p-55,
	0x1.93737b0cdc5e5p-1,	-0x1.75fc781b58p-58,
	0x1.95a44cbc8520fp-1,	-0x1.64b7c96a5fp-57,
	0x1.97d829fde4e5p-1,	-0x1.d185b7c1b86p-55,
	0x1.9a0f170ca07bap-1,	-0x1.173bd91cee6p-55,
	0x1.9c49182a3f09p-1,	 0x1.c7c46b071f2p-57,
	0x1.9e86319e32323p-1,	 0x1.824ca78e64cp-57,
	0x1.a0c667b5de565p-1,	-0x1.359495d1cd5p-55,
	0x1.a309bec4a2d33p-1,	 0x1.6305c7ddc368p-55,
	0x1.a5503b23e255dp-1,	-0x1.d2f6edb8d42p-55,
	0x1.a799e1330b358p-1,	 0x1.bcb7ecac564p-55,
	0x1.a9e6b5579fdbfp-1,	 0x1.0fac90ef7fdp-55,
	0x1.ac36bbfd3f37ap-1,	-0x1.f9234cae76dp-56,
	0x1.ae89f995ad3adp-1,	 0x1.7a1cd345dcc8p-55,
	0x1.b0e07298db666p-1,	-0x1.bdef54c80e4p-55,
	0x1.b33a2b84f15fbp-1,	-0x1.2805e3084d8p-58,
	0x1.b59728de5593ap-1,	-0x1.c71dfbbba6ep-55,
	0x1.b7f76f2fb5e47p-1,	-0x1.5584f7e54acp-57,
	0x1.ba5b030a1064ap-1,	-0x1.efcd30e5429p-55,
	0x1.bcc1e904bc1d2p-1,	 0x1.23dd07a2d9fp-56,
	0x1.bf2c25bd71e09p-1,	-0x1.efdca3f6b9c8p-55,
	0x1.c199bdd85529cp-1,	 0x1.11065895049p-56,
	0x1.c40ab5fffd07ap-1,	 0x1.b4537e083c6p-55,
	0x1.c67f12e57d14bp-1,	 0x1.2884dff483c8p-55,
	0x1.c8f6d9406e7b5p-1,	 0x1.1acbc48805cp-57,
	0x1.cb720dcef9069p-1,	 0x1.503cbd1e94ap-57,
	0x1.cdf0b555dc3fap-1,	-0x1.dd83b53829dp-56,
	0x1.d072d4a07897cp-1,	-0x1.cbc3743797a8p-55,
	0x1.d2f87080d89f2p-1,	-0x1.d487b719d858p-55,
	0x1.d5818dcfba487p-1,	 0x1.2ed02d75b37p-56,
	0x1.d80e316c98398p-1,	-0x1.11ec18bedep-55,
	0x1.da9e603db3285p-1,	 0x1.c2300696db5p-55,
	0x1.dd321f301b46p-1,	 0x1.2da5778f019p-55,
	0x1.dfc97337b9b5fp-1,	-0x1.1a5cd4f184b8p-55,
	0x1.e264614f5a129p-1,	-0x1.7b627817a148p-55,
	0x1.e502ee78b3ff6p-1,	 0x1.39e8980a9cdp-56,
	0x1.e7a51fbc74c83p-1,	 0x1.2d522ca0c8ep-55,
	0x1.ea4afa2a490dap-1,	-0x1.e9c23179c288p-55,
	0x1.ecf482d8e67f1p-1,	-0x1.c93f3b411ad8p-55,
	0x1.efa1bee615a27p-1,	 0x1.dc7f486a4b68p-55,
	0x1.f252b376bba97p-1,	 0x1.3a1a5bf0d8e8p-55,
	0x1.f50765b6e454p-1,	 0x1.9d3e12dd8a18p-55,
	0x1.f7bfdad9cbe14p-1,	-0x1.dbb12d00635p-55,
	0x1.fa7c1819e90d8p-1,	 0x1.74853f3a593p-56,
	0x1.fd3c22b8f71f1p-1,	 0x1.2eb74966578p-58,
	0x1p+0,	 0x0p+0,
	0x1.0163da9fb3335p+0,	 0x1.b61299ab8cd8p-54,
	0x1.02c9a3e778061p+0,	-0x1.19083535b08p-56,
	0x1.04315e86e7f85p+0,	-0x1.0a31c1977c98p-54,
	0x1.059b0d3158574p+0,	 0x1.d73e2a475b4p-55,
	0x1.0706b29ddf6dep+0,	-0x1.c91dfe2b13cp-55,
	0x1.0874518759bc8p+0,	 0x1.186be4bb284p-57,
	0x1.09e3ecac6f383p+0,	 0x1.14878183161p-54,
	0x1.0b5586cf9890fp+0,	 0x1.8a62e4adc61p-54,
	0x1.0cc922b7247f7p+0,	 0x1.01edc16e24f8p-54,
	0x1.0e3ec32d3d1a2p+0,	 0x1.03a1727c58p-59,
	0x1.0fb66affed31bp+0,	-0x1.b9bedc44ebcp-57,
	0x1.11301d0125b51p+0,	-0x1.6c51039449bp-54,
	0x1.12abdc06c31ccp+0,	-0x1.1b514b36ca8p-58,
	0x1.1429aaea92dep+0,	-0x1.32fbf9af1368p-54,
	0x1.15a98c8a58e51p+0,	 0x1.2406ab9eeabp-55,
	0x1.172b83c7d517bp+0,	-0x1.19041b9d78ap-55,
	0x1.18af9388c8deap+0,	-0x1.11023d1970f8p-54,
	0x1.1a35beb6fcb75p+0,	 0x1.e5b4c7b4969p-55,
	0x1.1bbe084045cd4p+0,	-0x1.95386352ef6p-54,
	0x1.1d4873168b9aap+0,	 0x1.e016e00a264p-54,
	0x1.1ed5022fcd91dp+0,	-0x1.1df98027bb78p-54,
	0x1.2063b88628cd6p+0,	 0x1.dc775814a85p-55,
	0x1.21f49917ddc96p+0,	 0x1.2a97e9494a6p-55,
	0x1.2387a6e756238p+0,	 0x1.9b07eb6c7058p-54,
	0x1.251ce4fb2a63fp+0,	 0x1.ac155bef4f5p-55,
	0x1.26b4565e27cddp+0,	 0x1.2bd339940eap-55,
	0x1.284dfe1f56381p+0,	-0x1.a4c3a8c3f0d8p-54,
	0x1.29e9df51fdee1p+0,	 0x1.612e8afad12p-55,
	0x1.2b87fd0dad99p+0,	-0x1.10adcd6382p-59,
	0x1.2d285a6e4030bp+0,	 0x1.0024754db42p-54,
	0x1.2ecafa93e2f56p+0,	 0x1.1ca0f45d524p-56,
	0x1.306fe0a31b715p+0,	 0x1.6f46ad23183p-55,
	0x1.32170fc4cd831p+0,	 0x1.a9ce78e1804p-55,
	0x1.33c08b26416ffp+0,	 0x1.327218436598p-54,
	0x1.356c55f929ff1p+0,	-0x1.b5cee5c4e46p-55,
	0x1.371a7373aa9cbp+0,	-0x1.63aeabf42ebp-54,
	0x1.38cae6d05d866p+0,	-0x1.e958d3c99048p-54,
	0x1.3a7db34e59ff7p+0,	-0x1.5e436d661f6p-56,
	0x1.3c32dc313a8e5p+0,	-0x1.efff8375d2ap-54,
	0x1.3dea64c123422p+0,	 0x1.ada0911f09fp-55,
	0x1.3fa4504ac801cp+0,	-0x1.7d023f956fap-54,
	0x1.4160a21f72e2ap+0,	-0x1.ef3691c309p-58,
	0x1.431f5d950a897p+0,	-0x1.1c7dde35f7ap-55,
	0x1.44e086061892dp+0,	 0x1.89b7a04ef8p-59,
	0x1.46a41ed1d0057p+0,	 0x1.c944bd1648a8p-54,
	0x1.486a2b5c13cdp+0,	 0x1.3c1a3b69062p-56,
	0x1.4a32af0d7d3dep+0,	 0x1.9cb62f3d1be8p-54,
	0x1.4bfdad5362a27p+0,	 0x1.d4397afec42p-56,
	0x1.4dcb299fddd0dp+0,	 0x1.8ecdbbc6a78p-54,
	0x1.4f9b2769d2ca7p+0,	-0x1.4b309d25958p-54,
	0x1.516daa2cf6642p+0,	-0x1.f768569bd94p-55,
	0x1.5342b569d4f82p+0,	-0x1.07abe1db13dp-55,
	0x1.551a4ca5d920fp+0,	-0x1.d689cefede6p-55,
	0x1.56f4736b527dap+0,	 0x1.9bb2c011d938p-54,
	0x1.58d12d497c7fdp+0,	 0x1.295e15b9a1ep-55,
	0x1.5ab07dd485429p+0,	 0x1.6324c0546478p-54,
	0x1.5c9268a5946b7p+0,	 0x1.c4b1b81698p-60,
	0x1.5e76f15ad2148p+0,	 0x1.ba6f93080e68p-54,
	0x1.605e1b976dc09p+0,	-0x1.3e2429b56de8p-54,
	0x1.6247eb03a5585p+0,	-0x1.383c17e40b48p-54,
	0x1.6434634ccc32p+0,	-0x1.c483c759d89p-55,
	0x1.6623882552225p+0,	-0x1.bb60987591cp-54,
	0x1.68155d44ca973p+0,	 0x1.038ae44f74p-57,
};

/*
 * exp2l(x): compute the base 2 exponential of x
 *
 * Accuracy: Peak error < 0.511 ulp.
 *
 * Method: (equally-spaced tables)
 *
 *   Reduce x:
 *     x = 2**k + y, for integer k and |y| <= 1/2.
 *     Thus we have exp2l(x) = 2**k * exp2(y).
 *
 *   Reduce y:
 *     y = i/TBLSIZE + z for integer i near y * TBLSIZE.
 *     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
 *     with |z| <= 2**-(TBLBITS+1).
 *
 *   We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
 *   degree-6 minimax polynomial with maximum error under 2**-69.
 *   The table entries each have 104 bits of accuracy, encoded as
 *   a pair of double precision values.
 */
long double
exp2l(long double x)
{
	union {
		long double e;
		struct ieee_ext bits;
	} u, v;
	long double r, twopk, twopkp10000, z, t_hi, t_lo;
	uint32_t es, hx, ix, i0;
	int k;

	/* Filter out exceptional cases. */
	u.e = x;
	hx = (u.bits.ext_sign << 15) | u.bits.ext_exp;
	ix = hx & EXPMASK;
	if (ix >= BIAS + 14) {		/* |x| >= 16384 or x is NaN */
		if (ix == BIAS + LDBL_MAX_EXP) {
			if ((u.bits.ext_frach != 1U << 31 &&
			    u.bits.ext_fracl != 0) || (hx & 0x8000) == 0)
				return (x + x);	/* x is +Inf or NaN */
			else
				return (0.0);	/* x is -Inf */
		}
		if (x >= 16384)
			return (huge * huge);	/* overflow */
		if (x <= -16446)
			return (twom10000 * twom10000);	/* underflow */
	} else if (ix <= BIAS - 66) {		/* |x| < 0x1p-66 */
		return (1.0 + x);
	}

	/*
	 * Reduce x, computing z, i0, and k. The low bits of x + redux
	 * contain the 16-bit integer part of the exponent (k) followed by
	 * TBLBITS fractional bits (i0). We use bit tricks to extract these
	 * as integers, then set z to the remainder.
	 *
	 * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
	 * Then the low-order word of x + redux is 0x000abc12,
	 * We split this into k = 0xabc and i0 = 0x12 (adjusted to
	 * index into the table), then we compute z = 0x0.003456p0.
	 *
	 * XXX If the exponent is negative, the computation of k depends on
	 *     '>>' doing sign extension.
	 */
	u.e = x + redux;
	i0 = u.bits.ext_fracl + TBLSIZE / 2;
	k = (int)i0 >> TBLBITS;
	i0 = (i0 & (TBLSIZE - 1)) << 1;
	u.e -= redux;
	z = x - u.e;

	v.bits.ext_frach = 1U << 31;
	v.bits.ext_fracl = 0;

	if (k >= LDBL_MIN_EXP) {
		es = LDBL_MAX_EXP - 1 + k;
		v.bits.ext_exp = es & 0x7fffffff;
		v.bits.ext_sign = u.bits.ext_sign >> 15;
		twopk = v.e;
	} else {
		es = LDBL_MAX_EXP - 1 + k + 10000;
		v.bits.ext_exp = es & 0x7fffffff;
		v.bits.ext_sign = u.bits.ext_sign >> 15;
		twopkp10000 = v.e;
	}

	/* Compute r = exp2l(y) = exp2lt[i0] * p(z). */
	t_hi = tbl[i0];
	t_lo = tbl[i0 + 1];
	/* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */
	r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4
	    + z * (P5 + z * P6))))) + t_hi;

	/* Scale by 2**k. */
	if (k >= LDBL_MIN_EXP) {
		if (k == LDBL_MAX_EXP)
			return (r * 2.0 * 0x1p16383L);
		return (r * twopk);
	} else {
		return (r * twopkp10000 * twom10000);
	}
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/* $OpenBSD: s_exp2l.c,v 1.3 2014/07/21 01:51:11 guenther Exp $ */
2			/*-
3			* Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG>
4			* All rights reserved.
5			*
6			* Redistribution and use in source and binary forms, with or without
7			* modification, are permitted provided that the following conditions
8			* are met:
9			* 1. Redistributions of source code must retain the above copyright
10			* notice, this list of conditions and the following disclaimer.
11			* 2. Redistributions in binary form must reproduce the above copyright
12			* notice, this list of conditions and the following disclaimer in the
13			* documentation and/or other materials provided with the distribution.
14			*
15			* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
16			* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18			* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19			* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20			* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21			* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22			* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23			* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24			* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25			* SUCH DAMAGE.
26			*/
27
28			#include <sys/types.h>
29			#include <machine/ieee.h>
30			#include <float.h>
31			#include <math.h>
32			#include <stdint.h>
33
34			#define TBLBITS 7
35			#define TBLSIZE (1 << TBLBITS)
36
37			#define BIAS (LDBL_MAX_EXP - 1)
38			#define EXPMASK (BIAS + LDBL_MAX_EXP)
39
40			static const long double huge = 0x1p10000L;
41			#if 0 /* XXX Prevent gcc from erroneously constant folding this. */
42			static const long double twom10000 = 0x1p-10000L;
43			#else
44			static volatile long double twom10000 = 0x1p-10000L;
45			#endif
46
47			static const double
48			redux = 0x1.8p63 / TBLSIZE,
49			P1 = 0x1.62e42fefa39efp-1,
50			P2 = 0x1.ebfbdff82c58fp-3,
51			P3 = 0x1.c6b08d7049fap-5,
52			P4 = 0x1.3b2ab6fba4da5p-7,
53			P5 = 0x1.5d8804780a736p-10,
54			P6 = 0x1.430918835e33dp-13;
55
56			static const double tbl[TBLSIZE * 2] = {
57			0x1.6a09e667f3bcdp-1, -0x1.bdd3413b2648p-55,
58			0x1.6c012750bdabfp-1, -0x1.2895667ff0cp-57,
59			0x1.6dfb23c651a2fp-1, -0x1.bbe3a683c88p-58,
60			0x1.6ff7df9519484p-1, -0x1.83c0f25860fp-56,
61			0x1.71f75e8ec5f74p-1, -0x1.16e4786887bp-56,
62			0x1.73f9a48a58174p-1, -0x1.0a8d96c65d5p-55,
63			0x1.75feb564267c9p-1, -0x1.0245957316ep-55,
64			0x1.780694fde5d3fp-1, 0x1.866b80a0216p-55,
65			0x1.7a11473eb0187p-1, -0x1.41577ee0499p-56,
66			0x1.7c1ed0130c132p-1, 0x1.f124cd1164ep-55,
67			0x1.7e2f336cf4e62p-1, 0x1.05d02ba157ap-57,
68			0x1.80427543e1a12p-1, -0x1.27c86626d97p-55,
69			0x1.82589994cce13p-1, -0x1.d4c1dd41533p-55,
70			0x1.8471a4623c7adp-1, -0x1.8d684a341cep-56,
71			0x1.868d99b4492edp-1, -0x1.fc6f89bd4f68p-55,
72			0x1.88ac7d98a6699p-1, 0x1.994c2f37cb5p-55,
73			0x1.8ace5422aa0dbp-1, 0x1.6e9f156864bp-55,
74			0x1.8cf3216b5448cp-1, -0x1.0d55e32e9e4p-57,
75			0x1.8f1ae99157736p-1, 0x1.5cc13a2e397p-56,
76			0x1.9145b0b91ffc6p-1, -0x1.dd6792e5825p-55,
77			0x1.93737b0cdc5e5p-1, -0x1.75fc781b58p-58,
78			0x1.95a44cbc8520fp-1, -0x1.64b7c96a5fp-57,
79			0x1.97d829fde4e5p-1, -0x1.d185b7c1b86p-55,
80			0x1.9a0f170ca07bap-1, -0x1.173bd91cee6p-55,
81			0x1.9c49182a3f09p-1, 0x1.c7c46b071f2p-57,
82			0x1.9e86319e32323p-1, 0x1.824ca78e64cp-57,
83			0x1.a0c667b5de565p-1, -0x1.359495d1cd5p-55,
84			0x1.a309bec4a2d33p-1, 0x1.6305c7ddc368p-55,
85			0x1.a5503b23e255dp-1, -0x1.d2f6edb8d42p-55,
86			0x1.a799e1330b358p-1, 0x1.bcb7ecac564p-55,
87			0x1.a9e6b5579fdbfp-1, 0x1.0fac90ef7fdp-55,
88			0x1.ac36bbfd3f37ap-1, -0x1.f9234cae76dp-56,
89			0x1.ae89f995ad3adp-1, 0x1.7a1cd345dcc8p-55,
90			0x1.b0e07298db666p-1, -0x1.bdef54c80e4p-55,
91			0x1.b33a2b84f15fbp-1, -0x1.2805e3084d8p-58,
92			0x1.b59728de5593ap-1, -0x1.c71dfbbba6ep-55,
93			0x1.b7f76f2fb5e47p-1, -0x1.5584f7e54acp-57,
94			0x1.ba5b030a1064ap-1, -0x1.efcd30e5429p-55,
95			0x1.bcc1e904bc1d2p-1, 0x1.23dd07a2d9fp-56,
96			0x1.bf2c25bd71e09p-1, -0x1.efdca3f6b9c8p-55,
97			0x1.c199bdd85529cp-1, 0x1.11065895049p-56,
98			0x1.c40ab5fffd07ap-1, 0x1.b4537e083c6p-55,
99			0x1.c67f12e57d14bp-1, 0x1.2884dff483c8p-55,
100			0x1.c8f6d9406e7b5p-1, 0x1.1acbc48805cp-57,
101			0x1.cb720dcef9069p-1, 0x1.503cbd1e94ap-57,
102			0x1.cdf0b555dc3fap-1, -0x1.dd83b53829dp-56,
103			0x1.d072d4a07897cp-1, -0x1.cbc3743797a8p-55,
104			0x1.d2f87080d89f2p-1, -0x1.d487b719d858p-55,
105			0x1.d5818dcfba487p-1, 0x1.2ed02d75b37p-56,
106			0x1.d80e316c98398p-1, -0x1.11ec18bedep-55,
107			0x1.da9e603db3285p-1, 0x1.c2300696db5p-55,
108			0x1.dd321f301b46p-1, 0x1.2da5778f019p-55,
109			0x1.dfc97337b9b5fp-1, -0x1.1a5cd4f184b8p-55,
110			0x1.e264614f5a129p-1, -0x1.7b627817a148p-55,
111			0x1.e502ee78b3ff6p-1, 0x1.39e8980a9cdp-56,
112			0x1.e7a51fbc74c83p-1, 0x1.2d522ca0c8ep-55,
113			0x1.ea4afa2a490dap-1, -0x1.e9c23179c288p-55,
114			0x1.ecf482d8e67f1p-1, -0x1.c93f3b411ad8p-55,
115			0x1.efa1bee615a27p-1, 0x1.dc7f486a4b68p-55,
116			0x1.f252b376bba97p-1, 0x1.3a1a5bf0d8e8p-55,
117			0x1.f50765b6e454p-1, 0x1.9d3e12dd8a18p-55,
118			0x1.f7bfdad9cbe14p-1, -0x1.dbb12d00635p-55,
119			0x1.fa7c1819e90d8p-1, 0x1.74853f3a593p-56,
120			0x1.fd3c22b8f71f1p-1, 0x1.2eb74966578p-58,
121			0x1p+0, 0x0p+0,
122			0x1.0163da9fb3335p+0, 0x1.b61299ab8cd8p-54,
123			0x1.02c9a3e778061p+0, -0x1.19083535b08p-56,
124			0x1.04315e86e7f85p+0, -0x1.0a31c1977c98p-54,
125			0x1.059b0d3158574p+0, 0x1.d73e2a475b4p-55,
126			0x1.0706b29ddf6dep+0, -0x1.c91dfe2b13cp-55,
127			0x1.0874518759bc8p+0, 0x1.186be4bb284p-57,
128			0x1.09e3ecac6f383p+0, 0x1.14878183161p-54,
129			0x1.0b5586cf9890fp+0, 0x1.8a62e4adc61p-54,
130			0x1.0cc922b7247f7p+0, 0x1.01edc16e24f8p-54,
131			0x1.0e3ec32d3d1a2p+0, 0x1.03a1727c58p-59,
132			0x1.0fb66affed31bp+0, -0x1.b9bedc44ebcp-57,
133			0x1.11301d0125b51p+0, -0x1.6c51039449bp-54,
134			0x1.12abdc06c31ccp+0, -0x1.1b514b36ca8p-58,
135			0x1.1429aaea92dep+0, -0x1.32fbf9af1368p-54,
136			0x1.15a98c8a58e51p+0, 0x1.2406ab9eeabp-55,
137			0x1.172b83c7d517bp+0, -0x1.19041b9d78ap-55,
138			0x1.18af9388c8deap+0, -0x1.11023d1970f8p-54,
139			0x1.1a35beb6fcb75p+0, 0x1.e5b4c7b4969p-55,
140			0x1.1bbe084045cd4p+0, -0x1.95386352ef6p-54,
141			0x1.1d4873168b9aap+0, 0x1.e016e00a264p-54,
142			0x1.1ed5022fcd91dp+0, -0x1.1df98027bb78p-54,
143			0x1.2063b88628cd6p+0, 0x1.dc775814a85p-55,
144			0x1.21f49917ddc96p+0, 0x1.2a97e9494a6p-55,
145			0x1.2387a6e756238p+0, 0x1.9b07eb6c7058p-54,
146			0x1.251ce4fb2a63fp+0, 0x1.ac155bef4f5p-55,
147			0x1.26b4565e27cddp+0, 0x1.2bd339940eap-55,
148			0x1.284dfe1f56381p+0, -0x1.a4c3a8c3f0d8p-54,
149			0x1.29e9df51fdee1p+0, 0x1.612e8afad12p-55,
150			0x1.2b87fd0dad99p+0, -0x1.10adcd6382p-59,
151			0x1.2d285a6e4030bp+0, 0x1.0024754db42p-54,
152			0x1.2ecafa93e2f56p+0, 0x1.1ca0f45d524p-56,
153			0x1.306fe0a31b715p+0, 0x1.6f46ad23183p-55,
154			0x1.32170fc4cd831p+0, 0x1.a9ce78e1804p-55,
155			0x1.33c08b26416ffp+0, 0x1.327218436598p-54,
156			0x1.356c55f929ff1p+0, -0x1.b5cee5c4e46p-55,
157			0x1.371a7373aa9cbp+0, -0x1.63aeabf42ebp-54,
158			0x1.38cae6d05d866p+0, -0x1.e958d3c99048p-54,
159			0x1.3a7db34e59ff7p+0, -0x1.5e436d661f6p-56,
160			0x1.3c32dc313a8e5p+0, -0x1.efff8375d2ap-54,
161			0x1.3dea64c123422p+0, 0x1.ada0911f09fp-55,
162			0x1.3fa4504ac801cp+0, -0x1.7d023f956fap-54,
163			0x1.4160a21f72e2ap+0, -0x1.ef3691c309p-58,
164			0x1.431f5d950a897p+0, -0x1.1c7dde35f7ap-55,
165			0x1.44e086061892dp+0, 0x1.89b7a04ef8p-59,
166			0x1.46a41ed1d0057p+0, 0x1.c944bd1648a8p-54,
167			0x1.486a2b5c13cdp+0, 0x1.3c1a3b69062p-56,
168			0x1.4a32af0d7d3dep+0, 0x1.9cb62f3d1be8p-54,
169			0x1.4bfdad5362a27p+0, 0x1.d4397afec42p-56,
170			0x1.4dcb299fddd0dp+0, 0x1.8ecdbbc6a78p-54,
171			0x1.4f9b2769d2ca7p+0, -0x1.4b309d25958p-54,
172			0x1.516daa2cf6642p+0, -0x1.f768569bd94p-55,
173			0x1.5342b569d4f82p+0, -0x1.07abe1db13dp-55,
174			0x1.551a4ca5d920fp+0, -0x1.d689cefede6p-55,
175			0x1.56f4736b527dap+0, 0x1.9bb2c011d938p-54,
176			0x1.58d12d497c7fdp+0, 0x1.295e15b9a1ep-55,
177			0x1.5ab07dd485429p+0, 0x1.6324c0546478p-54,
178			0x1.5c9268a5946b7p+0, 0x1.c4b1b81698p-60,
179			0x1.5e76f15ad2148p+0, 0x1.ba6f93080e68p-54,
180			0x1.605e1b976dc09p+0, -0x1.3e2429b56de8p-54,
181			0x1.6247eb03a5585p+0, -0x1.383c17e40b48p-54,
182			0x1.6434634ccc32p+0, -0x1.c483c759d89p-55,
183			0x1.6623882552225p+0, -0x1.bb60987591cp-54,
184			0x1.68155d44ca973p+0, 0x1.038ae44f74p-57,
185			};
186
187			/*
188			* exp2l(x): compute the base 2 exponential of x
189			*
190			* Accuracy: Peak error < 0.511 ulp.
191			*
192			* Method: (equally-spaced tables)
193			*
194			* Reduce x:
195			* x = 2**k + y, for integer k and \|y\| <= 1/2.
196			* Thus we have exp2l(x) = 2*k exp2(y).
197			*
198			* Reduce y:
199			* y = i/TBLSIZE + z for integer i near y * TBLSIZE.
200			* Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
201			* with \|z\| <= 2**-(TBLBITS+1).
202			*
203			* We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
204			* degree-6 minimax polynomial with maximum error under 2**-69.
205			* The table entries each have 104 bits of accuracy, encoded as
206			* a pair of double precision values.
207			*/
208			long double
209			exp2l(long double x)
210			{
211			union {
212			long double e;
213			struct ieee_ext bits;
214			} u, v;
215			long double r, twopk, twopkp10000, z, t_hi, t_lo;
216			uint32_t es, hx, ix, i0;
217			int k;
218
219			/* Filter out exceptional cases. */
220			u.e = x;
221			hx = (u.bits.ext_sign << 15) \| u.bits.ext_exp;
222			ix = hx & EXPMASK;
223			if (ix >= BIAS + 14) { /* \|x\| >= 16384 or x is NaN */
224			if (ix == BIAS + LDBL_MAX_EXP) {
225			if ((u.bits.ext_frach != 1U << 31 &&
226			u.bits.ext_fracl != 0) \|\| (hx & 0x8000) == 0)
227			return (x + x); /* x is +Inf or NaN */
228			else
229			return (0.0); /* x is -Inf */
230			}
231			if (x >= 16384)
232			return (huge * huge); /* overflow */
233			if (x <= -16446)
234			return (twom10000 * twom10000); /* underflow */
235			} else if (ix <= BIAS - 66) { /* \|x\| < 0x1p-66 */
236			return (1.0 + x);
237			}
238
239			/*
240			* Reduce x, computing z, i0, and k. The low bits of x + redux
241			* contain the 16-bit integer part of the exponent (k) followed by
242			* TBLBITS fractional bits (i0). We use bit tricks to extract these
243			* as integers, then set z to the remainder.
244			*
245			* Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
246			* Then the low-order word of x + redux is 0x000abc12,
247			* We split this into k = 0xabc and i0 = 0x12 (adjusted to
248			* index into the table), then we compute z = 0x0.003456p0.
249			*
250			* XXX If the exponent is negative, the computation of k depends on
251			* '>>' doing sign extension.
252			*/
253			u.e = x + redux;
254			i0 = u.bits.ext_fracl + TBLSIZE / 2;
255			k = (int)i0 >> TBLBITS;
256			i0 = (i0 & (TBLSIZE - 1)) << 1;
257			u.e -= redux;
258			z = x - u.e;
259
260			v.bits.ext_frach = 1U << 31;
261			v.bits.ext_fracl = 0;
262
263			if (k >= LDBL_MIN_EXP) {
264			es = LDBL_MAX_EXP - 1 + k;
265			v.bits.ext_exp = es & 0x7fffffff;
266			v.bits.ext_sign = u.bits.ext_sign >> 15;
267			twopk = v.e;
268			} else {
269			es = LDBL_MAX_EXP - 1 + k + 10000;
270			v.bits.ext_exp = es & 0x7fffffff;
271			v.bits.ext_sign = u.bits.ext_sign >> 15;
272			twopkp10000 = v.e;
273			}
274
275			/* Compute r = exp2l(y) = exp2lt[i0] * p(z). */
276			t_hi = tbl[i0];
277			t_lo = tbl[i0 + 1];
278			/* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */
279			r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4
280			+ z * (P5 + z * P6))))) + t_hi;
281
282			/* Scale by 2*k. /
283			if (k >= LDBL_MIN_EXP) {
284			if (k == LDBL_MAX_EXP)
285			return (r * 2.0 * 0x1p16383L);
286			return (r * twopk);
287			} else {
288			return (r * twopkp10000 * twom10000);
289			}
290			}