OpenE2K
/
glibc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

IBM long double mechanical changes to support little-endian 

http://sourceware.org/ml/libc-alpha/2013-07/msg00001.html

This patch starts the process of supporting powerpc64 little-endian
long double in glibc.  IBM long double is an array of two ieee
doubles, so making union ibm_extended_long_double reflect this fact is
the correct way to access fields of the doubles.

	* sysdeps/ieee754/ldbl-128ibm/ieee754.h
	(union ibm_extended_long_double): Define as an array of ieee754_double.
	(IBM_EXTENDED_LONG_DOUBLE_BIAS): Delete.
	* sysdeps/ieee754/ldbl-128ibm/printf_fphex.c: Update all references
	to ibm_extended_long_double and IBM_EXTENDED_LONG_DOUBLE_BIAS.
	* sysdeps/ieee754/ldbl-128ibm/e_exp10l.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_expl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/strtold_l.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c: Likewise.
This commit is contained in:
Alan Modra 2013-08-17 18:12:56 +09:30
parent d3d237560b
commit 9605ca6c08
11 changed files with 161 additions and 184 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
ChangeLog
View File

@ -1,3 +1,19 @@
2013-10-04 Alan Modra <amodra@gmail.com>
* sysdeps/ieee754/ldbl-128ibm/ieee754.h
(union ibm_extended_long_double): Define as an array of ieee754_double.
(IBM_EXTENDED_LONG_DOUBLE_BIAS): Delete.
* sysdeps/ieee754/ldbl-128ibm/printf_fphex.c: Update all references
to ibm_extended_long_double and IBM_EXTENDED_LONG_DOUBLE_BIAS.
* sysdeps/ieee754/ldbl-128ibm/e_exp10l.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_expl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h: Likewise.
* sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/strtold_l.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c: Likewise.
2013-10-03 Joseph Myers <joseph@codesourcery.com>
* locale/programs/locarchive.c (add_locale): Use constant 4096 for

6
sysdeps/ieee754/ldbl-128ibm/e_exp10l.c
View File

@ -36,9 +36,9 @@ __ieee754_exp10l (long double arg)
else if (arg > LDBL_MAX_10_EXP + 1)
return LDBL_MAX * LDBL_MAX;
u.d = arg;
arg_high = u.dd[0];
arg_low = u.dd[1];
u.ld = arg;
arg_high = u.d[0].d;
arg_low = u.d[1].d;
exp_high = arg_high * log10_high;
exp_low = arg_high * log10_low + arg_low * M_LN10l;
return __ieee754_expl (exp_high) * __ieee754_expl (exp_low);

34
sysdeps/ieee754/ldbl-128ibm/e_expl.c
View File

@ -162,39 +162,39 @@ __ieee754_expl (long double x)
x = x + xl;
/* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */
ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
* __expl_table[T_EXPL_RES2 + tval2];
ex2_u.ld = (__expl_table[T_EXPL_RES1 + tval1]
* __expl_table[T_EXPL_RES2 + tval2]);
n_i = (int)n;
/* 'unsafe' is 1 iff n_1 != 0. */
unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
ex2_u.ieee.exponent += n_i >> unsafe;
ex2_u.d[0].ieee.exponent += n_i >> unsafe;
/* Fortunately, there are no subnormal lowpart doubles in
__expl_table, only normal values and zeros.
But after scaling it can be subnormal. */
exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe);
if (ex2_u.ieee.exponent2 == 0)
/* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */;
exponent2 = ex2_u.d[1].ieee.exponent + (n_i >> unsafe);
if (ex2_u.d[1].ieee.exponent == 0)
/* assert ((ex2_u.d[1].ieee.mantissa0|ex2_u.d[1].ieee.mantissa1) == 0) */;
else if (exponent2 > 0)
ex2_u.ieee.exponent2 = exponent2;
ex2_u.d[1].ieee.exponent = exponent2;
else if (exponent2 <= -54)
{
ex2_u.ieee.exponent2 = 0;
ex2_u.ieee.mantissa2 = 0;
ex2_u.ieee.mantissa3 = 0;
ex2_u.d[1].ieee.exponent = 0;
ex2_u.d[1].ieee.mantissa0 = 0;
ex2_u.d[1].ieee.mantissa1 = 0;
}
else
{
static const double
two54 = 1.80143985094819840000e+16, /* 4350000000000000 */
twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */
ex2_u.dd[1] *= two54;
ex2_u.ieee.exponent2 += n_i >> unsafe;
ex2_u.dd[1] *= twom54;
ex2_u.d[1].d *= two54;
ex2_u.d[1].ieee.exponent += n_i >> unsafe;
ex2_u.d[1].d *= twom54;
}
/* Compute scale = 2^n_1. */
scale_u.d = 1.0L;
scale_u.ieee.exponent += n_i - (n_i >> unsafe);
scale_u.ld = 1.0L;
scale_u.d[0].ieee.exponent += n_i - (n_i >> unsafe);
/* Approximate e^x2 - 1, using a seventh-degree polynomial,
with maximum error in [-2^-16-2^-53,2^-16+2^-53]
@ -204,7 +204,7 @@ __ieee754_expl (long double x)
/* Return result. */
fesetenv (&oldenv);
result = x22 * ex2_u.d + ex2_u.d;
result = x22 * ex2_u.ld + ex2_u.ld;
/* Now we can test whether the result is ultimate or if we are unsure.
In the later case we should probably call a mpn based routine to give
@ -238,7 +238,7 @@ __ieee754_expl (long double x)
if (!unsafe)
return result;
else
return result * scale_u.d;
return result * scale_u.ld;
}
/* Exceptional cases: */
else if (isless (x, himark))

