generic-library/bionic
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge "Take upstream libm changes."

Browse Source
This commit is contained in:
Elliott Hughes 2013-06-13 00:29:21 +00:00 committed by Gerrit Code Review
commit 92e841d0aa
21 changed files with 383 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
libm/upstream-freebsd/lib/msun/src/e_acosh.c
View File

@ -29,6 +29,8 @@ __FBSDID("$FreeBSD$");
* acosh(NaN) is NaN without signal.
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
@ -60,3 +62,7 @@ __ieee754_acosh(double x)
return log1p(t+sqrt(2.0*t+t*t));
}
}
#if LDBL_MANT_DIG == 53
__weak_reference(acosh, acoshl);
#endif

6
libm/upstream-freebsd/lib/msun/src/e_atanh.c
View File

@ -33,6 +33,8 @@ __FBSDID("$FreeBSD$");
*
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
@ -60,3 +62,7 @@ __ieee754_atanh(double x)
t = 0.5*log1p((x+x)/(one-x));
if(hx>=0) return t; else return -t;
}
#if LDBL_MANT_DIG == 53
__weak_reference(atanh, atanhl);
#endif

2
libm/upstream-freebsd/lib/msun/src/e_exp.c
View File

@ -84,7 +84,6 @@ __FBSDID("$FreeBSD$");
static const double
one = 1.0,
halF[2] = {0.5,-0.5,},
huge = 1.0e+300,
o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
@ -99,6 +98,7 @@ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
static volatile double
huge = 1.0e+300,
twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/
double

5
libm/upstream-freebsd/lib/msun/src/e_expf.c
View File

@ -24,7 +24,6 @@ __FBSDID("$FreeBSD$");
static const float
one = 1.0,
halF[2] = {0.5,-0.5,},
huge = 1.0e+30,
o_threshold= 8.8721679688e+01, /* 0x42b17180 */
u_threshold= -1.0397208405e+02, /* 0xc2cff1b5 */
ln2HI[2] ={ 6.9314575195e-01, /* 0x3f317200 */
@ -39,7 +38,9 @@ invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
P1 = 1.6666625440e-1, /* 0xaaaa8f.0p-26 */
P2 = -2.7667332906e-3; /* -0xb55215.0p-32 */
static volatile float twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */
static volatile float
huge = 1.0e+30,
twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */
float
__ieee754_expf(float x)

9
libm/upstream-freebsd/lib/msun/src/e_log.c
View File

@ -65,6 +65,8 @@ __FBSDID("$FreeBSD$");
* to produce the hexadecimal values shown.
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
@ -81,6 +83,7 @@ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
static volatile double vzero = 0.0;
double
__ieee754_log(double x)
@ -94,7 +97,7 @@ __ieee754_log(double x)
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
return -two54/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx,x);
@ -138,3 +141,7 @@ __ieee754_log(double x)
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(log, logl);
#endif

9
libm/upstream-freebsd/lib/msun/src/e_log10.c
View File

@ -22,6 +22,8 @@ __FBSDID("$FreeBSD$");
* in not-quite-routine extra precision.
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
#include "k_log.h"
@ -34,6 +36,7 @@ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
static const double zero = 0.0;
static volatile double vzero = 0.0;
double
__ieee754_log10(double x)
@ -47,7 +50,7 @@ __ieee754_log10(double x)
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
return -two54/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx,x);
@ -85,3 +88,7 @@ __ieee754_log10(double x)
return val_lo + val_hi;
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(log10, log10l);
#endif

3
libm/upstream-freebsd/lib/msun/src/e_log10f.c
View File

@ -28,6 +28,7 @@ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */
static const float zero = 0.0;
static volatile float vzero = 0.0;
float
__ieee754_log10f(float x)
@ -40,7 +41,7 @@ __ieee754_log10f(float x)
k=0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff)==0)
return -two25/zero; /* log(+-0)=-inf */
return -two25/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(hx,x);

9
libm/upstream-freebsd/lib/msun/src/e_log2.c
View File

