aboutsummaryrefslogtreecommitdiff
path: root/stdlib/strtod.c
diff options
context:
space:
mode:
authorRoland McGrath <roland@gnu.org>1995-02-18 01:27:10 +0000
committerRoland McGrath <roland@gnu.org>1995-02-18 01:27:10 +0000
commit28f540f45bbacd939bfd07f213bcad2bf730b1bf (patch)
tree15f07c4c43d635959c6afee96bde71fb1b3614ee /stdlib/strtod.c
downloadglibc-28f540f45bbacd939bfd07f213bcad2bf730b1bf.tar
glibc-28f540f45bbacd939bfd07f213bcad2bf730b1bf.tar.gz
glibc-28f540f45bbacd939bfd07f213bcad2bf730b1bf.tar.bz2
glibc-28f540f45bbacd939bfd07f213bcad2bf730b1bf.zip
initial import
Diffstat (limited to 'stdlib/strtod.c')
-rw-r--r--stdlib/strtod.c1027
1 files changed, 1027 insertions, 0 deletions
diff --git a/stdlib/strtod.c b/stdlib/strtod.c
new file mode 100644
index 0000000000..d647753e79
--- /dev/null
+++ b/stdlib/strtod.c
@@ -0,0 +1,1027 @@
+/* Read decimal floating point numbers.
+Copyright (C) 1995 Free Software Foundation, Inc.
+Contributed by Ulrich Drepper.
+
+This file is part of the GNU C Library.
+
+The GNU C Library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Library General Public License as
+published by the Free Software Foundation; either version 2 of the
+License, or (at your option) any later version.
+
+The GNU C Library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with the GNU C Library; see the file COPYING.LIB. If
+not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+Cambridge, MA 02139, USA. */
+
+/* Configuration part. These macros are defined by `strtold.c' and `strtof.c'
+ to produce the `long double' and `float' versions of the reader. */
+#ifndef FLOAT
+#define FLOAT double
+#define FLT DBL
+#define STRTOF strtod
+#define MPN2FLOAT __mpn_construct_double
+#define FLOAT_HUGE_VAL HUGE_VAL
+#endif
+/* End of configuration part. */
+
+#include <ctype.h>
+#include <errno.h>
+#include <float.h>
+#include <localeinfo.h>
+#include <math.h>
+#include <stdlib.h>
+#include "../stdio/gmp.h"
+#include "../stdio/gmp-impl.h"
+#include <gmp-mparam.h>
+#include "../stdio/longlong.h"
+#include "../stdio/fpioconst.h"
+
+/* #define NDEBUG 1 */
+#include <assert.h>
+
+
+/* Constants we need from float.h; select the set for the FLOAT precision. */
+#define MANT_DIG FLT##_MANT_DIG
+#define MAX_EXP FLT##_MAX_EXP
+#define MIN_EXP FLT##_MIN_EXP
+#define MAX_10_EXP FLT##_MAX_10_EXP
+#define MIN_10_EXP FLT##_MIN_10_EXP
+#define MAX_10_EXP_LOG FLT##_MAX_10_EXP_LOG
+
+
+/* Function to construct a floating point number from an MP integer
+ containing the fraction bits, a base 2 exponent, and a sign flag. */
+extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
+
+/* Definitions according to limb size used. */
+#if BITS_PER_MP_LIMB == 32
+# define MAX_DIG_PER_LIMB 9
+# define MAX_FAC_PER_LIMB 1000000000L
+#elif BITS_PER_MP_LIMB == 64
+# define MAX_DIG_PER_LIMB 19
+# define MAX_FAC_PER_LIMB 10000000000000000000L
+#else
+# error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+
+/* Local data structure. */
+static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB] =
+{ 0, 10, 100,
+ 1000, 10000, 100000,
+ 1000000, 10000000, 100000000
+#if BITS_PER_MP_LIMB > 32
+ , 1000000000, 10000000000, 100000000000,
+ 1000000000000, 10000000000000, 100000000000000,
+ 1000000000000000, 10000000000000000, 100000000000000000,
+ 1000000000000000000
+#endif
+#if BITS_PER_MP_LIMB > 64
+ #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
+#endif
+};
+
+#ifndef howmany
+#define howmany(x,y) (((x)+((y)-1))/(y))
+#endif
+#define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
+
+#define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
+#define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
+
+#define RETURN(val,end) \
+ 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)
+/* 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. */
+#define MPN_ASSIGN(dst, src) \
+ memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb))
+
+
+/* Return a floating point number of the needed type according to the given
+ multi-precision number after possible rounding. */
+static inline FLOAT
+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)
+ {
+ mp_size_t shift = MIN_EXP - 1 - exponent;
+
+ if (shift >= MANT_DIG)
+ {
+ errno = EDOM;
+ return 0.0;
+ }
+
+ more_bits |= (round_limb & ((1 << round_bit) - 1)) != 0;
+ if (shift >= BITS_PER_MP_LIMB)
+ {
+ 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;
+ }
+ else
+ {
+ round_limb = retval[0];
+ round_bit = shift - 1;
+ }
+ (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))
+ {
+ 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)
+ {
+ ++exponent;
+ (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
+ retval[RETURN_LIMB_SIZE - 1] |= 1 << (MANT_DIG % BITS_PER_MP_LIMB);
+ }
+ }
+
+ if (exponent > MAX_EXP)
+ return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
+
+ return MPN2FLOAT (retval, exponent, negative);
+}
+
+
+/* Read a multi-precision integer starting at STR with exactly DIGCNT digits
+ into N. Return the size of the number limbs in NSIZE at the first
+ character od the string that is not part of the integer as the function
+ value. If the EXPONENT is small enough to be taken as an additional
+ factor for the resulting number (see code) multiply by it. */
+static inline const char *
+str_to_mpn (const char *str, int digcnt, mp_limb *n, mp_size_t *nsize,
+ int *exponent)
+{
+ /* Number of digits for actual limb. */
+ int cnt = 0;
+ mp_limb low = 0;
+ mp_limb base;
+
+ *nsize = 0;
+ assert (digcnt > 0);
+ do
+ {
+ if (cnt == MAX_DIG_PER_LIMB)
+ {
+ if (*nsize == 0)
+ n[0] = low;
+ else
+ {
+ mp_limb cy;
+ cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
+ cy += __mpn_add_1 (n, n, *nsize, low);
+ if (cy != 0)
+ n[*nsize] = cy;
+ }
+ ++(*nsize);
+ cnt = 0;
+ low = 0;
+ }
+
+ /* There might be thousands separators or radix characters in the string.
+ But these all can be ignored because we know the format of the number
+ is correct and we have an exact number of characters to read. */
+ while (!isdigit (*str))
+ ++str;
+ low = low * 10 + *str++ - '0';
+ ++cnt;
+ }
+ while (--digcnt > 0);
+
+ if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
+ {
+ low *= _tens_in_limb[*exponent];
+ base = _tens_in_limb[cnt + *exponent];
+ *exponent = 0;
+ }
+ else
+ base = _tens_in_limb[cnt];
+
+ if (*nsize == 0)
+ {
+ n[0] = low;
+ *nsize = 1;
+ }
+ else
+ {
+ mp_limb cy;
+ cy = __mpn_mul_1 (n, n, *nsize, base);
+ cy += __mpn_add_1 (n, n, *nsize, low);
+ if (cy != 0)
+ n[(*nsize)++] = cy;
+ }
+ return str;
+}
+
+
+/* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
+ with the COUNT most significant bits of LIMB.
+
+ Tege doesn't like this function so I have to write it here myself. :)
+ --drepper */
+static inline void
+__mpn_lshift_1 (mp_limb *ptr, mp_size_t size, unsigned int count, mp_limb limb)
+{
+ if (count == BITS_PER_MP_LIMB)
+ {
+ /* Optimize the case of shifting by exactly a word:
+ just copy words, with no actual bit-shifting. */
+ mp_size_t i;
+ for (i = size - 1; i > 0; --i)
+ ptr[i] = ptr[i - 1];
+ ptr[0] = limb;
+ }
+ else
+ {
+ (void) __mpn_lshift (ptr, ptr, size, count);
+ ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
+ }
+}
+
+
+/* Return a floating point number with the value of the given string NPTR.
+ Set *ENDPTR to the character after the last used one. If the number is
+ smaller than the smallest representable number, set `errno' to ERANGE and
+ return 0.0. If the number is too big to be represented, set `errno' to
+ ERANGE and return HUGE_VAL with the approriate sign. */
+FLOAT
+STRTOF (nptr, endptr)
+ const char *nptr;
+ char **endptr;
+{
+ int negative; /* The sign of the number. */
+ MPN_VAR (num); /* MP representation of the number. */
+ int exponent; /* Exponent of the number. */
+
+ /* When we have to compute fractional digits we form a fraction with a
+ second multi-precision number (and we sometimes need a second for
+ temporary results). */
+ MPN_VAR (den);
+
+ /* Representation for the return value. */
+ mp_limb retval[RETURN_LIMB_SIZE];
+ /* Number of bits currently in result value. */
+ int bits;
+
+ /* Running pointer after the last character processed in the string. */
+ const char *cp;
+ /* Start of significant part of the number. */
+ const char *startp;
+ /* 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;
+ /* Contains the last character read. */
+ char c;
+
+ /* The radix character of the current locale. */
+ wchar_t decimal;
+#ifdef USE_GROUPING
+ /* The thousands character of the current locale. */
+ wchar_t thousands;
+ /* The numeric grouping specification of the current locale,
+ in the format described in <locale.h>. */
+ const char *grouping;
+
+ /* Check the grouping of the integer part at [BEGIN,END).
+ Return zero iff a separator is found out of place. */
+ int grouping_ok (const char *begin, const char *end)
+ {
+ if (grouping)
+ while (end > begin)
+ {
+ const char *p = end;
+ do
+ --p;
+ while (*p != thousands && p > begin);
+ if (end - 1 - p != *grouping++)
+ return 0; /* Wrong number of digits in this group. */
+ end = p; /* Correct group; trim it off the end. */
+
+ if (*grouping == 0)
+ --grouping; /* Same grouping repeats in next iteration. */
+ else if (*grouping == CHAR_MAX || *grouping < 0)
+ {
+ /* No further grouping allowed. */
+ while (end > begin)
+ if (*--end == thousands)
+ return 0;
+ }
+ }
+ return 1;
+ }
+ /* Return with no conversion if the grouping of [STARTP,CP) is bad. */
+#define CHECK_GROUPING if (! grouping_ok (startp, cp)) RETURN (0.0, nptr); else
+
+ grouping = _numeric_info->grouping; /* Cache the grouping info array. */
+ if (*grouping <= 0 || *grouping == CHAR_MAX)
+ grouping = NULL;
+ else
+ {
+ /* Figure out the thousands seperator character. */
+ if (mbtowc (&thousands_sep, _numeric_info->thousands_sep,
+ strlen (_numeric_info->thousands_sep)) <= 0)
+ thousands = (wchar_t) *_numeric_info->thousands_sep;
+ if (thousands == L'\0')
+ grouping = NULL;
+ }
+#else
+#define grouping NULL
+#define thousands L'\0'
+#define CHECK_GROUPING ((void) 0)
+#endif
+
+ /* Find the locale's decimal point character. */
+ if (mbtowc (&decimal, _numeric_info->decimal_point,
+ strlen (_numeric_info->decimal_point)) <= 0)
+ decimal = (wchar_t) *_numeric_info->decimal_point;
+
+
+ /* Prepare number representation. */
+ exponent = 0;
+ negative = 0;
+ bits = 0;
+
+ /* Parse string to get maximal legal prefix. We need the number of
+ characters of the interger part, the fractional part and the exponent. */
+ cp = nptr - 1;
+ /* Ignore leading white space. */
+ do
+ c = *++cp;
+ while (isspace (c));
+
+ /* Get sign of the result. */
+ if (c == '-')
+ {
+ negative = 1;
+ c = *++cp;
+ }
+ else if (c == '+')
+ c = *++cp;
+
+ /* Return 0.0 if no legal string is found.
+ No character is used even if a sign was found. */
+ if (!isdigit (c))
+ RETURN (0.0, nptr);
+
+ /* Record the start of the digits, in case we will check their grouping. */
+ startp = cp;
+
+ /* Ignore leading zeroes. This helps us to avoid useless computations. */
+ while (c == '0' || (thousands != L'\0' && c == thousands))
+ c = *++cp;
+
+ CHECK_GROUPING;
+
+ /* If no other digit but a '0' is found the result is 0.0.
+ Return current read pointer. */
+ if (!isdigit (c) && c != decimal)
+ RETURN (0.0, cp);
+
+ /* Remember first significant digit and read following characters until the
+ decimal point, exponent character or any non-FP number character. */
+ startp = cp;
+ dig_no = 0;
+ while (dig_no < NDIG ||
+ /* If parsing grouping info, keep going past useful digits
+ so we can check all the grouping separators. */
+ grouping)
+ {
+ if (isdigit (c))
+ ++dig_no;
+ else if (thousands == L'\0' || c != thousands)
+ /* Not a digit or separator: end of the integer part. */
+ break;
+ c = *++cp;
+ }
+
+ CHECK_GROUPING;
+
+ if (dig_no >= NDIG)
+ /* Too many digits to be representable. Assigning this to EXPONENT
+ allows us to read the full number but return HUGE_VAL after parsing. */
+ exponent = MAX_10_EXP;
+
+ /* 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;
+
+ /* Read the fractional digits. */
+ if (c == decimal)
+ {
+ if (isdigit (cp[1]))
+ {
+ ++cp;
+ do
+ {
+ ++dig_no;
+ c = *++cp;
+ }
+ while (isdigit (c));
+ }
+ }
+
+ /* Remember start of exponent (if any). */
+ expp = cp;
+
+ /* Read exponent. */
+ if (tolower (c) == 'e')
+ {
+ int exp_negative = 0;
+
+ c = *++cp;
+ if (c == '-')
+ {
+ exp_negative = 1;
+ c = *++cp;
+ }
+ else if (c == '+')
+ c = *++cp;
+
+ if (isdigit (c))
+ {
+ do
+ {
+ if ((!exp_negative && exponent * 10 + int_no > MAX_10_EXP)
+ || (exp_negative
+ && exponent * 10 + int_no > -MIN_10_EXP + MANT_DIG))
+ /* The exponent is too large/small to represent a valid
+ number. */
+ {
+ FLOAT retval;
+
+ /* Overflow or underflow. */
+ errno = ERANGE;
+ retval = (exp_negative ? 0.0 :
+ negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
+
+ /* Accept all following digits as part of the exponent. */
+ do
+ ++cp;
+ while (isdigit (*cp));
+
+ RETURN (retval, cp);
+ /* NOTREACHED */
+ }
+
+ exponent *= 10;
+ exponent += c - '0';
+ c = *++cp;
+ }
+ while (isdigit (c));
+ }
+ else
+ cp = expp;
+
+ if (exp_negative)
+ exponent = -exponent;
+ }
+
+ /* We don't want to have to work with trailing zeroes after the radix. */
+ if (dig_no > int_no)
+ {
+ while (expp[-1] == '0')
+ {
+ --expp;
+ --dig_no;
+ }
+ assert (dig_no >= int_no);
+ }
+
+ /* The whole string is parsed. Store the address of the next character. */
+ if (endptr)
+ *endptr = (char *) cp;
+
+ if (dig_no == 0)
+ return 0.0;
+
+ /* Now we have the number of digits in total and the integer digits as well
+ as the exponent and its sign. We can decide whether the read digits are
+ really integer digits or belong to the fractional part; i.e. we normalize
+ 123e-2 to 1.23. */
+ {
+ register int incr = exponent < 0 ? MAX (-int_no, exponent)
+ : MIN (dig_no - int_no, exponent);
+ int_no += incr;
+ exponent -= incr;
+ }
+
+ if (int_no + exponent > MAX_10_EXP)
+ {
+ errno = ERANGE;
+ return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
+ }
+
+ if (int_no - dig_no + exponent < MIN_10_EXP - MANT_DIG)
+ {
+ errno = ERANGE;
+ return 0.0;
+ }
+
+ if (int_no > 0)
+ {
+ /* Read the integer part as a multi-precision number to NUM. */
+ startp = str_to_mpn (startp, int_no, num, &numsize, &exponent);
+
+ if (exponent > 0)
+ {
+ /* We now multiply the gained number by the given power of ten. */
+ mp_limb *psrc = num;
+ mp_limb *pdest = den;
+ int expbit = 1;
+ const struct mp_power *ttab = &_fpioconst_pow10[0];
+
+ assert (exponent < (1 << (MAX_10_EXP_LOG + 1)));
+ do
+ {
+ if ((exponent & expbit) != 0)
+ {
+ mp_limb cy;
+ exponent ^= expbit;
+
+ /* FIXME: not the whole multiplication has to be done.
+ If we have the needed number of bits we only need the
+ information whether more non-zero bits follow. */
+ if (numsize >= ttab->arraysize - 2)
+ cy = __mpn_mul (pdest, psrc, numsize,
+ &ttab->array[2], ttab->arraysize - 2);
+ else
+ cy = __mpn_mul (pdest, &ttab->array[2],
+ ttab->arraysize - 2,
+ psrc, numsize);
+ numsize += ttab->arraysize - 2;
+ if (cy == 0)
+ --numsize;
+ SWAP (psrc, pdest);
+ }
+ expbit <<= 1;
+ ++ttab;
+ }
+ while (exponent != 0);
+
+ if (psrc == den)
+ memcpy (num, den, numsize * sizeof (mp_limb));
+ }
+
+ /* Determine how many bits of the result we already have. */
+ count_leading_zeros (bits, num[numsize - 1]);
+ bits = numsize * BITS_PER_MP_LIMB - bits;
+
+ /* We have already the first BITS bits of the result. Together with
+ the information whether more non-zero bits follow this is enough
+ to determine the result. */
+ if (bits > MANT_DIG)
+ {
+ 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;
+
+ 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);
+
+ /* Check whether any limb beside the ones in RETVAL are non-zero. */
+ for (i = 0; num[i] == 0; ++i)
+ ;
+
+ return round_and_return (retval, bits - 1, negative,
+ num[round_idx], round_bit,
+ int_no < dig_no || i < round_idx);
+ /* NOTREACHED */
+ }
+ else if (dig_no == int_no)
+ {
+ const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
+ const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
+
+ if (target_bit == is_bit)
+ {
+ memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
+ numsize * sizeof (mp_limb));
+ /* FIXME: the following loop can be avoided if we assume a
+ maximal MANT_DIG value. */
+ MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
+ }
+ else if (target_bit > is_bit)
+ {
+ (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
+ num, numsize, target_bit - is_bit);
+ /* FIXME: the following loop can be avoided if we assume a
+ maximal MANT_DIG value. */
+ MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
+ }
+ else
+ {
+ mp_limb cy;
+ assert (numsize < RETURN_LIMB_SIZE);
+
+ cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
+ num, numsize, is_bit - target_bit);
+ retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
+ /* FIXME: the following loop can be avoided if we assume a
+ maximal MANT_DIG value. */
+ MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
+ }
+
+ return round_and_return (retval, bits - 1, negative, 0, 0, 0);
+ /* NOTREACHED */
+ }
+
+ /* Store the bits we already have. */
+ memcpy (retval, num, numsize * sizeof (mp_limb));
+#if RETURN_LIMB_SIZE > 1
+ if (numsize < RETURN_LIMB_SIZE)
+ retval[numsize] = 0;
+#endif
+ }
+
+ /* We have to compute at least some of the fractional digits. */
+ {
+ /* We construct a fraction and the result of the division gives us
+ the needed digits. The denominator is 1.0 multiplied by the
+ exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
+ 123e6 gives 123 / 1000000. */
+
+ int expbit;
+ int cnt;
+ 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);
+
+ /* Construct the denominator. */
+ densize = 0;
+ expbit = 1;
+ do
+ {
+ if ((neg_exp & expbit) != 0)
+ {
+ mp_limb cy;
+ neg_exp ^= expbit;
+
+ if (densize == 0)
+ memcpy (psrc, &ttab->array[2],
+ (densize = ttab->arraysize - 2) * sizeof (mp_limb));
+ else
+ {
+ cy = __mpn_mul (pdest, &ttab->array[2], ttab->arraysize - 2,
+ psrc, densize);
+ densize += ttab->arraysize - 2;
+ if (cy == 0)
+ --densize;
+ SWAP (psrc, pdest);
+ }
+ }
+ expbit <<= 1;
+ ++ttab;
+ }
+ while (neg_exp != 0);
+
+ if (psrc == num)
+ memcpy (den, num, densize * sizeof (mp_limb));
+
+ /* Read the fractional digits from the string. */
+ (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent);
+
+
+ /* We now have to shift both numbers so that the highest bit in the
+ denominator is set. In the same process we copy the numerator to
+ a high place in the array so that the division constructs the wanted
+ digits. This is done by a "quasi fix point" number representation.
+
+ num: ddddddddddd . 0000000000000000000000
+ |--- m ---|
+ den: ddddddddddd n >= m
+ |--- n ---|
+ */
+
+ count_leading_zeros (cnt, den[densize - 1]);
+
+ (void) __mpn_lshift (den, den, densize, cnt);
+ cy = __mpn_lshift (num, num, numsize, cnt);
+ if (cy != 0)
+ num[numsize++] = cy;
+
+ /* Now we are ready for the division. But it is not necessary to
+ do a full multi-precision division because we only need a small
+ number of bits for the result. So we do not use __mpn_divmod
+ here but instead do the division here by hand and stop whenever
+ the needed number of bits is reached. The code itself comes
+ from the GNU MP Library by Torbj\"orn Granlund. */
+
+ exponent = bits;
+
+ switch (densize)
+ {
+ case 1:
+ {
+ mp_limb d, n, quot;
+ int used = 0;
+
+ n = num[0];
+ d = den[0];
+ assert (numsize == 1 && n < d);
+
+ do
+ {
+ udiv_qrnnd (quot, n, n, 0, d);
+
+#define got_limb \
+ if (bits == 0) \
+ { \
+ register int cnt; \
+ if (quot == 0) \
+ cnt = BITS_PER_MP_LIMB; \
+ else \
+ count_leading_zeros (cnt, quot); \
+ exponent -= cnt; \
+ if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
+ { \
+ used = cnt + MANT_DIG; \
+ retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
+ bits -= BITS_PER_MP_LIMB - used; \
+ } \
+ else \
+ { \
+ /* Note that we only clear the second element. */ \
+ retval[1] = 0; \
+ retval[0] = quot; \
+ bits -= cnt; \
+ } \
+ } \
+ else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
+ __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); \
+ } \
+ bits += BITS_PER_MP_LIMB
+
+ got_limb;
+ }
+ while (bits <= MANT_DIG);
+
+ return round_and_return (retval, exponent - 1, negative,
+ quot, BITS_PER_MP_LIMB - 1 - used,
+ n != 0);
+ }
+ case 2:
+ {
+ mp_limb d0, d1, n0, n1;
+ mp_limb quot = 0;
+ int used = 0;
+
+ d0 = den[0];
+ d1 = den[1];
+
+ if (numsize < densize)
+ {
+ if (bits <= 0)
+ exponent -= BITS_PER_MP_LIMB;
+ else
+ {
+ if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
+ __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
+ BITS_PER_MP_LIMB, 0);
+ else
+ {
+ used = MANT_DIG - bits;
+ if (used > 0)
+ __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
+ }
+ bits += BITS_PER_MP_LIMB;
+ }
+ n1 = num[0];
+ n0 = 0;
+ }
+ else
+ {
+ n1 = num[1];
+ n0 = num[0];
+ }
+
+ while (bits <= MANT_DIG)
+ {
+ mp_limb r;
+
+ if (n1 == d1)
+ {
+ /* QUOT should be either 111..111 or 111..110. We need
+ special treatment of this rare case as normal division
+ would give overflow. */
+ quot = ~(mp_limb) 0;
+
+ r = n0 + d1;
+ if (r < d1) /* Carry in the addition? */
+ {
+ add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
+ goto have_quot;
+ }
+ n1 = d0 - (d0 != 0);
+ n0 = -d0;
+ }
+ else
+ {
+ udiv_qrnnd (quot, r, n1, n0, d1);
+ umul_ppmm (n1, n0, d0, quot);
+ }
+
+ q_test:
+ if (n1 > r || (n1 == r && n0 > 0))
+ {
+ /* The estimated QUOT was too large. */
+ --quot;
+
+ sub_ddmmss (n1, n0, n1, n0, 0, d0);
+ r += d1;
+ if (r >= d1) /* If not carry, test QUOT again. */
+ goto q_test;
+ }
+ sub_ddmmss (n1, n0, r, 0, n1, n0);
+
+ have_quot:
+ got_limb;
+ }
+
+ return round_and_return (retval, exponent - 1, negative,
+ quot, BITS_PER_MP_LIMB - 1 - used,
+ n1 != 0 || n0 != 0);
+ }
+ default:
+ {
+ int i;
+ mp_limb cy, dX, d1, n0, n1;
+ mp_limb quot = 0;
+ int used = 0;
+
+ dX = den[densize - 1];
+ d1 = den[densize - 2];
+
+ /* The division does not work if the upper limb of the two-limb
+ numerator is greater than the denominator. */
+ if (num[numsize - 1] > dX)
+ num[numsize++] = 0;
+
+ if (numsize < densize)
+ {
+ mp_size_t empty = densize - numsize;
+
+ if (bits <= 0)
+ {
+ register int i;
+ for (i = numsize; i > 0; --i)
+ num[i + empty] = num[i - 1];
+ MPN_ZERO (num, empty + 1);
+ exponent -= empty * BITS_PER_MP_LIMB;
+ }
+ else
+ {
+ if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
+ {
+ /* We make a difference here because the compiler
+ cannot optimize the `else' case that good and
+ this reflects all currently used FLOAT types
+ and GMP implementations. */
+ register int i;
+#if RETURN_LIMB_SIZE <= 2
+ assert (empty == 1);
+ __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
+ BITS_PER_MP_LIMB, 0);
+#else
+ for (i = RETURN_LIMB_SIZE; i > empty; --i)
+ retval[i] = retval[i - empty];
+#endif
+ retval[1] = 0;
+ for (i = numsize; i > 0; --i)
+ num[i + empty] = num[i - 1];
+ MPN_ZERO (num, empty + 1);
+ }
+ else
+ {
+ used = MANT_DIG - bits;
+ if (used >= BITS_PER_MP_LIMB)
+ {
+ register int i;
+ (void) __mpn_lshift (&retval[used
+ / BITS_PER_MP_LIMB],
+ retval, RETURN_LIMB_SIZE,
+ used % BITS_PER_MP_LIMB);
+ for (i = used / BITS_PER_MP_LIMB; i >= 0; --i)
+ retval[i] = 0;
+ }
+ else if (used > 0)
+ __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
+ }
+ bits += empty * BITS_PER_MP_LIMB;
+ }
+ }
+ else
+ {
+ int i;
+ assert (numsize == densize);
+ for (i = numsize; i > 0; --i)
+ num[i] = num[i - 1];
+ }
+
+ den[densize] = 0;
+ n0 = num[densize];
+
+ while (bits <= MANT_DIG)
+ {
+ if (n0 == dX)
+ /* This might over-estimate QUOT, but it's probably not
+ worth the extra code here to find out. */
+ quot = ~(mp_limb) 0;
+ else
+ {
+ mp_limb r;
+
+ udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
+ umul_ppmm (n1, n0, d1, quot);
+
+ while (n1 > r || (n1 == r && n0 > num[densize - 2]))
+ {
+ --quot;
+ r += dX;
+ if (r < dX) /* I.e. "carry in previous addition?" */
+ break;
+ n1 -= n0 < d1;
+ n0 -= d1;
+ }
+ }
+
+ /* Possible optimization: We already have (q * n0) and (1 * n1)
+ after the calculation of QUOT. Taking advantage of this, we
+ could make this loop make two iterations less. */
+
+ cy = __mpn_submul_1 (num, den, densize + 1, quot);
+
+ if (num[densize] != cy)
+ {
+ cy = __mpn_add_n (num, num, den, densize);
+ assert (cy != 0);
+ --quot;
+ }
+ n0 = num[densize] = num[densize - 1];
+ for (i = densize - 1; i > 0; --i)
+ num[i] = num[i - 1];
+
+ got_limb;
+ }
+
+ for (i = densize - 1; num[i] != 0 && i >= 0; --i)
+ ;
+ return round_and_return (retval, exponent - 1, negative,
+ quot, BITS_PER_MP_LIMB - 1 - used,
+ i >= 0);
+ }
+ }
+ }
+
+ /* NOTREACHED */
+}