43
sysdeps/ieee754/ldbl-128ibm/ieee754.h
View File

@ -179,49 +179,10 @@ union ieee854_long_double
union ibm_extended_long_double
{
long double d;
double dd[2];
/* This is the IBM extended format long double. */
struct
{ /* Big endian. There is no other. */
unsigned int negative:1;
unsigned int exponent:11;
/* Together Mantissa0-3 comprise the mantissa. */
unsigned int mantissa0:20;
unsigned int mantissa1:32;
unsigned int negative2:1;
unsigned int exponent2:11;
/* There is an implied 1 here? */
/* Together these comprise the mantissa. */
unsigned int mantissa2:20;
unsigned int mantissa3:32;
} ieee;
/* This format makes it easier to see if a NaN is a signalling NaN. */
struct
{ /* Big endian. There is no other. */
unsigned int negative:1;
unsigned int exponent:11;
unsigned int quiet_nan:1;
/* Together Mantissa0-3 comprise the mantissa. */
unsigned int mantissa0:19;
unsigned int mantissa1:32;
unsigned int negative2:1;
unsigned int exponent2:11;
/* There is an implied 1 here? */
/* Together these comprise the mantissa. */
unsigned int mantissa2:20;
unsigned int mantissa3:32;
} ieee_nan;
long double ld;
union ieee754_double d[2];
};
#define IBM_EXTENDED_LONG_DOUBLE_BIAS 0x3ff /* Added to exponent. */
__END_DECLS
#endif /* ieee754.h */

20
sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c
View File