@ -24,6 +24,8 @@ __FBSDID("$FreeBSD$");
* in not-quite-routine extra precision.
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
#include "k_log.h"
@ -34,6 +36,7 @@ ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
static const double zero = 0.0;
static volatile double vzero = 0.0;
double
__ieee754_log2(double x)
@ -47,7 +50,7 @@ __ieee754_log2(double x)
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
return -two54/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx,x);
@ -108,3 +111,7 @@ __ieee754_log2(double x)
return val_lo + val_hi;
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(log2, log2l);
#endif

3
libm/upstream-freebsd/lib/msun/src/e_log2f.c
View File

@ -26,6 +26,7 @@ ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
static const float zero = 0.0;
static volatile float vzero = 0.0;
float
__ieee754_log2f(float x)
@ -38,7 +39,7 @@ __ieee754_log2f(float x)
k=0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff)==0)
return -two25/zero; /* log(+-0)=-inf */
return -two25/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(hx,x);

3
libm/upstream-freebsd/lib/msun/src/e_logf.c
View File

@ -30,6 +30,7 @@ Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
static const float zero = 0.0;
static volatile float vzero = 0.0;
float
__ieee754_logf(float x)
@ -42,7 +43,7 @@ __ieee754_logf(float x)
k=0;
if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix&0x7fffffff)==0)
return -two25/zero; /* log(+-0)=-inf */
return -two25/vzero; /* log(+-0)=-inf */
if (ix<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(ix,x);

292
libm/upstream-freebsd/lib/msun/src/math_private.h
View File

@ -188,6 +188,33 @@ do { \
(d) = sf_u.value; \
} while (0)
/*
* Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
* double.
*/
#define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \
do { \
union IEEEl2bits ew_u; \
ew_u.e = (d); \
(ix0) = ew_u.xbits.expsign; \
(ix1) = ew_u.xbits.man; \
} while (0)
/*
* Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
* long double.
*/
#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \
do { \
union IEEEl2bits ew_u; \
ew_u.e = (d); \
(ix0) = ew_u.xbits.expsign; \
(ix1) = ew_u.xbits.manh; \
(ix2) = ew_u.xbits.manl; \
} while (0)
/* Get expsign as a 16 bit int from a long double. */
#define GET_LDBL_EXPSIGN(i,d) \
@ -197,6 +224,33 @@ do { \
(i) = ge_u.xbits.expsign; \
} while (0)
/*
* Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
* mantissa.
*/
#define INSERT_LDBL80_WORDS(d,ix0,ix1) \
do { \
union IEEEl2bits iw_u; \
iw_u.xbits.expsign = (ix0); \
iw_u.xbits.man = (ix1); \
(d) = iw_u.e; \
} while (0)
/*
* Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
* comprising the mantissa.
*/
#define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \
do { \
union IEEEl2bits iw_u; \
iw_u.xbits.expsign = (ix0); \
iw_u.xbits.manh = (ix1); \
iw_u.xbits.manl = (ix2); \
(d) = iw_u.e; \
} while (0)
/* Set expsign of a long double from a 16 bit int. */
#define SET_LDBL_EXPSIGN(d,v) \
@ -260,6 +314,110 @@ do { \
/* Default return statement if hack*_t() is not used. */
#define RETURNF(v) return (v)
/*
* 2sum gives the same result as 2sumF without requiring |a| >= |b| or
* a == 0, but is slower.
*/
#define _2sum(a, b) do { \
__typeof(a) __s, __w; \
\
__w = (a) + (b); \
__s = __w - (a); \
(b) = ((a) - (__w - __s)) + ((b) - __s); \
(a) = __w; \
} while (0)
/*
* 2sumF algorithm.
*
* "Normalize" the terms in the infinite-precision expression a + b for
* the sum of 2 floating point values so that b is as small as possible
* relative to 'a'. (The resulting 'a' is the value of the expression in
* the same precision as 'a' and the resulting b is the rounding error.)
* |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
* exponent overflow or underflow must not occur. This uses a Theorem of
* Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum"
* is apparently due to Skewchuk (1997).
*
* For this to always work, assignment of a + b to 'a' must not retain any
* extra precision in a + b. This is required by C standards but broken
* in many compilers. The brokenness cannot be worked around using
* STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
* algorithm would be destroyed by non-null strict assignments. (The
* compilers are correct to be broken -- the efficiency of all floating
* point code calculations would be destroyed similarly if they forced the
* conversions.)
*
* Fortunately, a case that works well can usually be arranged by building
* any extra precision into the type of 'a' -- 'a' should have type float_t,
* double_t or long double. b's type should be no larger than 'a's type.
* Callers should use these types with scopes as large as possible, to
* reduce their own extra-precision and efficiciency problems. In
* particular, they shouldn't convert back and forth just to call here.
*/
#ifdef DEBUG
#define _2sumF(a, b) do { \
__typeof(a) __w; \
volatile __typeof(a) __ia, __ib, __r, __vw; \
\
__ia = (a); \
__ib = (b); \
assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \
\
__w = (a) + (b); \
(b) = ((a) - __w) + (b); \
(a) = __w; \
\
/* The next 2 assertions are weak if (a) is already long double. */ \
assert((long double)__ia + __ib == (long double)(a) + (b)); \
__vw = __ia + __ib; \
__r = __ia - __vw; \
__r += __ib; \
assert(__vw == (a) && __r == (b)); \
} while (0)
#else /* !DEBUG */
#define _2sumF(a, b) do { \
__typeof(a) __w; \
\
__w = (a) + (b); \
(b) = ((a) - __w) + (b); \
(a) = __w; \
} while (0)
#endif /* DEBUG */
/*
* Set x += c, where x is represented in extra precision as a + b.
* x must be sufficiently normalized and sufficiently larger than c,
* and the result is then sufficiently normalized.
*
* The details of ordering are that |a| must be >= |c| (so that (a, c)
* can be normalized without extra work to swap 'a' with c). The details of
* the normalization are that b must be small relative to the normalized 'a'.
* Normalization of (a, c) makes the normalized c tiny relative to the
* normalized a, so b remains small relative to 'a' in the result. However,
* b need not ever be tiny relative to 'a'. For example, b might be about
* 2**20 times smaller than 'a' to give about 20 extra bits of precision.
* That is usually enough, and adding c (which by normalization is about
* 2**53 times smaller than a) cannot change b significantly. However,
* cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
* significantly relative to b. The caller must ensure that significant
* cancellation doesn't occur, either by having c of the same sign as 'a',
* or by having |c| a few percent smaller than |a|. Pre-normalization of
* (a, b) may help.
*
* This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
* exercise 19). We gain considerable efficiency by requiring the terms to
* be sufficiently normalized and sufficiently increasing.
*/
#define _3sumF(a, b, c) do { \
__typeof(a) __tmp; \
\
__tmp = (c); \
_2sumF(__tmp, (a)); \
(b) += (a); \
(a) = __tmp; \
} while (0)
/*
* Common routine to process the arguments to nan(), nanf(), and nanl().
*/
@ -370,6 +528,140 @@ irintl(long double x)
#endif /* __GNUCLIKE_ASM */
#ifdef DEBUG
#if defined(__amd64__) || defined(__i386__)
#define breakpoint() asm("int $3")
#else
#include <signal.h>
#define breakpoint() raise(SIGTRAP)
#endif
#endif
/* Write a pari script to test things externally. */
#ifdef DOPRINT
#include <stdio.h>
#ifndef DOPRINT_SWIZZLE
#define DOPRINT_SWIZZLE 0
#endif
#ifdef DOPRINT_LD80
#define DOPRINT_START(xp) do { \
uint64_t __lx; \
uint16_t __hx; \
\
/* Hack to give more-problematic args. */ \
EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
__lx ^= DOPRINT_SWIZZLE; \
INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
printf("x = %.21Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#elif defined(DOPRINT_D64)
#define DOPRINT_START(xp) do { \
uint32_t __hx, __lx; \
\
EXTRACT_WORDS(__hx, __lx, *xp); \
__lx ^= DOPRINT_SWIZZLE; \
INSERT_WORDS(*xp, __hx, __lx); \
printf("x = %.21Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#elif defined(DOPRINT_F32)
#define DOPRINT_START(xp) do { \
uint32_t __hx; \
\
GET_FLOAT_WORD(__hx, *xp); \
__hx ^= DOPRINT_SWIZZLE; \
SET_FLOAT_WORD(*xp, __hx); \
printf("x = %.21Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
#ifndef DOPRINT_SWIZZLE_HIGH
#define DOPRINT_SWIZZLE_HIGH 0
#endif
#define DOPRINT_START(xp) do { \
uint64_t __lx, __llx; \
uint16_t __hx; \
\
EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
__llx ^= DOPRINT_SWIZZLE; \
__lx ^= DOPRINT_SWIZZLE_HIGH; \
INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
printf("x = %.36Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#endif /* DOPRINT_LD80 */
#else /* !DOPRINT */
#define DOPRINT_START(xp)
#define DOPRINT_END1(v)
#define DOPRINT_END2(hi, lo)
#endif /* DOPRINT */
#define RETURNP(x) do { \
DOPRINT_END1(x); \
RETURNF(x); \
} while (0)
#define RETURNPI(x) do { \
DOPRINT_END1(x); \
RETURNI(x); \
} while (0)
#define RETURN2P(x, y) do { \
DOPRINT_END2((x), (y)); \
RETURNF((x) + (y)); \
} while (0)
#define RETURN2PI(x, y) do { \
DOPRINT_END2((x), (y)); \
RETURNI((x) + (y)); \
} while (0)
#ifdef STRUCT_RETURN
#define RETURNSP(rp) do { \
if (!(rp)->lo_set) \
RETURNP((rp)->hi); \
RETURN2P((rp)->hi, (rp)->lo); \
} while (0)
#define RETURNSPI(rp) do { \
if (!(rp)->lo_set) \
RETURNPI((rp)->hi); \
RETURN2PI((rp)->hi, (rp)->lo); \
} while (0)
#endif
#define SUM2P(x, y) ({ \
const __typeof (x) __x = (x); \
const __typeof (y) __y = (y); \
\
DOPRINT_END2(__x, __y); \
__x + __y; \
})
/*
* ieee style elementary functions
*

6
libm/upstream-freebsd/lib/msun/src/s_asinh.c
View File

@ -24,6 +24,8 @@ __FBSDID("$FreeBSD$");
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
@ -54,3 +56,7 @@ asinh(double x)
}
if(hx>0) return w; else return -w;
}
#if LDBL_MANT_DIG == 53
__weak_reference(asinh, asinhl);
#endif

5
libm/upstream-freebsd/lib/msun/src/s_exp2.c
View File

@ -36,7 +36,6 @@ __FBSDID("$FreeBSD$");
#define TBLSIZE (1 << TBLBITS)
static const double
huge = 0x1p1000,
redux = 0x1.8p52 / TBLSIZE,
P1 = 0x1.62e42fefa39efp-1,
P2 = 0x1.ebfbdff82c575p-3,
@ -44,7 +43,9 @@ static const double
P4 = 0x1.3b2ab88f70400p-7,
P5 = 0x1.5d88003875c74p-10;
static volatile double twom1000 = 0x1p-1000;
static volatile double
huge = 0x1p1000,
twom1000 = 0x1p-1000;
static const double tbl[TBLSIZE * 2] = {
/* exp2(z + eps) eps */

