initial import

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 */
+}