standards/152415: [libm] implementation of expl()
Steven G. Kargl
kargl at troutmask.apl.washington.edu
Fri Nov 19 23:10:09 UTC 2010
>Number: 152415
>Category: standards
>Synopsis: [libm] implementation of expl()
>Confidential: no
>Severity: non-critical
>Priority: low
>Responsible: freebsd-standards
>State: open
>Quarter:
>Keywords:
>Date-Required:
>Class: change-request
>Submitter-Id: current-users
>Arrival-Date: Fri Nov 19 23:10:08 UTC 2010
>Closed-Date:
>Last-Modified:
>Originator: Steven G. Kargl
>Release: FreeBSD 9.0-CURRENT amd64
>Organization:
apl/uw
>Environment:
System: FreeBSD troutmask.apl.washington.edu 9.0-CURRENT FreeBSD 9.0-CURRENT #6 r213691M: Tue Oct 26 11:01:30 PDT 2010 kargl at troutmask.apl.washington.edu:/usr/obj/usr/src/sys/SPEW amd64
>Description:
FreeBSD currently lacks a long double version of the exponential
function, also known as expl(). An algorithm for both float
and double versions of this function, expf() and exp(), has been
given by Tang in
PTP Tang, "Table-driven implementation of the exponential function
in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
144-157 (1989).
I have implemented Tang's algorithm for expf(), exp(), and expl().
My testing on amd64 has found
Tang Libm
Range Max ULP sec Max ULP sec # of calls
expf: [-18.0:88.72] 0.525668 1.11 0.901775 1.33 23124111
exp: [-37.0:709.7] 0.526008 1.47 0.895724 1.60 23000001
expl: [-44.0:11356] 0.506063 2.81 23000001
where "Range" is the range over which the function will return
a finite normal result. The 'sec' column is the self-seconds
reported by gprof for the "# of calls" to the function. Thus,
for over 23M function calls expf() and exp() based on Tang
are more accurate and actually faster than the corresponding
functions in libm.
In the included code, the functions are named t_expf(), t_exp(),
and t_expl(), which permits inclusion into libm for ease of
testing. If these implementations are to supercede the versions
available in libm, then rename the functions accordingly.
>How-To-Repeat:
>Fix:
/*-
* Copyright (c) 2009-2010 Steven G. Kargl
* 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 unmodified, 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 ``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 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.
*/
/*
* Single precision implemenation of exp(x) as outlined in
*
* PTP Tang, "Table-driven implementation of the exponential function
* in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
* 144-157 (1989).
*/
#include <sys/cdefs.h>
#include <math.h>
#include "math_private.h"
/*
* Tang states that THRESHOLD_1 is 341*log(2), which is about 236.36. But,
* the constant given in the Appendix is 0x435c6bba, which is about 220.42.
*/
#define THRESHOLD_1 0x1.d8b9f4p+7f /* 341 * log(2) */
#define THRESHOLD_2 0x1.000000p-25f /* 2**(-25) */
#define INV_L 0x1.715476p+5f /* 32 / log(2) */
#define L1 0x1.62e400p-6f /* L1+L2=log(2)/32 */
#define L2 0x1.7f7d1cp-25f /* L1+L2=log(2)/32 */
/* Remes polynomial coefficients. */
#define A1 0x1.000088p-1f
#define A2 0x1.5555d8p-3f
/* Table of values for 2**(j/32) with j = 0, ..., 31. */
static const struct {
float head;
float tail;
} s[32] = {
0x1.00000p+0f, 0x0.000000p+00f,
0x1.059b0p+0f, 0x1.a62b0ap-21f,
0x1.0b558p+0f, 0x1.b3e624p-22f,
0x1.11300p+0f, 0x1.d0125cp-20f,
0x1.172b8p+0f, 0x1.e3ea8cp-23f,
0x1.1d480p+0f, 0x1.cc5a2ep-18f,
0x1.23878p+0f, 0x1.373ab2p-19f,
0x1.29e98p+0f, 0x1.7d47f8p-18f,
0x1.306f8p+0f, 0x1.828c6ep-18f,
0x1.371a0p+0f, 0x1.cdceaap-18f,
0x1.3dea0p+0f, 0x1.93048ep-18f,
0x1.44e08p+0f, 0x1.818624p-22f,
0x1.4bfd8p+0f, 0x1.6a9b16p-19f,
0x1.53428p+0f, 0x1.ab4ea8p-19f,
0x1.5ab00p+0f, 0x1.f75216p-18f,
0x1.62478p+0f, 0x1.ac0e96p-18f,
0x1.6a098p+0f, 0x1.999fcep-18f,
0x1.71f70p+0f, 0x1.7a3b18p-18f,
0x1.7a110p+0f, 0x1.1cfac0p-18f,
0x1.82588p+0f, 0x1.994ccep-20f,
0x1.8ace0p+0f, 0x1.508aa8p-18f,
0x1.93730p+0f, 0x1.ec3372p-18f,
0x1.9c490p+0f, 0x1.82a3f0p-20f,
0x1.a5500p+0f, 0x1.d91f12p-19f,
0x1.ae898p+0f, 0x1.e656b4p-18f,
0x1.b7f70p+0f, 0x1.bcbed8p-18f,
0x1.c1998p+0f, 0x1.eec2aap-19f,
0x1.cb720p+0f, 0x1.b9df20p-21f,
0x1.d5818p+0f, 0x1.b9f74ap-21f,
0x1.dfc90p+0f, 0x1.ccdee6p-18f,
0x1.ea4a8p+0f, 0x1.e8a924p-18f,
0x1.f5070p+0f, 0x1.96db92p-18f
};
static const float tiny = 0x1.p-100f;
/*
* From Tang's paper: INTRND rounds a floating-point number to the nearest
* integer in the manner prescribed by the IEEE standard. It is crucial that
* the default round-to-nearest mode, not any other rounding mode, is in
* effect here.
*/
#define INTRND(x) roundf((x))
float
t_expf(float x)
{
int32_t xsb;
u_int32_t hx;
int n, n1, n2;
float ax, r1, r2, r, p, q, t;
GET_FLOAT_WORD(hx,x);
/* Determine the sign bit of x. */
xsb = (hx >> 31) & 1;
/* Set high word to |x|. */
hx &= 0x7fffffff;
/* Test for NaN and +-Inf. */
if (hx >= 0x7f800000) {
/* exp(NaN) = NaN. */
if (hx > 0x7f800000)
return (x + x);
/* exp(+inf) = inf (no exception!) or exp(-inf) = 0. */
return (xsb ? 0 : x);
}
SET_FLOAT_WORD(ax, hx);
if (ax > THRESHOLD_1)
return (x > 0 ? (ax / (ax - ax)) : (0 + tiny * tiny));
if (ax < THRESHOLD_2)
return (1 + x);
n = INTRND(x * INV_L);
n2 = n % 32;
if (n2 < 0) n2 += 32; /* Tang's j. */
n1 = n - n2;
if (n >= 512 || n <= -512)
r1 = (x - n1 * L1) - n2 * L1;
else
r1 = x - n * L1;
r2 = - n * L2;
n1 /= 32; /* Tang's m. */
r = r1 + r2;
q = r * r * (A1 + r * A2);
p = r1 + (r2 + q);
t = s[n2].head + s[n2].tail;
t = s[n2].head + (s[n2].tail + t * p);
return (ldexpf(t, n1));
}
/*-
* Copyright (c) 2009-2010 Steven G. Kargl
* 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 unmodified, 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 ``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 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.
*/
/*
* Double precision implemenation of exp(x) as outlined in
*
* PTP Tang, "Table-driven implementation of the exponential function
* in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
* 144-157 (1989).
*
*/
#include <sys/cdefs.h>
#include "math.h"
#include "math_private.h"
#define THRESHOLD_1 0x1.c4474e1726454p+10 /* 2610 * log(2) */
#define THRESHOLD_2 0x1.0000000000000p-54 /* 2**(-54) */
#define INV_L 0x1.71547652b82fep+5 /* 32 / log(2) */
#define L1 0x1.62e42fef00000p-6 /* L1+L2=log(2)/32 */
#define L2 0x1.473de6af278ecp-39 /* L1+L2=log(2)/32 */
/* Remes polynomial coefficients. */
#define A1 0x1.0000000000000p-1
#define A2 0x1.5555555548f7cp-3
#define A3 0x1.5555555545d4ep-5
#define A4 0x1.11115b7aa905ep-7
#define A5 0x1.6c1728d739764p-10
static const struct {
double head;
double tail;
} s[32] = {
0x1p+0, 0x0p+0,
0x1.059b0d315854p+0, 0x1.a1d73e2a475b4p-47,
0x1.0b5586cf9890p+0, 0x1.ec5317256e308p-49,
0x1.11301d0125b4p+0, 0x1.0a4ebbf1aed93p-48,
0x1.172b83c7d514p+0, 0x1.d6e6fbe462876p-47,
0x1.1d4873168b98p+0, 0x1.53c02dc0144c8p-47,
0x1.2387a6e75620p+0, 0x1.c3360fd6d8e0bp-47,
0x1.29e9df51fdecp+0, 0x1.09612e8afad12p-47,
0x1.306fe0a31b70p+0, 0x1.52de8d5a46306p-48,
0x1.371a7373aa9cp+0, 0x1.54e28aa05e8a9p-49,
0x1.3dea64c12340p+0, 0x1.11ada0911f09fp-47,
0x1.44e086061890p+0, 0x1.68189b7a04ef8p-47,
0x1.4bfdad5362a0p+0, 0x1.38ea1cbd7f621p-47,
0x1.5342b569d4f8p+0, 0x1.df0a83c49d86ap-52,
0x1.5ab07dd48540p+0, 0x1.4ac64980a8c8fp-47,
0x1.6247eb03a558p+0, 0x1.2c7c3e81bf4b7p-50,
0x1.6a09e667f3bcp+0, 0x1.921165f626cddp-49,
0x1.71f75e8ec5f4p+0, 0x1.9ee91b8797785p-47,
0x1.7a11473eb018p+0, 0x1.b5f54408fdb37p-50,
0x1.82589994cce0p+0, 0x1.28acf88afab35p-48,
0x1.8ace5422aa0cp+0, 0x1.b5ba7c55a192dp-48,
0x1.93737b0cdc5cp+0, 0x1.27a280e1f92a0p-47,
0x1.9c49182a3f08p+0, 0x1.01c7c46b071f3p-48,
0x1.a5503b23e254p+0, 0x1.c8b424491caf8p-48,
0x1.ae89f995ad38p+0, 0x1.6af439a68bb99p-47,
0x1.b7f76f2fb5e4p+0, 0x1.baa9ec206ad4fp-50,
0x1.c199bdd85528p+0, 0x1.c2220cb12a092p-48,
0x1.cb720dcef904p+0, 0x1.48a81e5e8f4a5p-47,
0x1.d5818dcfba48p+0, 0x1.c976816bad9b8p-50,
0x1.dfc97337b9b4p+0, 0x1.eb968cac39ed3p-48,
0x1.ea4afa2a490cp+0, 0x1.9858f73a18f5ep-48,
0x1.f50765b6e454p+0, 0x1.9d3e12dd8a18bp-54
};
static const double tiny = 0x1.p-1000;
/*
* From Tang's paper: INTRND rounds a floating-point number to the nearest
* integer in the manner prescribed by the IEEE standard. It is crucial that
* the default round-to-nearest mode, not any other rounding mode, is in
* effect here.
*/
#define INTRND(x) round((x))
double
t_exp(double x)
{
u_int32_t hx, lx, xsb;
int n, n1, n2;
double ax, r1, r2, r, p, q, t;
EXTRACT_WORDS(hx, lx, x);
/* Detemine the sign bit of x. */
xsb = (hx >> 31) & 1;
/* Set high word of |x| */
hx &= 0x7fffffff;
/* Test for NaN and +-Inf. */
if (hx >= 0x7ff00000) {
/* exp(NaN) = NaN. */
if (((hx & 0xfffff) | lx) != 0)
return (x + x);
/* exp(+inf) = inf (no exception!) or exp(-inf) = 0. */
return (xsb ? 0 : x);
}
ax = x;
SET_HIGH_WORD(ax, hx);
if (ax > THRESHOLD_1)
return (x > 0 ?(ax / (ax - ax)) : (0 + tiny * tiny));
if (ax < THRESHOLD_2)
return (1 + x);
n = INTRND(x * INV_L);
n2 = n % 32;
if (n2 < 0) n2 += 32; /* Tang's j. */
n1 = n - n2;
if (n >= 512 || n <= -512)
r1 = (x - n1 * L1) - n2 * L1;
else
r1 = x - n * L1;
r2 = - n * L2;
n1 /= 32; /* Tang's m. */
r = r1 + r2;
q = r * r * (A1 + r * (A2 + r * (A3 + r * (A4 + r * A5))));
p = r1 + (r2 + q);
t = s[n2].head + s[n2].tail;
t = s[n2].head + (s[n2].tail + t * p);
return (ldexp(t, n1));
}
/*-
* Copyright (c) 2009-2010 Steven G. Kargl
* 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 unmodified, 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 ``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 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.
*/
/*
* Implemenation of expl(x) in Intel 80-bit long double format.
* The significand is 64 bits. This is based on
*
* PTP Tang, "Table-driven implementation of the exponential function
* in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
* 144-157 (1989).
*
* where the 32 table entries has been expanded to NUM (see below).
*
*/
#include <float.h>
#include "math.h"
#include "fpmath.h"
#define BIAS (LDBL_MAX_EXP - 1)
/*
* The range of overflow to +Inf and of underflow to +0 is given by
* 41021 * log(2). 41021 is obtained from -(-16382-64+1)+24576 where
* -16382 is the minimum exponent, 64 is the precision and 24576 is
* a bias adjustment given in IEEE 754 as 3*2^(e-2) with e being the
* width of the exponent in bits.
*/
static const long double
THRESHOLD_1 = 2.84335904937495165381e+04L, /* 41021 * log(2) */
THRESHOLD_2 = 2.71050543121376108502e-20L, /* 2**(-65) */
L1 = 5.4152123481245872938e-03L, /* L1+L2=log(2)/128 */
L2 = -1.4563974013471118145e-17L, /* L1+L2=log(2)/128 */
INV_L = 1.84664965233787316146e+02L; /* 128 / log(2) */
/*
* The polynomial coefficients are the initial rational approximation
* of the minimax polynomial coefficients determined by a Remes algorithm
* coded in a much higher working precision. For details, see
*
* N. Brisebarre, et al, "Computing Machine-Efficient Polynomial
* Approximations," ACM Tran. Math. Soft., 32, 236-256 (2006).
*/
static const long double
P1 = 5.000000000000000000e-01L,
P2 = 1.666666666666666666e-01L,
P3 = 4.166666666666666677e-02L,
P4 = 8.333333333333284767e-03L,
P5 = 1.388888888835006555e-03L,
P6 = 1.984127249252972877e-04L,
P7 = 2.480718929383586384e-05L;
/*
* 2^(i/NUM) for i in [0,NUM] is represented by two values where the
* first 47-bits of the significand is stored in head and the next 53
* bits are in tail.
*/
#define NUM 128
static const struct {
double head;
double tail;
} s[NUM] = {
0x1p+0, 0x0p+0,
0x1.0163da9fb330p+0, 0x1.ab6c25335719bp-47,
0x1.02c9a3e77804p+0, 0x1.07737be56527cp-47,
0x1.04315e86e7f8p+0, 0x1.2f5ce3e688369p-50,
0x1.059b0d315854p+0, 0x1.a1d73e2a475b4p-47,
0x1.0706b29ddf6cp+0, 0x1.dc6dc403a9d88p-48,
0x1.0874518759bcp+0, 0x1.01186be4bb285p-49,
0x1.09e3ecac6f38p+0, 0x1.a290f03062c27p-51,
0x1.0b5586cf9890p+0, 0x1.ec5317256e308p-49,
0x1.0cc922b7247cp+0, 0x1.ba03db82dc49fp-47,
0x1.0e3ec32d3d18p+0, 0x1.10103a1727c58p-47,
0x1.0fb66affed30p+0, 0x1.af232091dd8a1p-48,
0x1.11301d0125b4p+0, 0x1.0a4ebbf1aed93p-48,
0x1.12abdc06c31cp+0, 0x1.7f72575a649adp-49,
0x1.1429aaea92dcp+0, 0x1.fb34101943b26p-48,
0x1.15a98c8a58e4p+0, 0x1.12480d573dd56p-48,
0x1.172b83c7d514p+0, 0x1.d6e6fbe462876p-47,
0x1.18af9388c8dcp+0, 0x1.4dddfb85cd1e1p-47,
0x1.1a35beb6fcb4p+0, 0x1.a9e5b4c7b4969p-47,
0x1.1bbe084045ccp+0, 0x1.39ab1e72b4428p-48,
0x1.1d4873168b98p+0, 0x1.53c02dc0144c8p-47,
0x1.1ed5022fcd90p+0, 0x1.cb8819ff61122p-48,
0x1.2063b88628ccp+0, 0x1.63b8eeb029509p-48,
0x1.21f49917ddc8p+0, 0x1.62552fd29294cp-48,
0x1.2387a6e75620p+0, 0x1.c3360fd6d8e0bp-47,
0x1.251ce4fb2a60p+0, 0x1.f9ac155bef4f5p-47,
0x1.26b4565e27ccp+0, 0x1.d257a673281d4p-48,
0x1.284dfe1f5638p+0, 0x1.2d9e2b9e07941p-53,
0x1.29e9df51fdecp+0, 0x1.09612e8afad12p-47,
0x1.2b87fd0dad98p+0, 0x1.ffbbd48ca71f9p-49,
0x1.2d285a6e4030p+0, 0x1.680123aa6da0fp-49,
0x1.2ecafa93e2f4p+0, 0x1.611ca0f45d524p-48,
0x1.306fe0a31b70p+0, 0x1.52de8d5a46306p-48,
0x1.32170fc4cd80p+0, 0x1.89a9ce78e1804p-47,
0x1.33c08b26416cp+0, 0x1.fa64e43086cb3p-47,
0x1.356c55f929fcp+0, 0x1.864a311a3b1bap-47,
0x1.371a7373aa9cp+0, 0x1.54e28aa05e8a9p-49,
0x1.38cae6d05d84p+0, 0x1.2c2d4e586cdf7p-47,
0x1.3a7db34e59fcp+0, 0x1.b750de494cf05p-47,
0x1.3c32dc313a8cp+0, 0x1.242000f9145acp-47,
0x1.3dea64c12340p+0, 0x1.11ada0911f09fp-47,
0x1.3fa4504ac800p+0, 0x1.ba0bf701aa418p-48,
0x1.4160a21f72e0p+0, 0x1.4fc2192dc79eep-47,
0x1.431f5d950a88p+0, 0x1.6dc704439410dp-48,
0x1.44e086061890p+0, 0x1.68189b7a04ef8p-47,
0x1.46a41ed1d004p+0, 0x1.772512f45922ap-48,
0x1.486a2b5c13ccp+0, 0x1.013c1a3b69063p-48,
0x1.4a32af0d7d3cp+0, 0x1.e672d8bcf46f9p-48,
0x1.4bfdad5362a0p+0, 0x1.38ea1cbd7f621p-47,
0x1.4dcb299fddd0p+0, 0x1.ac766dde353c2p-49,
0x1.4f9b2769d2c8p+0, 0x1.35699ec5b4d50p-47,
0x1.516daa2cf664p+0, 0x1.c112f52c84d82p-52,
0x1.5342b569d4f8p+0, 0x1.df0a83c49d86ap-52,
0x1.551a4ca5d920p+0, 0x1.d8a5d8c40486ap-49,
0x1.56f4736b527cp+0, 0x1.a66ecb004764fp-48,
0x1.58d12d497c7cp+0, 0x1.e9295e15b9a1ep-47,
0x1.5ab07dd48540p+0, 0x1.4ac64980a8c8fp-47,
0x1.5c9268a59468p+0, 0x1.b80e258dc0b4cp-47,
0x1.5e76f15ad214p+0, 0x1.0dd37c9840733p-49,
0x1.605e1b976dc0p+0, 0x1.160edeb25490ep-49,
0x1.6247eb03a558p+0, 0x1.2c7c3e81bf4b7p-50,
0x1.6434634ccc30p+0, 0x1.fc76f8714c4eep-48,
0x1.662388255220p+0, 0x1.24893ecf14dc8p-47,
0x1.68155d44ca94p+0, 0x1.9840e2b913dd0p-47,
0x1.6a09e667f3bcp+0, 0x1.921165f626cddp-49,
0x1.6c012750bda8p+0, 0x1.f76bb54cc007ap-47,
0x1.6dfb23c651a0p+0, 0x1.779107165f0dep-47,
0x1.6ff7df951948p+0, 0x1.e7c3f0da79f11p-51,
0x1.71f75e8ec5f4p+0, 0x1.9ee91b8797785p-47,
0x1.73f9a48a5814p+0, 0x1.9deae4d273456p-47,
0x1.75feb564267cp+0, 0x1.17edd35467491p-49,
0x1.780694fde5d0p+0, 0x1.fb0cd7014042cp-47,
0x1.7a11473eb018p+0, 0x1.b5f54408fdb37p-50,
0x1.7c1ed0130c10p+0, 0x1.93e2499a22c9cp-47,
0x1.7e2f336cf4e4p+0, 0x1.1082e815d0abdp-47,
0x1.80427543e1a0p+0, 0x1.1b60de67649a3p-48,
0x1.82589994cce0p+0, 0x1.28acf88afab35p-48,
0x1.8471a4623c78p+0, 0x1.667297b5cbe32p-47,
0x1.868d99b4492cp+0, 0x1.640720ec85613p-47,
0x1.88ac7d98a668p+0, 0x1.966530bcdf2d5p-48,
0x1.8ace5422aa0cp+0, 0x1.b5ba7c55a192dp-48,
0x1.8cf3216b5448p+0, 0x1.7de55439a2c39p-49,
0x1.8f1ae9915770p+0, 0x1.b15cc13a2e397p-47,
0x1.9145b0b91ffcp+0, 0x1.622986d1a7daep-50,
0x1.93737b0cdc5cp+0, 0x1.27a280e1f92a0p-47,
0x1.95a44cbc8520p+0, 0x1.dd36906d2b420p-49,
0x1.97d829fde4e4p+0, 0x1.f173d241f23d1p-49,
0x1.9a0f170ca078p+0, 0x1.cdd1884dc6234p-47,
0x1.9c49182a3f08p+0, 0x1.01c7c46b071f3p-48,
0x1.9e86319e3230p+0, 0x1.18c12653c7326p-47,
0x1.a0c667b5de54p+0, 0x1.2594d6d45c656p-47,
0x1.a309bec4a2d0p+0, 0x1.9ac60b8fbb86dp-47,
0x1.a5503b23e254p+0, 0x1.c8b424491caf8p-48,
0x1.a799e1330b34p+0, 0x1.86f2dfb2b158fp-48,
0x1.a9e6b5579fd8p+0, 0x1.fa1f5921deffap-47,
0x1.ac36bbfd3f34p+0, 0x1.ce06dcb351893p-47,
0x1.ae89f995ad38p+0, 0x1.6af439a68bb99p-47,
0x1.b0e07298db64p+0, 0x1.2c8421566fe38p-47,
0x1.b33a2b84f15cp+0, 0x1.d7b5fe873decap-47,
0x1.b59728de5590p+0, 0x1.cc71c40888b24p-47,
0x1.b7f76f2fb5e4p+0, 0x1.baa9ec206ad4fp-50,
0x1.ba5b030a1064p+0, 0x1.30819678d5eb7p-49,
0x1.bcc1e904bc1cp+0, 0x1.2247ba0f45b3dp-48,
0x1.bf2c25bd71e0p+0, 0x1.10811ae04a31cp-49,
0x1.c199bdd85528p+0, 0x1.c2220cb12a092p-48,
0x1.c40ab5fffd04p+0, 0x1.d368a6fc1078cp-47,
0x1.c67f12e57d14p+0, 0x1.694426ffa41e5p-49,
0x1.c8f6d9406e78p+0, 0x1.a88d65e24402ep-47,
0x1.cb720dcef904p+0, 0x1.48a81e5e8f4a5p-47,
0x1.cdf0b555dc3cp+0, 0x1.ce227c4ac7d63p-47,
0x1.d072d4a07894p+0, 0x1.dc68791790d0bp-47,
0x1.d2f87080d89cp+0, 0x1.8c56f091cc4f5p-47,
0x1.d5818dcfba48p+0, 0x1.c976816bad9b8p-50,
0x1.d80e316c9838p+0, 0x1.7bb84f9d04880p-48,
0x1.da9e603db328p+0, 0x1.5c2300696db53p-50,
0x1.dd321f301b44p+0, 0x1.025b4aef1e032p-47,
0x1.dfc97337b9b4p+0, 0x1.eb968cac39ed3p-48,
0x1.e264614f5a10p+0, 0x1.45093b0fd0bd7p-47,
0x1.e502ee78b3fcp+0, 0x1.b139e8980a9cdp-47,
0x1.e7a51fbc74c8p+0, 0x1.a5aa4594191bcp-51,
0x1.ea4afa2a490cp+0, 0x1.9858f73a18f5ep-48,
0x1.ecf482d8e67cp+0, 0x1.846d81897dca5p-47,
0x1.efa1bee615a0p+0, 0x1.3bb8fe90d496dp-47,
0x1.f252b376bba8p+0, 0x1.74e8696fc3639p-48,
0x1.f50765b6e454p+0, 0x1.9d3e12dd8a18bp-54,
0x1.f7bfdad9cbe0p+0, 0x1.38913b4bfe72cp-48,
0x1.fa7c1819e90cp+0, 0x1.82e90a7e74b26p-48,
0x1.fd3c22b8f71cp+0, 0x1.884badd25995ep-47
};
static const long double tiny = 0x1.p-10000L;
/*
* From Tang's paper: INTRND rounds a floating-point number to the nearest
* integer in the manner prescribed by the IEEE standard. It is crucial that
* the default round-to-nearest mode, not any other rounding mode, is in
* effect here.
*/
#define INTRND(x) roundl((x))
long double
t_expl(long double x)
{
union IEEEl2bits z;
int k, n, n1, n2, sgn;
long double r1, r2, r, p, q, t;
z.e = x;
sgn = z.bits.sign;
z.bits.sign = 0;
/*
* If x = +Inf, then exp(x) = Inf.
* If x = -Inf, then exp(x) = 0.
* If x = NaN, then exp(x) = NaN.
*/
if (z.bits.exp == LDBL_MAX_EXP + BIAS) {
mask_nbit_l(z);
if (!(z.bits.manh | z.bits.manl))
return (sgn ? 0 : x);
return (x + x);
}
if (z.e > THRESHOLD_1)
return (x > 0 ? (z.e / (z.e - z.e)) : (0 + tiny * tiny));
if (z.e < THRESHOLD_2)
return (1 + x);
n = INTRND(x * INV_L);
n2 = n % NUM;
if (n2 < 0) n2 += NUM; /* Tang's j. */
n1 = n - n2;
if (n >= 512 || n <= -512)
r1 = (x - n1 * L1) - n2 * L1;
else
r1 = x - n * L1;
r2 = - n * L2;
n1 /= NUM; /* Tang's m. */
r = r1 + r2;
q = r * r * (P1 + r * (P2 + r * (P3 + r * (P4 + r * (P5 +
r * (P6 + r * P7))))));
p = r1 + (r2 + q);
t = (long double)s[n2].tail + s[n2].head;
t = (long double)s[n2].head + (s[n2].tail + t * p);
return (ldexpl(t, n1));
}
>Release-Note:
>Audit-Trail:
>Unformatted:
More information about the freebsd-standards
mailing list