5
libm/upstream-freebsd/lib/msun/src/s_exp2f.c
View File

@ -36,14 +36,15 @@ __FBSDID("$FreeBSD$");
#define TBLSIZE (1 << TBLBITS)
static const float
huge = 0x1p100f,
redux = 0x1.8p23f / TBLSIZE,
P1 = 0x1.62e430p-1f,
P2 = 0x1.ebfbe0p-3f,
P3 = 0x1.c6b348p-5f,
P4 = 0x1.3b2c9cp-7f;
static volatile float twom100 = 0x1p-100f;
static volatile float
huge = 0x1p100f,
twom100 = 0x1p-100f;
static const double exp2ft[TBLSIZE] = {
0x1.6a09e667f3bcdp-1,

7
libm/upstream-freebsd/lib/msun/src/s_expm1.c
View File

@ -115,7 +115,6 @@ __FBSDID("$FreeBSD$");
static const double
one = 1.0,
huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
@ -128,6 +127,8 @@ Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
static volatile double huge = 1.0e+300;
double
expm1(double x)
{
@ -215,3 +216,7 @@ expm1(double x)
}
return y;
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(expm1, expm1l);
#endif

3
libm/upstream-freebsd/lib/msun/src/s_expm1f.c
View File

@ -23,7 +23,6 @@ __FBSDID("$FreeBSD$");
static const float
one = 1.0,
huge = 1.0e+30,
tiny = 1.0e-30,
o_threshold = 8.8721679688e+01,/* 0x42b17180 */
ln2_hi = 6.9313812256e-01,/* 0x3f317180 */
@ -37,6 +36,8 @@ invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
Q1 = -3.3333212137e-2, /* -0x888868.0p-28 */
Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */
static volatile float huge = 1.0e+30;
float
expm1f(float x)
{

8
libm/upstream-freebsd/lib/msun/src/s_fma.c
View File

@ -238,6 +238,8 @@ fma(double x, double y, double z)
zs = copysign(DBL_MIN, zs);
fesetround(FE_TONEAREST);
/* work around clang bug 8100 */
volatile double vxs = xs;
/*
* Basic approach for round-to-nearest:
@ -247,7 +249,7 @@ fma(double x, double y, double z)
* adj = xy.lo + r.lo (inexact; low bit is sticky)
* result = r.hi + adj (correctly rounded)
*/
xy = dd_mul(xs, ys);
xy = dd_mul(vxs, ys);
r = dd_add(xy.hi, zs);
spread = ex + ey;
@ -268,7 +270,9 @@ fma(double x, double y, double z)
* rounding modes.
*/
fesetround(oround);
adj = r.lo + xy.lo;
/* work around clang bug 8100 */
volatile double vrlo = r.lo;
adj = vrlo + xy.lo;
return (ldexp(r.hi + adj, spread));
}

8
libm/upstream-freebsd/lib/msun/src/s_fmal.c
View File

@ -226,6 +226,8 @@ fmal(long double x, long double y, long double z)
zs = copysignl(LDBL_MIN, zs);
fesetround(FE_TONEAREST);
/* work around clang bug 8100 */
volatile long double vxs = xs;
/*
* Basic approach for round-to-nearest:
@ -235,7 +237,7 @@ fmal(long double x, long double y, long double z)
* adj = xy.lo + r.lo (inexact; low bit is sticky)
* result = r.hi + adj (correctly rounded)
*/
xy = dd_mul(xs, ys);
xy = dd_mul(vxs, ys);
r = dd_add(xy.hi, zs);
spread = ex + ey;
@ -256,7 +258,9 @@ fmal(long double x, long double y, long double z)
* rounding modes.
*/
fesetround(oround);
adj = r.lo + xy.lo;
/* work around clang bug 8100 */
volatile long double vrlo = r.lo;
adj = vrlo + xy.lo;
return (ldexpl(r.hi + adj, spread));
}

