Skip to content

Commit

Permalink
[BZ #2423]
Browse files Browse the repository at this point in the history 
2006-03-07  Jakub Jelinek  <[email protected]>
	[BZ #2423]
	* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
	round_test, trunc_test): Only run some of the new tests if
	LDBL_MANT_DIG > 100.

2006-03-03  Steven Munroe  <[email protected]>
	    Alan Modra  <[email protected]>

	* sysdeps/powerpc/fpu/fenv_libc.h (__fegetround, __fesetround):
	Define inline implementations.
	* sysdeps/powerpc/fpu/fegetround.c: Use __fegetround.
	* sysdeps/powerpc/fpu/fesetround.c: Use __fesetround.

	* sysdeps/powerpc/fpu/math_ldbl.h: New file.

	[BZ #2423]
	* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
	round_test, trunc_test): Add new tests.
	* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
	(EXTRACT_IBM_EXTENDED_MANTISSA, INSERT_IBM_EXTENDED_MANTISSA):
	Removed, replaced with ...
	(ldbl_extract_mantissa, ldbl_insert_mantissa, ldbl_pack, ldbl_unpack,
	ldbl_canonicalise, ldbl_nearbyint): New functions.
	* sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Replace
	EXTRACT_IBM_EXTENDED_MANTISSA and INSERT_IBM_EXTENDED_MANTISSA
	with ldbl_extract_mantissa and ldbl_insert_mantissa.
	* sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c (__ieee754_rem_pio2l):
	Replace EXTRACT_IBM_EXTENDED_MANTISSA with ldbl_extract_mantissa.
	(ldbl_extract_mantissa, ldbl_insert_mantissa): New inline functions.
	* sysdeps/ieee754/ldbl-128ibm/s_ceill.c (__ceill): Handle rounding
	that spans doubles in IBM long double format.
	* sysdeps/ieee754/ldbl-128ibm/s_floorl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_roundl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_truncl.c: Likewise.
	* sysdeps/powerpc/powerpc64/fpu/s_rintl.S: File removed.
  • Loading branch information
Roland McGrath committed Mar 16, 2006
1 parent 671ca69 commit 5c68d40
Show file tree
Hide file tree
Showing 15 changed files with 1,047 additions and 623 deletions.
39 changes: 39 additions & 0 deletions ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,42 @@
2006-03-07 Jakub Jelinek <[email protected]>

[BZ #2423]
* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
round_test, trunc_test): Only run some of the new tests if
LDBL_MANT_DIG > 100.

2006-03-03 Steven Munroe <[email protected]>
Alan Modra <[email protected]>

* sysdeps/powerpc/fpu/fenv_libc.h (__fegetround, __fesetround):
Define inline implementations.
* sysdeps/powerpc/fpu/fegetround.c: Use __fegetround.
* sysdeps/powerpc/fpu/fesetround.c: Use __fesetround.

* sysdeps/powerpc/fpu/math_ldbl.h: New file.

[BZ #2423]
* math/libm-test.inc [TEST_LDOUBLE] (ceil_test, floor_test, rint_test,
round_test, trunc_test): Add new tests.
* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
(EXTRACT_IBM_EXTENDED_MANTISSA, INSERT_IBM_EXTENDED_MANTISSA):
Removed, replaced with ...
(ldbl_extract_mantissa, ldbl_insert_mantissa, ldbl_pack, ldbl_unpack,
ldbl_canonicalise, ldbl_nearbyint): New functions.
* sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Replace
EXTRACT_IBM_EXTENDED_MANTISSA and INSERT_IBM_EXTENDED_MANTISSA
with ldbl_extract_mantissa and ldbl_insert_mantissa.
* sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c (__ieee754_rem_pio2l):
Replace EXTRACT_IBM_EXTENDED_MANTISSA with ldbl_extract_mantissa.
(ldbl_extract_mantissa, ldbl_insert_mantissa): New inline functions.
* sysdeps/ieee754/ldbl-128ibm/s_ceill.c (__ceill): Handle rounding
that spans doubles in IBM long double format.
* sysdeps/ieee754/ldbl-128ibm/s_floorl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_roundl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_truncl.c: Likewise.
* sysdeps/powerpc/powerpc64/fpu/s_rintl.S: File removed.

2004-12-09 Randolph Chung <[email protected]>

* sysdeps/unix/sysv/linux/kernel-features.h (__ASSUME_UTIMES): Don't
Expand Down
334 changes: 332 additions & 2 deletions math/libm-test.inc
View file

Large diffs are not rendered by default.

8 changes: 4 additions & 4 deletions sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,8 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
/* Make the IBM extended format 105 bit mantissa look like the ieee854 112
bit mantissa so the following operatations will give the correct
result. */
EXTRACT_IBM_EXTENDED_MANTISSA(hx, lx, temp, x);
EXTRACT_IBM_EXTENDED_MANTISSA(hy, ly, temp, y);
ldbl_extract_mantissa(&hx, &lx, &temp, x);
ldbl_extract_mantissa(&hy, &ly, &temp, y);

/* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -1022)
Expand Down Expand Up @@ -127,7 +127,7 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
iy -= 1;
}
if(iy>= -1022) { /* normalize output */
INSERT_IBM_EXTENDED_MANTISSA(x, (sx>>63), iy, hx, lx);
x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
} else { /* subnormal output */
n = -1022 - iy;
if(n<=48) {
Expand All @@ -138,7 +138,7 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
} else {
lx = hx>>(n-64); hx = sx;
}
INSERT_IBM_EXTENDED_MANTISSA(x, (sx>>63), iy, hx, lx);
x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
Expand Down
5 changes: 3 additions & 2 deletions sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,8 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
{
long double z, w, t;
double tx[8];
int64_t exp, n, ix, hx, ixd;
int exp;
int64_t n, ix, hx, ixd;
u_int64_t lx, lxd;

GET_LDOUBLE_WORDS64 (hx, lx, x);
Expand Down Expand Up @@ -243,7 +244,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
stored in a double array. */
/* Make the IBM extended format 105 bit mantissa look like the ieee854 112
bit mantissa so the next operatation will give the correct result. */
EXTRACT_IBM_EXTENDED_MANTISSA (ixd, lxd, exp, x);
ldbl_extract_mantissa (&ixd, &lxd, &exp, x);
exp = exp - 23;
/* This is faster than doing this in floating point, because we
have to convert it to integers anyway and like this we can keep
Expand Down
291 changes: 174 additions & 117 deletions sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,122 +3,179 @@
#endif

#include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
#include <ieee754.h>

static inline void
ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x)
{
/* We have 105 bits of mantissa plus one implicit digit. Since
106 bits are representable we use the first implicit digit for
the number before the decimal point and the second implicit bit
as bit 53 of the mantissa. */
unsigned long long hi, lo;
int ediff;
union ibm_extended_long_double eldbl;
eldbl.d = x;
*exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;

#define EXTRACT_IBM_EXTENDED_MANTISSA(hi64, lo64, expnt, ibm_ext_ldbl) \
do \
{ \
/* We have 105 bits of mantissa plus one implicit digit. Since \
106 bits are representable without the rest using hexadecimal \
digits we use only the implicit digits for the number before \
the decimal point. */ \
unsigned long long hi, lo; \
int ediff; \
union ibm_extended_long_double eldbl; \
eldbl.d = ibm_ext_ldbl; \
expnt = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS; \
\
lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3; \
hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1; \
/* If the lower double is not a denomal or zero then set the hidden \
53rd bit. */ \
if (eldbl.ieee.exponent2 > 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; \
if (ediff > 53) \
lo = lo >> (ediff-53); \
} \
hi |= (1ULL << 52); \
\
if ((eldbl.ieee.negative != eldbl.ieee.negative2) \
&& ((eldbl.ieee.exponent2 != 0) && (lo != 0LL))) \
{ \
hi--; \
lo = (1ULL << 60) - lo; \
if (hi < (1ULL << 52)) \
{ \
/* we have a borrow from the hidden bit, so shift left 1. */ \
hi = (hi << 1) | (lo >> 59); \
lo = 0xfffffffffffffffLL & (lo << 1); \
expnt--; \
} \
} \
lo64 = (hi << 60) | lo; \
hi64 = hi >> 4; \
} \
while (0)
lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3;
hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1;
/* If the lower double is not a denomal or zero then set the hidden
53rd bit. */
if (eldbl.ieee.exponent2 > 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;
if (ediff > 53)
lo = lo >> (ediff-53);
}
hi |= (1ULL << 52);

if ((eldbl.ieee.negative != eldbl.ieee.negative2)
&& ((eldbl.ieee.exponent2 != 0) && (lo != 0LL)))
{
hi--;
lo = (1ULL << 60) - lo;
if (hi < (1ULL << 52))
{
/* we have a borrow from the hidden bit, so shift left 1. */
hi = (hi << 1) | (lo >> 59);
lo = 0xfffffffffffffffLL & (lo << 1);
*exp = *exp - 1;
}
}
*lo64 = (hi << 60) | lo;
*hi64 = hi >> 4;
}

#define INSERT_IBM_EXTENDED_MANTISSA(ibm_ext_ldbl, sign, expnt, hi64, lo64) \
do \
{ \
union ibm_extended_long_double u; \
unsigned long hidden2, lzcount; \
unsigned long long hi, lo; \
\
u.ieee.negative = sign; \
u.ieee.negative2 = sign; \
u.ieee.exponent = expnt + IBM_EXTENDED_LONG_DOUBLE_BIAS; \
u.ieee.exponent2 = expnt-53 + IBM_EXTENDED_LONG_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); \
hidden2 = (lo64 >> 59) & 1ULL; \
/* The high order 53 bits (52 + hidden) go into the upper double */ \
hi = (lo64 >> 60) & ((1ULL << 11) - 1); \
hi |= (hi64 << 4); \
\
if (lo != 0LL) \
{ \
/* hidden2 bit of low double controls rounding of the high double. \
If hidden2 is '1' then round up hi and adjust lo (2nd mantissa) \
plus change the sign of the low double to compensate. */ \
if (hidden2) \
{ \
hi++; \
u.ieee.negative2 = !sign; \
lo = (1ULL << 53) - lo; \
} \
/* The hidden bit of the lo mantissa is zero so we need to \
normalize the it for the low double. Shift it left until the \
hidden bit is '1' then adjust the 2nd exponent accordingly. */ \
\
if (sizeof (lo) == sizeof (long)) \
lzcount = __builtin_clzl (lo); \
else if ((lo >> 32) != 0) \
lzcount = __builtin_clzl ((long) (lo >> 32)); \
else \
lzcount = __builtin_clzl ((long) lo) + 32; \
lzcount = lzcount - 11; \
if (lzcount > 0) \
{ \
int expnt2 = u.ieee.exponent2 - lzcount; \
if (expnt2 >= 1) \
{ \
/* Not denormal. Normalize and set low exponent. */ \
lo = lo << lzcount; \
u.ieee.exponent2 = expnt2; \
} \
else \
{ \
/* Is denormal. */ \
lo = lo << (lzcount + expnt2); \
u.ieee.exponent2 = 0; \
} \
} \
} \
else \
{ \
u.ieee.negative2 = 0; \
u.ieee.exponent2 = 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); \
ibm_ext_ldbl = u.d; \
} \
while (0)
static inline long double
ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
{
union ibm_extended_long_double u;
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;
/* 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);
hidden2 = (lo64 >> 59) & 1ULL;
/* The high order 53 bits (52 + hidden) go into the upper double */
hi = (lo64 >> 60) & ((1ULL << 11) - 1);
hi |= (hi64 << 4);

if (lo != 0LL)
{
/* hidden2 bit of low double controls rounding of the high double.
If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)
plus change the sign of the low double to compensate. */
if (hidden2)
{
hi++;
u.ieee.negative2 = !sign;
lo = (1ULL << 53) - lo;
}
/* The hidden bit of the lo mantissa is zero so we need to
normalize the it for the low double. Shift it left until the
hidden bit is '1' then adjust the 2nd exponent accordingly. */

if (sizeof (lo) == sizeof (long))
lzcount = __builtin_clzl (lo);
else if ((lo >> 32) != 0)
lzcount = __builtin_clzl ((long) (lo >> 32));
else
lzcount = __builtin_clzl ((long) lo) + 32;
lzcount = lzcount - 11;
if (lzcount > 0)
{
int expnt2 = u.ieee.exponent2 - lzcount;
if (expnt2 >= 1)
{
/* Not denormal. Normalize and set low exponent. */
lo = lo << lzcount;
u.ieee.exponent2 = expnt2;
}
else
{
/* Is denormal. */
lo = lo << (lzcount + expnt2);
u.ieee.exponent2 = 0;
}
}
}
else
{
u.ieee.negative2 = 0;
u.ieee.exponent2 = 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;
}

/* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
of long double implemented as double double. */
static inline long double
ldbl_pack (double a, double aa)
{
union ibm_extended_long_double u;
u.dd[0] = a;
u.dd[1] = aa;
return u.d;
}

static inline void
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];
}


/* Convert a finite long double to canonical form.
Does not handle +/-Inf properly. */
static inline void
ldbl_canonicalize (double *a, double *aa)
{
double xh, xl;

xh = *a + *aa;
xl = (*a - xh) + *aa;
*a = xh;
*aa = xl;
}

/* Simple inline nearbyint (double) function .
Only works in the default rounding mode
but is useful in long double rounding functions. */
static inline double
ldbl_nearbyint (double a)
{
double two52 = 0x10000000000000LL;

if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
{
if (__builtin_expect ((a > 0.0), 1))
{
a += two52;
a -= two52;
}
else if (__builtin_expect ((a < 0.0), 1))
{
a = two52 - a;
a = -(a - two52);
}
}
return a;
}
Loading

0 comments on commit 5c68d40

Please sign in to comment.