@ -36,22 +36,22 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
union ibm_extended_long_double u;
unsigned long long hi, lo;
int ediff;
u.d = value;
u.ld = value;
*is_neg = u.ieee.negative;
*expt = (int) u.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;
*is_neg = u.d[0].ieee.negative;
*expt = (int) u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
lo = ((long long) u.ieee.mantissa2 << 32) | u.ieee.mantissa3;
hi = ((long long) u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
lo = ((long long) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
hi = ((long long) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
/* If the lower double is not a denomal or zero then set the hidden
53rd bit. */
if (u.ieee.exponent2 > 0)
if (u.d[1].ieee.exponent > 0)
{
lo |= 1LL << 52;
/* The lower double is normalized separately from the upper. We may
need to adjust the lower manitissa to reflect this. */
ediff = u.ieee.exponent - u.ieee.exponent2;
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;
if (ediff > 53)
lo = lo >> (ediff-53);
}
@ -59,8 +59,8 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
difference between the long double and the rounded high double
value. This is indicated by a differnce between the signs of the
high and low doubles. */
if ((u.ieee.negative != u.ieee.negative2)
&& ((u.ieee.exponent2 != 0) && (lo != 0L)))
if ((u.d[0].ieee.negative != u.d[1].ieee.negative)
&& ((u.d[1].ieee.exponent != 0) && (lo != 0L)))
{
lo = (1ULL << 53) - lo;
if (hi == 0LL)
@ -92,7 +92,7 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
- (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
if (u.ieee.exponent == 0)
if (u.d[0].ieee.exponent == 0)
{
/* A biased exponent of zero is a special case.
Either it is a zero or it is a denormal number. */

60
sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
View File

@ -15,28 +15,28 @@ ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x)
as bit 53 of the mantissa. */
uint64_t hi, lo;
int ediff;
union ibm_extended_long_double eldbl;
eldbl.d = x;
*exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;
union ibm_extended_long_double u;
u.ld = x;
*exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
lo = ((int64_t)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3;
hi = ((int64_t)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1;
lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
/* If the lower double is not a denomal or zero then set the hidden
53rd bit. */
if (eldbl.ieee.exponent2 > 0x001)
if (u.d[1].ieee.exponent > 0x001)
{
lo |= (1ULL << 52);
lo = lo << 7; /* pre-shift lo to match ieee854. */
/* The lower double is normalized separately from the upper. We
may need to adjust the lower manitissa to reflect this. */
ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2;
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;
if (ediff > 53)
lo = lo >> (ediff-53);
hi |= (1ULL << 52);
}
if ((eldbl.ieee.negative != eldbl.ieee.negative2)
&& ((eldbl.ieee.exponent2 != 0) && (lo != 0LL)))
if ((u.d[0].ieee.negative != u.d[1].ieee.negative)
&& ((u.d[1].ieee.exponent != 0) && (lo != 0LL)))
{
hi--;
lo = (1ULL << 60) - lo;
@ -59,10 +59,10 @@ ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
unsigned long hidden2, lzcount;
unsigned long long hi, lo;
u.ieee.negative = sign;
u.ieee.negative2 = sign;
u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS;
u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
u.d[0].ieee.negative = sign;
u.d[1].ieee.negative = sign;
u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS;
u.d[1].ieee.exponent = exp-53 + IEEE754_DOUBLE_BIAS;
/* Expect 113 bits (112 bits + hidden) right justified in two longs.
The low order 53 bits (52 + hidden) go into the lower double */
lo = (lo64 >> 7)& ((1ULL << 53) - 1);
@ -79,7 +79,7 @@ ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
if (hidden2)
{
hi++;
u.ieee.negative2 = !sign;
u.d[1].ieee.negative = !sign;
lo = (1ULL << 53) - lo;
}
/* The hidden bit of the lo mantissa is zero so we need to
@ -95,32 +95,32 @@ ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
lzcount = lzcount - 11;
if (lzcount > 0)
{
int expnt2 = u.ieee.exponent2 - lzcount;
int expnt2 = u.d[1].ieee.exponent - lzcount;
if (expnt2 >= 1)
{
/* Not denormal. Normalize and set low exponent. */
lo = lo << lzcount;
u.ieee.exponent2 = expnt2;
u.d[1].ieee.exponent = expnt2;
}
else
{
/* Is denormal. */
lo = lo << (lzcount + expnt2);
u.ieee.exponent2 = 0;
u.d[1].ieee.exponent = 0;
}
}
}
else
{
u.ieee.negative2 = 0;
u.ieee.exponent2 = 0;
u.d[1].ieee.negative = 0;
u.d[1].ieee.exponent = 0;
}
u.ieee.mantissa3 = lo & ((1ULL << 32) - 1);
u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1);
u.ieee.mantissa1 = hi & ((1ULL << 32) - 1);
u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);
return u.d;
u.d[1].ieee.mantissa1 = lo & ((1ULL << 32) - 1);
u.d[1].ieee.mantissa0 = (lo >> 32) & ((1ULL << 20) - 1);
u.d[0].ieee.mantissa1 = hi & ((1ULL << 32) - 1);
u.d[0].ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);
return u.ld;
}
/* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
@ -129,18 +129,18 @@ static inline long double
default_ldbl_pack (double a, double aa)
{
union ibm_extended_long_double u;
u.dd[0] = a;
u.dd[1] = aa;
return u.d;
u.d[0].d = a;
u.d[1].d = aa;
return u.ld;
}
static inline void
default_ldbl_unpack (long double l, double *a, double *aa)
{
union ibm_extended_long_double u;
u.d = l;
*a = u.dd[0];
*aa = u.dd[1];
u.ld = l;
*a = u.d[0].d;
*aa = u.d[1].d;
}
#ifndef ldbl_pack

36
sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
View File

@ -33,11 +33,11 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
unsigned long long hi, lo;
int exponent2;
u.ieee.negative = sign;
u.ieee.negative2 = sign;
u.ieee.exponent = expt + IBM_EXTENDED_LONG_DOUBLE_BIAS;
u.ieee.exponent2 = 0;
exponent2 = expt - 53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
u.d[0].ieee.negative = sign;
u.d[1].ieee.negative = sign;
u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
u.d[1].ieee.exponent = 0;
exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS;
#if BITS_PER_MP_LIMB == 32
/* The low order 53 bits (52 + hidden) go into the lower double */
@ -73,15 +73,15 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
else
lzcount = lzcount + 42;
if (lzcount > u.ieee.exponent)
if (lzcount > u.d[0].ieee.exponent)
{
lzcount = u.ieee.exponent;
u.ieee.exponent = 0;
lzcount = u.d[0].ieee.exponent;
u.d[0].ieee.exponent = 0;
exponent2 -= lzcount;
}
else
{
u.ieee.exponent -= (lzcount - 1);
u.d[0].ieee.exponent -= (lzcount - 1);
exponent2 -= (lzcount - 1);
}
@ -111,9 +111,9 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
{
if ((hi & (1LL << 53)) != 0)
hi -= 1LL << 52;
u.ieee.exponent++;
u.d[0].ieee.exponent++;
}
u.ieee.negative2 = !sign;
u.d[1].ieee.negative = !sign;
lo = (1LL << 53) - lo;
}
@ -134,17 +134,17 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
exponent2 = exponent2 - lzcount;
}
if (exponent2 > 0)
u.ieee.exponent2 = exponent2;
u.d[1].ieee.exponent = exponent2;
else
lo >>= 1 - exponent2;
}
else
u.ieee.negative2 = 0;
u.d[1].ieee.negative = 0;
u.ieee.mantissa3 = lo & 0xffffffffLL;
u.ieee.mantissa2 = (lo >> 32) & 0xfffff;
u.ieee.mantissa1 = hi & 0xffffffffLL;
u.ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1);
u.d[1].ieee.mantissa1 = lo & 0xffffffffLL;
u.d[1].ieee.mantissa0 = (lo >> 32) & 0xfffff;
u.d[0].ieee.mantissa1 = hi & 0xffffffffLL;
u.d[0].ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1);
return u.d;
return u.ld;
}