7
libm/upstream-freebsd/lib/msun/src/s_log1p.c
View File

@ -96,6 +96,7 @@ Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
static volatile double vzero = 0.0;
double
log1p(double x)
@ -109,7 +110,7 @@ log1p(double x)
k = 1;
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
if(ax>=0x3ff00000) { /* x <= -1.0 */
if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */
if(x==-1.0) return -two54/vzero; /* log1p(-1)=+inf */
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if(ax<0x3e200000) { /* |x| < 2**-29 */
@ -173,3 +174,7 @@ log1p(double x)
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(log1p, log1pl);
#endif

3
libm/upstream-freebsd/lib/msun/src/s_log1pf.c
View File

@ -34,6 +34,7 @@ Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */
static const float zero = 0.0;
static volatile float vzero = 0.0;
float
log1pf(float x)
@ -47,7 +48,7 @@ log1pf(float x)
k = 1;
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
if(ax>=0x3f800000) { /* x <= -1.0 */
if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
if(x==(float)-1.0) return -two25/vzero; /* log1p(-1)=+inf */
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if(ax<0x38000000) { /* |x| < 2**-15 */

6
libm/upstream-freebsd/lib/msun/src/s_nearbyint.c
View File

@ -36,12 +36,16 @@ __FBSDID("$FreeBSD$");
* instead of feclearexcept()/feupdateenv() to restore the environment
* because the only exception defined for rint() is overflow, and
* rounding can't overflow as long as emax >= p.
*
* The volatile keyword is needed below because clang incorrectly assumes
* that rint won't raise any floating-point exceptions. Declaring ret volatile
* is sufficient to trick the compiler into doing the right thing.
*/
#define DECL(type, fn, rint) \
type \
fn(type x) \
{ \
type ret; \
volatile type ret; \
fenv_t env; \
\
fegetenv(&env); \