* ld80/s_expl.c:
. Guard a comment from reformatting by indent(1). . Re-order variables in declarations to alphabetical order. . Remove a banal comment. * ld128/s_expl.c: . Add a comment to point to ld80/s_expl.c for implementation details. . Move the #define of INTERVAL to reduce the diff with ld80/s_expl.c. . twom10000 does not need to be volatile, so move its declaration. . Re-order variables in declarations to alphabetical order. . Add a comment that describes the argument reduction. . Remove the same banal comment found in ld80/s_expl.c. Reviewed by: bde Approved by: das (mentor)
This commit is contained in:
parent
0ad5179dee
commit
8f647ffd7f
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=240864
@ -27,24 +27,29 @@
|
||||
#include <sys/cdefs.h>
|
||||
__FBSDID("$FreeBSD$");
|
||||
|
||||
/*
|
||||
* ld128 version of s_expl.c. See ../ld80/s_expl.c for most comments.
|
||||
*/
|
||||
|
||||
#include <float.h>
|
||||
|
||||
#include "fpmath.h"
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#define INTERVALS 128
|
||||
#define BIAS (LDBL_MAX_EXP - 1)
|
||||
|
||||
/* XXX Prevent gcc from erroneously constant folding this: */
|
||||
static volatile const long double twom10000 = 0x1p-10000L, tiny = 0x1p-10000L;
|
||||
static volatile const long double tiny = 0x1p-10000L;
|
||||
|
||||
static const long double
|
||||
huge = 0x1p10000L,
|
||||
o_threshold = 11356.523406294143949491931077970763428L,
|
||||
u_threshold = -11433.462743336297878837243843452621503L,
|
||||
INV_L = 1.84664965233787316142070359168242182e+02L,
|
||||
L1 = 5.41521234812457272982212595914567508e-03L,
|
||||
L2 = -1.02536706388947310094527932552595546e-29L,
|
||||
INV_L = 1.84664965233787316142070359168242182e+02L;
|
||||
huge = 0x1p10000L,
|
||||
o_threshold = 11356.523406294143949491931077970763428L,
|
||||
twom10000 = 0x1p-10000L,
|
||||
u_threshold = -11433.462743336297878837243843452621503L;
|
||||
|
||||
static const long double
|
||||
P2 = 5.00000000000000000000000000000000000e-1L,
|
||||
@ -58,8 +63,6 @@ P9 = 2.75573192240103867817876199544468806e-6L,
|
||||
P10 = 2.75573236172670046201884000197885520e-7L,
|
||||
P11 = 2.50517544183909126492878226167697856e-8L;
|
||||
|
||||
#define INTERVALS 128
|
||||
|
||||
static const struct {
|
||||
long double hi;
|
||||
long double lo;
|
||||
@ -225,9 +228,10 @@ expl(long double x)
|
||||
return (1.0L + x);
|
||||
}
|
||||
|
||||
/* Reduce x to (k*ln2 + midpoint[n2] + r1 + r2). */
|
||||
fn = x * INV_L + 0x1.8p112 - 0x1.8p112;
|
||||
n = (int)fn;
|
||||
n2 = (unsigned)n % INTERVALS; /* Tang's j. */
|
||||
n2 = (unsigned)n % INTERVALS;
|
||||
k = (n - n2) / INTERVALS;
|
||||
r1 = x - fn * L1;
|
||||
r2 = -fn * L2;
|
||||
|
@ -29,7 +29,7 @@
|
||||
#include <sys/cdefs.h>
|
||||
__FBSDID("$FreeBSD$");
|
||||
|
||||
/*
|
||||
/*-
|
||||
* Compute the exponential of x for Intel 80-bit format. This is based on:
|
||||
*
|
||||
* PTP Tang, "Table-driven implementation of the exponential function
|
||||
@ -70,9 +70,9 @@ static const double __aligned(64)
|
||||
* have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
|
||||
* bits zero so that multiplication of it by n is exact.
|
||||
*/
|
||||
INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */
|
||||
L1 = 5.4152123484527692e-3, /* 0x162e42ff000000.0p-60 */
|
||||
L2 = -3.2819649005320973e-13, /* -0x1718432a1b0e26.0p-94 */
|
||||
INV_L = 1.8466496523378731e+2, /* 0x171547652b82fe.0p-45 */
|
||||
/*
|
||||
* Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]:
|
||||
* |exp(x) - p(x)| < 2**-77.2
|
||||
@ -231,7 +231,7 @@ long double
|
||||
expl(long double x)
|
||||
{
|
||||
union IEEEl2bits u, v;
|
||||
long double fn, r, r1, r2, q, t, t23, t45, twopk, twopkp10000, z;
|
||||
long double fn, q, r, r1, r2, t, t23, t45, twopk, twopkp10000, z;
|
||||
int k, n, n2;
|
||||
uint16_t hx, ix;
|
||||
|
||||
@ -268,7 +268,7 @@ expl(long double x)
|
||||
#else
|
||||
n = (int)fn;
|
||||
#endif
|
||||
n2 = (unsigned)n % INTERVALS; /* Tang's j. */
|
||||
n2 = (unsigned)n % INTERVALS;
|
||||
k = (n - n2) / INTERVALS;
|
||||
r1 = x - fn * L1;
|
||||
r2 = -fn * L2;
|
||||
|
Loading…
Reference in New Issue
Block a user