From: Alan Modra Date: Sat, 17 Aug 2013 09:41:11 +0000 (+0930) Subject: A little more precise %Lf for IBM long double X-Git-Url: http://git.ipfire.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c325a11f7fa409423514de4af55ba63cc559cbf0;p=thirdparty%2Fglibc.git A little more precise %Lf for IBM long double http://sourceware.org/ml/libc-alpha/2013-08/msg00106.html IBM long double has variable precision and we often have values with 107 bits during the normal course of calculations. It's possible to have an IBM long double with 2k bits of precision, and while some might argue we ought to print to the full precision, I'm not inclined to try to implement that. So this patch gives the mpn values returned from ldbl2mpn() an extra 11 bits of precision. To do that it's necessary to inform the caller of __mpn_extract_long_double() exactly how many bits are returned rather than assuming a fixed LDBL_MANT_DIG bits. I did that indirectly by returning the number of zero bits in the last mp_limb, which turns out to be convenient both in __mpn_extract_long_double(), and it's caller. Given the extra parameter it then becomes possible to omit shifting in __mpn_extract_long_double(), since the caller does that anyway. In the following, note that 1 - IEEE854_LONG_DOUBLE_BIAS is equal to LDBL_MIN_EXP - 1. I prefer the former since we already use IEEE854_LONG_DOUBLE_BIAS in the exponent calculation for normalised values, and it reinforces the fact that denormals are treated as if their unbiased exponent was 1. [BZ #5268] * include/gmp.h (__mpn_extract_double, __mpn_extract_long_double): Update prototypes. * stdio-common/printf_fp.c: Likewise. (__printf_fp): Use returned zero_bits. * stdlib/dbl2mpn.c (__mpn_extract_double): Update funtion arguments. * sysdeps/i386/ldbl2mpn.c (__mpn_extract_long_double): Don't perform shifts for denormals here. Instead return leading zero bit count and adjust return value of function. * sysdeps/ieee754/dbl-64/dbl2mpn.c (__mpn_extract_long_double): Likewise. * sysdeps/ieee754/ldbl-128/ldbl2mpn.c (__mpn_extract_long_double): Likewise. * sysdeps/ieee754/ldbl-96/ldbl2mpn.c (__mpn_extract_long_double): Likewise. * sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c (__mpn_extract_long_double): Likewise. Return an extra 11 bits of precision. --- diff --git a/ChangeLog b/ChangeLog index 0bea4a90c42..b18cdaacd89 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,22 @@ +2013-11-16 Alan Modra + + [BZ #5268] + * include/gmp.h (__mpn_extract_double, __mpn_extract_long_double): + Update prototypes. + * stdio-common/printf_fp.c: Likewise. + (__printf_fp): Use returned zero_bits. + * stdlib/dbl2mpn.c (__mpn_extract_double): Update funtion arguments. + * sysdeps/i386/ldbl2mpn.c (__mpn_extract_long_double): Don't perform + shifts for denormals here. Instead return leading zero bit count + and adjust return value of function. + * sysdeps/ieee754/dbl-64/dbl2mpn.c (__mpn_extract_long_double): Likewise. + * sysdeps/ieee754/ldbl-128/ldbl2mpn.c (__mpn_extract_long_double): + Likewise. + * sysdeps/ieee754/ldbl-96/ldbl2mpn.c (__mpn_extract_long_double): + Likewise. + * sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c (__mpn_extract_long_double): + Likewise. Return an extra 11 bits of precision. + 2013-11-16 Alan Modra * Makerules (abilist): Define and use var in abilist rules. diff --git a/include/gmp.h b/include/gmp.h index b74167097d9..bcc9360c16a 100644 --- a/include/gmp.h +++ b/include/gmp.h @@ -8,12 +8,12 @@ /* Now define the internal interfaces. */ extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, - double value); + int *expt, int *zero_bits, + int *is_neg, double value); extern mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, - long double value); + int *expt, int *zero_bits, + int *is_neg, long double value); extern float __mpn_construct_float (mp_srcptr frac_ptr, int expt, int sign); diff --git a/stdio-common/printf_fp.c b/stdio-common/printf_fp.c index 2b93e6c57ac..c0a7a82454f 100644 --- a/stdio-common/printf_fp.c +++ b/stdio-common/printf_fp.c @@ -134,11 +134,11 @@ (u##size > v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) >= 0)) extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, - double value); + int *expt, int *zero_bits, + int *is_neg, double value); extern mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, - long double value); + int *expt, int *zero_bits, + int *is_neg, long double value); extern unsigned int __guess_grouping (unsigned int intdig_max, const char *grouping); @@ -360,12 +360,13 @@ ___printf_fp (FILE *fp, } else { + int zero_bits; fracsize = __mpn_extract_long_double (fp_input, (sizeof (fp_input) / sizeof (fp_input[0])), - &exponent, &is_neg, + &exponent, &zero_bits, &is_neg, fpnum.ldbl); - to_shift = 1 + fracsize * BITS_PER_MP_LIMB - LDBL_MANT_DIG; + to_shift = 1 + zero_bits; } } else @@ -406,11 +407,13 @@ ___printf_fp (FILE *fp, } else { + int zero_bits; fracsize = __mpn_extract_double (fp_input, (sizeof (fp_input) / sizeof (fp_input[0])), - &exponent, &is_neg, fpnum.dbl); - to_shift = 1 + fracsize * BITS_PER_MP_LIMB - DBL_MANT_DIG; + &exponent, &zero_bits, &is_neg, + fpnum.dbl); + to_shift = 1 + zero_bits; } } diff --git a/stdlib/dbl2mpn.c b/stdlib/dbl2mpn.c index 429e20aa7bb..dd835f36ab1 100644 --- a/stdlib/dbl2mpn.c +++ b/stdlib/dbl2mpn.c @@ -24,7 +24,7 @@ mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, + int *expt, int *zero_bits, int *is_neg, double value) { #error "__mpn_extract_double is not implemented for this floating point format" diff --git a/sysdeps/i386/ldbl2mpn.c b/sysdeps/i386/ldbl2mpn.c index c7b322b4528..968aea5ec16 100644 --- a/sysdeps/i386/ldbl2mpn.c +++ b/sysdeps/i386/ldbl2mpn.c @@ -29,7 +29,7 @@ mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, + int *expt, int *zero_bits, int *is_neg, long double value) { union ieee854_long_double u; @@ -51,18 +51,26 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" #endif +/* The format does not fill the last limb. There are some zeros. */ +#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - LDBL_MANT_DIG) + *zero_bits = NUM_LEADING_ZEROS; + if (u.ieee.exponent == 0) { /* A biased exponent of zero is a special case. Either it is a zero or it is a denormal number. */ if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2. */ - /* It's zero. */ - *expt = 0; + { + /* It's zero. */ + *expt = 0; + return 1; + } else { /* It is a denormal number, meaning it has no implicit leading one bit, and its exponent is in fact the format minimum. */ int cnt; + int exp = 1 - IEEE854_LONG_DOUBLE_BIAS; /* One problem with Intel's 80-bit format is that the explicit leading one in the normalized representation has to be zero @@ -74,24 +82,17 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, if (res_ptr[N - 1] != 0) { count_leading_zeros (cnt, res_ptr[N - 1]); - if (cnt != 0) - { -#if N == 2 - res_ptr[N - 1] = res_ptr[N - 1] << cnt - | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); - res_ptr[0] <<= cnt; -#else - res_ptr[N - 1] <<= cnt; -#endif - } - *expt = LDBL_MIN_EXP - 1 - cnt; + *zero_bits = cnt; + *expt = exp + NUM_LEADING_ZEROS - cnt; + return N; } else if (res_ptr[0] != 0) { count_leading_zeros (cnt, res_ptr[0]); - res_ptr[N - 1] = res_ptr[0] << cnt; - res_ptr[0] = 0; - *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt; + exp -= BITS_PER_MP_LIMB; + *zero_bits = cnt; + *expt = exp + NUM_LEADING_ZEROS - cnt; + return 1; } else { @@ -104,7 +105,7 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, #else res_ptr[0] = 0x8000000000000000ul; #endif - *expt = LDBL_MIN_EXP - 1; + *expt = 1 - IEEE854_LONG_DOUBLE_BIAS; } } } diff --git a/sysdeps/ieee754/dbl-64/dbl2mpn.c b/sysdeps/ieee754/dbl-64/dbl2mpn.c index f2294de5aaa..cae062a7e04 100644 --- a/sysdeps/ieee754/dbl-64/dbl2mpn.c +++ b/sysdeps/ieee754/dbl-64/dbl2mpn.c @@ -28,7 +28,7 @@ mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, + int *expt, int *zero_bits, int *is_neg, double value) { union ieee754_double u; @@ -49,9 +49,10 @@ __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, #else #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" #endif + /* The format does not fill the last limb. There are some zeros. */ -#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ - - (DBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) +#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - DBL_MANT_DIG) + *zero_bits = NUM_LEADING_ZEROS; if (u.ieee.exponent == 0) { @@ -65,37 +66,20 @@ __mpn_extract_double (mp_ptr res_ptr, mp_size_t size, /* It is a denormal number, meaning it has no implicit leading one bit, and its exponent is in fact the format minimum. */ int cnt; + int exp = 1 - IEEE754_DOUBLE_BIAS; + int n = N; if (res_ptr[N - 1] != 0) - { - count_leading_zeros (cnt, res_ptr[N - 1]); - cnt -= NUM_LEADING_ZEROS; -#if N == 2 - res_ptr[N - 1] = res_ptr[1] << cnt - | (N - 1) - * (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); - res_ptr[0] <<= cnt; -#else - res_ptr[N - 1] <<= cnt; -#endif - *expt = DBL_MIN_EXP - 1 - cnt; - } + count_leading_zeros (cnt, res_ptr[N - 1]); else { count_leading_zeros (cnt, res_ptr[0]); - if (cnt >= NUM_LEADING_ZEROS) - { - res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS); - res_ptr[0] = 0; - } - else - { - res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt); - res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt); - } - *expt = DBL_MIN_EXP - 1 - - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt; + exp -= BITS_PER_MP_LIMB; + n = 1; } + *zero_bits = cnt; + *expt = exp + NUM_LEADING_ZEROS - cnt; + return n; } } else diff --git a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c index 64f98a59d2b..8e2145bb323 100644 --- a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c @@ -30,7 +30,7 @@ mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, + int *expt, int *zero_bits, int *is_neg, long double value) { union ieee854_long_double u; @@ -54,9 +54,10 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, #else #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" #endif + /* The format does not fill the last limb. There are some zeros. */ -#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ - - (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) +#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - LDBL_MANT_DIG) + *zero_bits = NUM_LEADING_ZEROS; if (u.ieee.exponent == 0) { @@ -64,70 +65,42 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, Either it is a zero or it is a denormal number. */ if (res_ptr[0] == 0 && res_ptr[1] == 0 && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4. */ - /* It's zero. */ - *expt = 0; + { + /* It's zero. */ + *expt = 0; + return 1; + } else { /* It is a denormal number, meaning it has no implicit leading one bit, and its exponent is in fact the format minimum. */ int cnt; + int exp = 1 - IEEE854_LONG_DOUBLE_BIAS; + int n = N; #if N == 2 if (res_ptr[N - 1] != 0) - { - count_leading_zeros (cnt, res_ptr[N - 1]); - cnt -= NUM_LEADING_ZEROS; - res_ptr[N - 1] = res_ptr[N - 1] << cnt - | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); - res_ptr[0] <<= cnt; - *expt = LDBL_MIN_EXP - 1 - cnt; - } + count_leading_zeros (cnt, res_ptr[N - 1]); else { count_leading_zeros (cnt, res_ptr[0]); - if (cnt >= NUM_LEADING_ZEROS) - { - res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS); - res_ptr[0] = 0; - } - else - { - res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt); - res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt); - } - *expt = LDBL_MIN_EXP - 1 - - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt; + exp -= BITS_PER_MP_LIMB; + n = 1; } #else - int j, k, l; + int j; for (j = N - 1; j > 0; j--) if (res_ptr[j] != 0) break; count_leading_zeros (cnt, res_ptr[j]); - cnt -= NUM_LEADING_ZEROS; - l = N - 1 - j; - if (cnt < 0) - { - cnt += BITS_PER_MP_LIMB; - l--; - } - if (!cnt) - for (k = N - 1; k >= l; k--) - res_ptr[k] = res_ptr[k-l]; - else - { - for (k = N - 1; k > l; k--) - res_ptr[k] = res_ptr[k-l] << cnt - | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt); - res_ptr[k--] = res_ptr[0] << cnt; - } - - for (; k >= 0; k--) - res_ptr[k] = 0; - *expt = LDBL_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt; + exp -= (N - 1 - j) * BITS_PER_MP_LIMB; + n = j + 1; #endif + *zero_bits = cnt; + *expt = exp + NUM_LEADING_ZEROS - cnt; + return n; } } else diff --git a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c index e46fde74fcf..14f9ed36f18 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c @@ -30,157 +30,121 @@ mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, + int *expt, int *zero_bits, int *is_neg, long double value) { union ibm_extended_long_double u; - unsigned long long hi, lo; + uint64_t hi, lo; int ediff; u.ld = value; *is_neg = u.d[0].ieee.negative; *expt = (int) u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS; +#define NUM_LEADING_ZEROS (128 - (LDBL_MANT_DIG + 11)) + *zero_bits = NUM_LEADING_ZEROS; - 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 = ((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 denormal or zero then set the hidden - 53rd bit. */ - if (u.d[1].ieee.exponent != 0) - lo |= 1ULL << 52; - else - lo = lo << 1; - - /* The lower double is normalized separately from the upper. We may - need to adjust the lower manitissa to reflect this. */ - ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; - if (ediff > 0) + if (u.d[0].ieee.exponent != 0) { - if (ediff < 64) - lo = lo >> ediff; + /* If the high double is not a denormal or zero then set the hidden + 53rd bit. */ + hi |= (uint64_t) 1 << 52; + + /* If the lower double is not a denormal or zero then set the hidden + 53rd bit. */ + if (u.d[1].ieee.exponent != 0) + lo |= (uint64_t) 1 << 52; else - lo = 0; - } - else if (ediff < 0) - lo = lo << -ediff; - - /* The high double may be rounded and the low double reflects the - 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.d[0].ieee.negative != u.d[1].ieee.negative - && lo != 0) - { - lo = (1ULL << 53) - lo; - if (hi == 0) + lo = lo << 1; + + /* We currently only have 53 bits in lo. Gain a few more bits + of precision. */ + lo = lo << 11; + + /* The lower double is normalized separately from the upper. We may + need to adjust the lower manitissa to reflect this. */ + ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; + if (ediff > 0) { - /* we have a borrow from the hidden bit, so shift left 1. */ - hi = 0x0ffffffffffffeLL | (lo >> 51); - lo = 0x1fffffffffffffLL & (lo << 1); - (*expt)--; + if (ediff < 64) + lo = lo >> ediff; + else + lo = 0; } - else - hi--; - } + else if (ediff < 0) + lo = lo << -ediff; + + /* The high double may be rounded and the low double reflects the + 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.d[0].ieee.negative != u.d[1].ieee.negative + && lo != 0) + { + hi--; + lo = -lo; + if (hi < (uint64_t) 1 << 52) + { + /* We have a borrow from the hidden bit, so shift left 1. */ + hi = (hi << 1) | (lo >> 63); + lo = lo << 1; + (*expt)--; + } + } + #if BITS_PER_MP_LIMB == 32 - /* Combine the mantissas to be contiguous. */ - res_ptr[0] = lo; - res_ptr[1] = (hi << (53 - 32)) | (lo >> 32); - res_ptr[2] = hi >> 11; - res_ptr[3] = hi >> (32 + 11); - #define N 4 + res_ptr[0] = lo; + res_ptr[1] = lo >> 32; + res_ptr[2] = hi; + res_ptr[3] = hi >> 32; + return 4; #elif BITS_PER_MP_LIMB == 64 - /* Combine the two mantissas to be contiguous. */ - res_ptr[0] = (hi << 53) | lo; - res_ptr[1] = hi >> 11; - #define N 2 + res_ptr[0] = lo; + res_ptr[1] = hi; + return 2; #else - #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" +# error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" #endif -/* The format does not fill the last limb. There are some zeros. */ -#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ - - (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) + } - if (u.d[0].ieee.exponent == 0) + /* The high double is a denormal or zero. The low double must + be zero. A denormal is interpreted as having a biased + exponent of 1. */ + res_ptr[0] = hi; +#if BITS_PER_MP_LIMB == 32 + res_ptr[1] = hi >> 32; +#endif + if (hi == 0) { - /* A biased exponent of zero is a special case. - Either it is a zero or it is a denormal number. */ - if (res_ptr[0] == 0 && res_ptr[1] == 0 - && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4. */ - /* It's zero. */ - *expt = 0; + /* It's zero. */ + *expt = 0; + return 1; + } + else + { + int cnt; + int exp = 1 - IEEE754_DOUBLE_BIAS; + int n = 1; + +#if BITS_PER_MP_LIMB == 32 + if (res_ptr[1] != 0) + { + count_leading_zeros (cnt, res_ptr[1]); + n = 2; + } else { - /* It is a denormal number, meaning it has no implicit leading - one bit, and its exponent is in fact the format minimum. We - use DBL_MIN_EXP instead of LDBL_MIN_EXP below because the - latter describes the properties of both parts together, but - the exponent is computed from the high part only. */ - int cnt; - -#if N == 2 - if (res_ptr[N - 1] != 0) - { - count_leading_zeros (cnt, res_ptr[N - 1]); - cnt -= NUM_LEADING_ZEROS; - res_ptr[N - 1] = res_ptr[N - 1] << cnt - | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); - res_ptr[0] <<= cnt; - *expt = DBL_MIN_EXP - 1 - cnt; - } - else - { - count_leading_zeros (cnt, res_ptr[0]); - if (cnt >= NUM_LEADING_ZEROS) - { - res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS); - res_ptr[0] = 0; - } - else - { - res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt); - res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt); - } - *expt = DBL_MIN_EXP - 1 - - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt; - } + count_leading_zeros (cnt, res_ptr[0]); + exp -= BITS_PER_MP_LIMB; + } #else - int j, k, l; - - for (j = N - 1; j > 0; j--) - if (res_ptr[j] != 0) - break; - - count_leading_zeros (cnt, res_ptr[j]); - cnt -= NUM_LEADING_ZEROS; - l = N - 1 - j; - if (cnt < 0) - { - cnt += BITS_PER_MP_LIMB; - l--; - } - if (!cnt) - for (k = N - 1; k >= l; k--) - res_ptr[k] = res_ptr[k-l]; - else - { - for (k = N - 1; k > l; k--) - res_ptr[k] = res_ptr[k-l] << cnt - | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt); - res_ptr[k--] = res_ptr[0] << cnt; - } - - for (; k >= 0; k--) - res_ptr[k] = 0; - *expt = DBL_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt; + count_leading_zeros (cnt, hi); #endif - } + *zero_bits = cnt; + *expt = exp + NUM_LEADING_ZEROS - cnt; + return n; } - else - /* Add the implicit leading one bit for a normalized number. */ - res_ptr[N - 1] |= (mp_limb_t) 1 << (LDBL_MANT_DIG - 1 - - ((N - 1) * BITS_PER_MP_LIMB)); - - return N; } diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c index 3f85e0ae904..e83171b0287 100644 --- a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c @@ -30,7 +30,7 @@ mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, - int *expt, int *is_neg, + int *expt, int *zero_bits, int *is_neg, long double value) { union ieee854_long_double u; @@ -52,41 +52,38 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" #endif +/* The format does not fill the last limb. There are some zeros. */ +#define NUM_LEADING_ZEROS (N * BITS_PER_MP_LIMB - LDBL_MANT_DIG) + *zero_bits = NUM_LEADING_ZEROS; + if (u.ieee.exponent == 0) { /* A biased exponent of zero is a special case. Either it is a zero or it is a denormal number. */ if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2. */ - /* It's zero. */ - *expt = 0; + { + /* It's zero. */ + *expt = 0; + return 1; + } else { /* It is a denormal number, meaning it has no implicit leading one bit, and its exponent is in fact the format minimum. */ int cnt; + int exp = 1 - IEEE854_LONG_DOUBLE_BIAS; + int n = N; if (res_ptr[N - 1] != 0) - { - count_leading_zeros (cnt, res_ptr[N - 1]); - if (cnt != 0) - { -#if N == 2 - res_ptr[N - 1] = res_ptr[N - 1] << cnt - | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); - res_ptr[0] <<= cnt; -#else - res_ptr[N - 1] <<= cnt; -#endif - } - *expt = LDBL_MIN_EXP - 1 - cnt; - } + count_leading_zeros (cnt, res_ptr[N - 1]); else { count_leading_zeros (cnt, res_ptr[0]); - res_ptr[N - 1] = res_ptr[0] << cnt; - res_ptr[0] = 0; - *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt; + exp -= BITS_PER_MP_LIMB; + n = 1; } + *zero_bits = cnt; + *expt = exp + NUM_LEADING_ZEROS - cnt; } }