32
sysdeps/ieee754/ldbl-128ibm/printf_fphex.c
View File

@ -26,31 +26,31 @@ do { \
unsigned long long int num0, num1; \
unsigned long long hi, lo; \
int ediff; \
union ibm_extended_long_double eldbl; \
eldbl.d = fpnum.ldbl.d; \
union ibm_extended_long_double u; \
u.ld = fpnum.ldbl.d; \
\
assert (sizeof (long double) == 16); \
\
lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3; \
hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1; \
lo = ((long long)u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1; \
hi = ((long long)u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1; \
lo <<= 7; /* pre-shift lo to match ieee854. */ \
/* If the lower double is not a denomal or zero then set the hidden \
53rd bit. */ \
if (eldbl.ieee.exponent2 != 0) \
if (u.d[1].ieee.exponent != 0) \
lo |= (1ULL << (52 + 7)); \
else \
lo <<= 1; \
/* The lower double is normalized separately from the upper. We \
may need to adjust the lower manitissa to reflect this. */ \
ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2; \
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent; \
if (ediff > 53 + 63) \
lo = 0; \
else if (ediff > 53) \
lo = lo >> (ediff - 53); \
else if (eldbl.ieee.exponent2 == 0 && ediff < 53) \
else if (u.d[1].ieee.exponent == 0 && ediff < 53) \
lo = lo << (53 - ediff); \
if (eldbl.ieee.negative != eldbl.ieee.negative2 \
&& (eldbl.ieee.exponent2 != 0 || lo != 0L)) \
if (u.d[0].ieee.negative != u.d[1].ieee.negative \
&& (u.d[1].ieee.exponent != 0 || lo != 0L)) \
{ \
lo = (1ULL << 60) - lo; \
if (hi == 0L) \
@ -58,7 +58,7 @@ do { \
/* we have a borrow from the hidden bit, so shift left 1. */ \
hi = 0xffffffffffffeLL | (lo >> 59); \
lo = 0xfffffffffffffffLL & (lo << 1); \
eldbl.ieee.exponent--; \
u.d[0].ieee.exponent--; \
} \
else \
hi--; \
@ -109,9 +109,9 @@ do { \
*--wnumstr = L'0'; \
} \
\
leading = eldbl.ieee.exponent == 0 ? '0' : '1'; \
leading = u.d[0].ieee.exponent == 0 ? '0' : '1'; \
\
exponent = eldbl.ieee.exponent; \
exponent = u.d[0].ieee.exponent; \
\
if (exponent == 0) \
{ \
@ -121,18 +121,18 @@ do { \
{ \
/* This is a denormalized number. */ \
expnegative = 1; \
exponent = IBM_EXTENDED_LONG_DOUBLE_BIAS - 1; \
exponent = IEEE754_DOUBLE_BIAS - 1; \
} \
} \
else if (exponent >= IBM_EXTENDED_LONG_DOUBLE_BIAS) \
else if (exponent >= IEEE754_DOUBLE_BIAS) \
{ \
expnegative = 0; \
exponent -= IBM_EXTENDED_LONG_DOUBLE_BIAS; \
exponent -= IEEE754_DOUBLE_BIAS; \
} \
else \
{ \
expnegative = 1; \
exponent = -(exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS); \
exponent = -(exponent - IEEE754_DOUBLE_BIAS); \
} \
} while (0)

62
sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c
View File

@ -34,11 +34,11 @@ __nearbyintl (long double x)
fenv_t env;
static const long double TWO52 = 4503599627370496.0L;
union ibm_extended_long_double u;
u.d = x;
u.ld = x;
if (fabs (u.dd[0]) < TWO52)
if (fabs (u.d[0].d) < TWO52)
{
double high = u.dd[0];
double high = u.d[0].d;
feholdexcept (&env);
if (high > 0.0)
{
@ -52,13 +52,13 @@ __nearbyintl (long double x)
high += TWO52;
if (high == 0.0) high = -0.0;
}
u.dd[0] = high;
u.dd[1] = 0.0;
math_force_eval (u.dd[0]);
math_force_eval (u.dd[1]);
u.d[0].d = high;
u.d[1].d = 0.0;
math_force_eval (u.d[0]);
math_force_eval (u.d[1]);
fesetenv (&env);
}
else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0)
else if (fabs (u.d[1].d) < TWO52 && u.d[1].d != 0.0)
{
double high, low, tau;
/* In this case we have to round the low double and handle any
@ -67,57 +67,57 @@ __nearbyintl (long double x)
may already be rounded and the low double may have the
opposite sign to compensate. */
feholdexcept (&env);
if (u.dd[0] > 0.0)
if (u.d[0].d > 0.0)
{
if (u.dd[1] > 0.0)
if (u.d[1].d > 0.0)
{
/* If the high/low doubles are the same sign then simply
round the low double. */
high = u.dd[0];
low = u.dd[1];
high = u.d[0].d;
low = u.d[1].d;
}
else if (u.dd[1] < 0.0)
else if (u.d[1].d < 0.0)
{
/* Else the high double is pre rounded and we need to
adjust for that. */
tau = __nextafter (u.dd[0], 0.0);
tau = (u.dd[0] - tau) * 2.0;
high = u.dd[0] - tau;
low = u.dd[1] + tau;
tau = __nextafter (u.d[0].d, 0.0);
tau = (u.d[0].d - tau) * 2.0;
high = u.d[0].d - tau;
low = u.d[1].d + tau;
}
low += TWO52;
low -= TWO52;
}
else if (u.dd[0] < 0.0)
else if (u.d[0].d < 0.0)
{
if (u.dd[1] < 0.0)
if (u.d[1].d < 0.0)
{
/* If the high/low doubles are the same sign then simply
round the low double. */
high = u.dd[0];
low = u.dd[1];
high = u.d[0].d;
low = u.d[1].d;
}
else if (u.dd[1] > 0.0)
else if (u.d[1].d > 0.0)
{
/* Else the high double is pre rounded and we need to
adjust for that. */
tau = __nextafter (u.dd[0], 0.0);
tau = (u.dd[0] - tau) * 2.0;
high = u.dd[0] - tau;
low = u.dd[1] + tau;
tau = __nextafter (u.d[0].d, 0.0);
tau = (u.d[0].d - tau) * 2.0;
high = u.d[0].d - tau;
low = u.d[1].d + tau;
}
low = TWO52 - low;
low = -(low - TWO52);
}
u.dd[0] = high + low;
u.dd[1] = high - u.dd[0] + low;
math_force_eval (u.dd[0]);
math_force_eval (u.dd[1]);
u.d[0].d = high + low;
u.d[1].d = high - u.d[0].d + low;
math_force_eval (u.d[0]);
math_force_eval (u.d[1]);
fesetenv (&env);
}
return u.d;
return u.ld;
}
long_double_symbol (libm, __nearbyintl, nearbyintl);

10
sysdeps/ieee754/ldbl-128ibm/strtold_l.c
View File

@ -43,11 +43,11 @@ libc_hidden_proto (STRTOF)
#define FLOAT_HUGE_VAL HUGE_VALL
# define SET_MANTISSA(flt, mant) \
do { union ibm_extended_long_double u; \
u.d = (flt); \
u.ieee_nan.mantissa0 = (mant) >> 32; \
u.ieee_nan.mantissa1 = (mant); \
if ((u.ieee.mantissa0 | u.ieee.mantissa1) != 0) \
(flt) = u.d; \
u.ld = (flt); \
u.d[0].ieee_nan.mantissa0 = (mant) >> 32; \
u.d[0].ieee_nan.mantissa1 = (mant); \
if ((u.d[0].ieee.mantissa0 | u.d[0].ieee.mantissa1) != 0) \
(flt) = u.ld; \
} while (0)
#include <strtod_l.c>

26
sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c
View File

@ -89,23 +89,23 @@ __x2y2m1l (long double x, long double y)
double vals[12];
SET_RESTORE_ROUND (FE_TONEAREST);
union ibm_extended_long_double xu, yu;
xu.d = x;
yu.d = y;
if (fabs (xu.dd[1]) < 0x1p-500)
xu.dd[1] = 0.0;
if (fabs (yu.dd[1]) < 0x1p-500)
yu.dd[1] = 0.0;
mul_split (&vals[1], &vals[0], xu.dd[0], xu.dd[0]);
mul_split (&vals[3], &vals[2], xu.dd[0], xu.dd[1]);
xu.ld = x;
yu.ld = y;
if (fabs (xu.d[1].d) < 0x1p-500)
xu.d[1].d = 0.0;
if (fabs (yu.d[1].d) < 0x1p-500)
yu.d[1].d = 0.0;
mul_split (&vals[1], &vals[0], xu.d[0].d, xu.d[0].d);
mul_split (&vals[3], &vals[2], xu.d[0].d, xu.d[1].d);
vals[2] *= 2.0;
vals[3] *= 2.0;
mul_split (&vals[5], &vals[4], xu.dd[1], xu.dd[1]);
mul_split (&vals[7], &vals[6], yu.dd[0], yu.dd[0]);
mul_split (&vals[9], &vals[8], yu.dd[0], yu.dd[1]);
mul_split (&vals[5], &vals[4], xu.d[1].d, xu.d[1].d);
mul_split (&vals[7], &vals[6], yu.d[0].d, yu.d[0].d);
mul_split (&vals[9], &vals[8], yu.d[0].d, yu.d[1].d);
vals[8] *= 2.0;
vals[9] *= 2.0;
mul_split (&vals[11], &vals[10], yu.dd[1], yu.dd[1]);
if (xu.dd[0] >= 0.75)
mul_split (&vals[11], &vals[10], yu.d[1].d, yu.d[1].d);
if (xu.d[0].d >= 0.75)
vals[1] -= 1.0;
else
{