#define STRTOF strtod
#define MPN2FLOAT __mpn_construct_double
#define FLOAT_HUGE_VAL HUGE_VAL
+#define IMPLICIT_ONE 1
#endif
/* End of configuration part. */
\f
#include <localeinfo.h>
#include <math.h>
#include <stdlib.h>
-#include "../stdio/gmp.h"
-#include "../stdio/gmp-impl.h"
+#include "gmp.h"
+#include "gmp-impl.h"
#include <gmp-mparam.h>
-#include "../stdio/longlong.h"
-#include "../stdio/fpioconst.h"
+#include "longlong.h"
+#include "fpioconst.h"
-/* #define NDEBUG 1 */
+#define NDEBUG 1
#include <assert.h>
/* Constants we need from float.h; select the set for the FLOAT precision. */
#define MANT_DIG PASTE(FLT,_MANT_DIG)
+#define DIG PASTE(FLT,_DIG)
#define MAX_EXP PASTE(FLT,_MAX_EXP)
#define MIN_EXP PASTE(FLT,_MIN_EXP)
#define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
/* Local data structure. */
-static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB] =
-{ 0, 10, 100,
- 1000, 10000, 100000,
- 1000000, 10000000, 100000000
+static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB + 1] =
+{ 0, 10, 100,
+ 1000, 10000, 100000,
+ 1000000, 10000000, 100000000,
+ 1000000000
#if BITS_PER_MP_LIMB > 32
- , 1000000000, 10000000000, 100000000000,
- 1000000000000, 10000000000000, 100000000000000,
- 1000000000000000, 10000000000000000, 100000000000000000,
- 1000000000000000000
+ , 10000000000, 100000000000,
+ 1000000000000, 10000000000000, 100000000000000,
+ 1000000000000000, 10000000000000000, 100000000000000000,
+ 1000000000000000000, 10000000000000000000
#endif
#if BITS_PER_MP_LIMB > 64
#error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
do { if (endptr != 0) *endptr = (char *) end; return val; } while (0)
/* Maximum size necessary for mpn integers to hold floating point numbers. */
-#define MPNSIZE (howmany (MAX_EXP + MANT_DIG, BITS_PER_MP_LIMB) + 1)
+#define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) + 2)
/* Declare an mpn integer variable that big. */
#define MPN_VAR(name) mp_limb name[MPNSIZE]; mp_size_t name##size
/* Copy an mpn integer value. */
round_and_return (mp_limb *retval, int exponent, int negative,
mp_limb round_limb, mp_size_t round_bit, int more_bits)
{
- if (exponent < MIN_EXP)
+ if (exponent < MIN_EXP - 2 + IMPLICIT_ONE)
{
- mp_size_t shift = MIN_EXP - 1 - exponent;
+ mp_size_t shift = MIN_EXP - 2 + IMPLICIT_ONE - exponent;
if (shift >= MANT_DIG)
{
{
round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
round_bit = (shift - 1) % BITS_PER_MP_LIMB;
-#if RETURN_LIMB_SIZE <= 2
- assert (RETURN_LIMB_SIZE == 2);
- more_bits |= retval[0] != 0;
- retval[0] = retval[1];
- retval[1] = 0;
-#else
- int disp = shift / BITS_PER_MP_LIMB;
- int i = 0;
- while (retval[i] == 0 && i < disp)
- ++i;
- more_bits |= i < disp;
- for (i = disp; i < RETURN_LIMB_SIZE; ++i)
- retval[i - disp] = retval[i];
- MPN_ZERO (&retval[RETURN_LIMB_SIZE - disp], disp);
-#endif
- shift %= BITS_PER_MP_LIMB;
+
+ (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
+ RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
+ shift % BITS_PER_MP_LIMB);
+ MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
+ shift / BITS_PER_MP_LIMB);
}
- else
+ else if (shift > 0)
{
round_limb = retval[0];
round_bit = shift - 1;
+ (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
}
- (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
exponent = MIN_EXP - 2;
}
- if ((round_limb & (1 << round_bit)) != 0 &&
- (more_bits || (retval[0] & 1) != 0 ||
- (round_limb & ((1 << round_bit) - 1)) != 0))
+ if ((round_limb & (1 << round_bit)) != 0
+ && (more_bits || (retval[0] & 1) != 0
+ || (round_limb & ((1 << round_bit) - 1)) != 0))
{
mp_limb cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
- if (cy || (retval[RETURN_LIMB_SIZE - 1]
- & (1 << (MANT_DIG % BITS_PER_MP_LIMB))) != 0)
+
+ if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
+ ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
+ (retval[RETURN_LIMB_SIZE - 1]
+ & (1 << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
{
++exponent;
(void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
- retval[RETURN_LIMB_SIZE - 1] |= 1 << (MANT_DIG % BITS_PER_MP_LIMB);
+ retval[RETURN_LIMB_SIZE - 1] |= 1 << ((MANT_DIG - 1)
+ % BITS_PER_MP_LIMB);
}
+ else if (IMPLICIT_ONE && exponent == MIN_EXP - 2
+ && (retval[RETURN_LIMB_SIZE - 1]
+ & (1 << ((MANT_DIG - 1) % BITS_PER_MP_LIMB))) != 0)
+ /* The number was denormalized but now normalized. */
+ exponent = MIN_EXP - 1;
}
if (exponent > MAX_EXP)
/* Points at the character following the integer and fractional digits. */
const char *expp;
/* Total number of digit and number of digits in integer part. */
- int dig_no, int_no;
+ int dig_no, int_no, lead_zero;
/* Contains the last character read. */
char c;
/* We have the number digits in the integer part. Whether these are all or
any is really a fractional digit will be decided later. */
int_no = dig_no;
+ lead_zero = int_no == 0 ? -1 : 0;
/* Read the fractional digits. */
if (c == decimal)
{
if (isdigit (cp[1]))
{
- ++cp;
+ c = *++cp;
do
{
+ if (c != '0' && lead_zero == -1)
+ lead_zero = dig_no - int_no;
++dig_no;
c = *++cp;
}
return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
}
- if (int_no - dig_no + exponent < MIN_10_EXP - MANT_DIG)
+ if (exponent - MAX(0, lead_zero) < MIN_10_EXP - (DIG + 1))
{
errno = ERANGE;
return 0.0;
to determine the result. */
if (bits > MANT_DIG)
{
+ int i;
const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
: least_idx;
const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
- : least_idx - 1;
- int i;
+ : least_bit - 1;
if (least_bit == 0)
memcpy (retval, &num[least_idx],
RETURN_LIMB_SIZE * sizeof (mp_limb));
else
- (void) __mpn_rshift (retval, &num[least_idx],
- numsize - least_idx + 1, least_bit);
+ {
+ for (i = least_idx; i < numsize - 1; ++i)
+ retval[i - least_idx] = (num[i] >> least_bit)
+ | (num[i + 1]
+ << (BITS_PER_MP_LIMB - least_bit));
+ if (i - least_idx < RETURN_LIMB_SIZE)
+ retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
+ }
/* Check whether any limb beside the ones in RETVAL are non-zero. */
for (i = 0; num[i] == 0; ++i)
int expbit;
int cnt;
+ int neg_exp;
+ int more_bits;
mp_limb cy;
mp_limb *psrc = den;
mp_limb *pdest = num;
- int neg_exp = dig_no - int_no - exponent;
const struct mp_power *ttab = &_fpioconst_pow10[0];
assert (dig_no > int_no && exponent <= 0);
+
+ /* For the fractional part we need not process too much digits. One
+ decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
+ ceil(BITS / 3) =: N
+ digits we should have enough bits for the result. The remaining
+ decimal digits give us the information that more bits are following.
+ This can be used while rounding. (One added as a safety margin.) */
+ if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
+ {
+ dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
+ more_bits = 1;
+ }
+ else
+ more_bits = 0;
+
+ neg_exp = dig_no - int_no - exponent;
+
/* Construct the denominator. */
densize = 0;
expbit = 1;
exponent -= cnt; \
if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
{ \
- used = cnt + MANT_DIG; \
+ used = MANT_DIG + cnt; \
retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
- bits -= BITS_PER_MP_LIMB - used; \
+ bits = MANT_DIG + 1; \
} \
else \
{ \
/* Note that we only clear the second element. */ \
- retval[1] = 0; \
+ /* The conditional is determined at compile time. */ \
+ if (RETURN_LIMB_SIZE > 1) \
+ retval[1] = 0; \
retval[0] = quot; \
- bits -= cnt; \
+ bits = -cnt; \
} \
} \
else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
- __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
+ __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
quot); \
else \
{ \
used = MANT_DIG - bits; \
if (used > 0) \
- __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
+ __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
} \
bits += BITS_PER_MP_LIMB
- got_limb;
+ got_limb;
}
while (bits <= MANT_DIG);
return round_and_return (retval, exponent - 1, negative,
quot, BITS_PER_MP_LIMB - 1 - used,
- n != 0);
+ more_bits || n != 0);
}
case 2:
{
return round_and_return (retval, exponent - 1, negative,
quot, BITS_PER_MP_LIMB - 1 - used,
- n1 != 0 || n0 != 0);
+ more_bits || n1 != 0 || n0 != 0);
}
default:
{
/* The division does not work if the upper limb of the two-limb
numerator is greater than the denominator. */
- if (num[numsize - 1] > dX)
+ if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
num[numsize++] = 0;
if (numsize < densize)
got_limb;
}
- for (i = densize - 1; num[i] != 0 && i >= 0; --i)
+ for (i = densize; num[i] == 0 && i >= 0; --i)
;
return round_and_return (retval, exponent - 1, negative,
quot, BITS_PER_MP_LIMB - 1 - used,
- i >= 0);
+ more_bits || i >= 0);
}
}
}
projects / external / glibc.git / commitdiff