stdlib/strtod.c

   1 /* Read decimal floating point numbers.
   2 Copyright (C) 1995, 1996 Free Software Foundation, Inc.
   3 Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
   4
   5 This file is part of the GNU C Library.
   6
   7 The GNU C Library is free software; you can redistribute it and/or
   8 modify it under the terms of the GNU Library General Public License as
   9 published by the Free Software Foundation; either version 2 of the
  10 License, or (at your option) any later version.
  11
  12 The GNU C Library is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 Library General Public License for more details.
  16
  17 You should have received a copy of the GNU Library General Public
  18 License along with the GNU C Library; see the file COPYING.LIB.  If
  19 not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  20 Boston, MA 02111-1307, USA.  */
  21
  22 /* Configuration part.  These macros are defined by `strtold.c',
  23    `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the
  24    `long double' and `float' versions of the reader.  */
  25 #ifndef FLOAT
  26 # define FLOAT          double
  27 # define FLT            DBL
  28 # ifdef USE_WIDE_CHAR
  29 #  define STRTOF        wcstod
  30 # else
  31 #  define STRTOF        strtod
  32 # endif
  33 # define MPN2FLOAT      __mpn_construct_double
  34 # define FLOAT_HUGE_VAL HUGE_VAL
  35 #endif
  36
  37 #ifdef USE_WIDE_CHAR
  38 # include <wctype.h>
  39 # include <wchar.h>
  40 # define STRING_TYPE wchar_t
  41 # define CHAR_TYPE wint_t
  42 # define L_(Ch) L##Ch
  43 # define ISSPACE(Ch) iswspace (Ch)
  44 # define TOLOWER(Ch) towlower (Ch)
  45 #else
  46 # define STRING_TYPE char
  47 # define CHAR_TYPE char
  48 # define L_(Ch) Ch
  49 # define ISSPACE(Ch) isspace (Ch)
  50 # define TOLOWER(Ch) tolower (Ch)
  51 #endif
  52 /* End of configuration part.  */
  53 \f
  54 #include <ctype.h>
  55 #include <errno.h>
  56 #include <float.h>
  57 #include "../locale/localeinfo.h"
  58 #include <math.h>
  59 #include <stdlib.h>
  60
  61 /* The gmp headers need some configuration frobs.  */
  62 #define HAVE_ALLOCA 1
  63
  64 #include "gmp.h"
  65 #include "gmp-impl.h"
  66 #include <gmp-mparam.h>
  67 #include "longlong.h"
  68 #include "fpioconst.h"
  69
  70 #define NDEBUG 1
  71 #include <assert.h>
  72
  73
  74 /* Constants we need from float.h; select the set for the FLOAT precision.  */
  75 #define MANT_DIG        PASTE(FLT,_MANT_DIG)
  76 #define DIG             PASTE(FLT,_DIG)
  77 #define MAX_EXP         PASTE(FLT,_MAX_EXP)
  78 #define MIN_EXP         PASTE(FLT,_MIN_EXP)
  79 #define MAX_10_EXP      PASTE(FLT,_MAX_10_EXP)
  80 #define MIN_10_EXP      PASTE(FLT,_MIN_10_EXP)
  81
  82 /* Extra macros required to get FLT expanded before the pasting.  */
  83 #define PASTE(a,b)      PASTE1(a,b)
  84 #define PASTE1(a,b)     a##b
  85
  86 /* Function to construct a floating point number from an MP integer
  87    containing the fraction bits, a base 2 exponent, and a sign flag.  */
  88 extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
  89 \f
  90 /* Definitions according to limb size used.  */
  91 #if     BITS_PER_MP_LIMB == 32
  92 #  define MAX_DIG_PER_LIMB      9
  93 #  define MAX_FAC_PER_LIMB      1000000000UL
  94 #elif   BITS_PER_MP_LIMB == 64
  95 #  define MAX_DIG_PER_LIMB      19
  96 #  define MAX_FAC_PER_LIMB      10000000000000000000UL
  97 #else
  98 #  error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
  99 #endif
 100
 101
 102 /* Local data structure.  */
 103 static const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] =
 104 {    0,                   10,                   100,
 105      1000,                10000,                100000,
 106      1000000,             10000000,             100000000,
 107      1000000000
 108 #if BITS_PER_MP_LIMB > 32
 109                ,          10000000000,          100000000000,
 110      1000000000000,       10000000000000,       100000000000000,
 111      1000000000000000,    10000000000000000,    100000000000000000,
 112      1000000000000000000, 10000000000000000000U
 113 #endif
 114 #if BITS_PER_MP_LIMB > 64
 115   #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
 116 #endif
 117 };
 118 \f
 119 #ifndef howmany
 120 #define howmany(x,y)            (((x)+((y)-1))/(y))
 121 #endif
 122 #define SWAP(x, y)              ({ typeof(x) _tmp = x; x = y; y = _tmp; })
 123
 124 #define NDIG                    (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
 125 #define RETURN_LIMB_SIZE                howmany (MANT_DIG, BITS_PER_MP_LIMB)
 126
 127 #define RETURN(val,end)                                                       \
 128     do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end);                 \
 129          return val; } while (0)
 130
 131 /* Maximum size necessary for mpn integers to hold floating point numbers.  */
 132 #define MPNSIZE         (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
 133                          + 2)
 134 /* Declare an mpn integer variable that big.  */
 135 #define MPN_VAR(name)   mp_limb_t name[MPNSIZE]; mp_size_t name##size
 136 /* Copy an mpn integer value.  */
 137 #define MPN_ASSIGN(dst, src) \
 138         memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
 139
 140
 141 /* Return a floating point number of the needed type according to the given
 142    multi-precision number after possible rounding.  */
 143 static inline FLOAT
 144 round_and_return (mp_limb_t *retval, int exponent, int negative,
 145                   mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
 146 {
 147   if (exponent < MIN_EXP - 1)
 148     {
 149       mp_size_t shift = MIN_EXP - 1 - exponent;
 150
 151       if (shift > MANT_DIG)
 152         {
 153           errno = EDOM;
 154           return 0.0;
 155         }
 156
 157       more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
 158       if (shift == MANT_DIG)
 159         /* This is a special case to handle the very seldom case where
 160            the mantissa will be empty after the shift.  */
 161         {
 162           int i;
 163
 164           round_limb = retval[RETURN_LIMB_SIZE - 1];
 165           round_bit = BITS_PER_MP_LIMB - 1;
 166           for (i = 0; i < RETURN_LIMB_SIZE; ++i)
 167             more_bits |= retval[i] != 0;
 168           MPN_ZERO (retval, RETURN_LIMB_SIZE);
 169         }
 170       else if (shift >= BITS_PER_MP_LIMB)
 171         {
 172           int i;
 173
 174           round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
 175           round_bit = (shift - 1) % BITS_PER_MP_LIMB;
 176           for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
 177             more_bits |= retval[i] != 0;
 178           more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
 179                         != 0);
 180
 181           (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
 182                                RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
 183                                shift % BITS_PER_MP_LIMB);
 184           MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
 185                     shift / BITS_PER_MP_LIMB);
 186         }
 187       else if (shift > 0)
 188         {
 189           round_limb = retval[0];
 190           round_bit = shift - 1;
 191           (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
 192         }
 193       exponent = MIN_EXP - 2;
 194     }
 195
 196   if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
 197       && (more_bits || (retval[0] & 1) != 0
 198           || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
 199     {
 200       mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
 201
 202       if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
 203           ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
 204            (retval[RETURN_LIMB_SIZE - 1]
 205             & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
 206         {
 207           ++exponent;
 208           (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
 209           retval[RETURN_LIMB_SIZE - 1]
 210             |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
 211         }
 212       else if (exponent == MIN_EXP - 2
 213                && (retval[RETURN_LIMB_SIZE - 1]
 214                    & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
 215                != 0)
 216           /* The number was denormalized but now normalized.  */
 217         exponent = MIN_EXP - 1;
 218     }
 219
 220   if (exponent > MAX_EXP)
 221     return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
 222
 223   return MPN2FLOAT (retval, exponent, negative);
 224 }
 225
 226
 227 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
 228    into N.  Return the size of the number limbs in NSIZE at the first
 229    character od the string that is not part of the integer as the function
 230    value.  If the EXPONENT is small enough to be taken as an additional
 231    factor for the resulting number (see code) multiply by it.  */
 232 static inline const STRING_TYPE *
 233 str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
 234             int *exponent)
 235 {
 236   /* Number of digits for actual limb.  */
 237   int cnt = 0;
 238   mp_limb_t low = 0;
 239   mp_limb_t base;
 240
 241   *nsize = 0;
 242   assert (digcnt > 0);
 243   do
 244     {
 245       if (cnt == MAX_DIG_PER_LIMB)
 246         {
 247           if (*nsize == 0)
 248             n[0] = low;
 249           else
 250             {
 251               mp_limb_t cy;
 252               cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
 253               cy += __mpn_add_1 (n, n, *nsize, low);
 254               if (cy != 0)
 255                 n[*nsize] = cy;
 256             }
 257           ++(*nsize);
 258           cnt = 0;
 259           low = 0;
 260         }
 261
 262       /* There might be thousands separators or radix characters in the string.
 263          But these all can be ignored because we know the format of the number
 264          is correct and we have an exact number of characters to read.  */
 265       while (*str < L_('0') || *str > L_('9'))
 266         ++str;
 267       low = low * 10 + *str++ - L_('0');
 268       ++cnt;
 269     }
 270   while (--digcnt > 0);
 271
 272   if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
 273     {
 274       low *= _tens_in_limb[*exponent];
 275       base = _tens_in_limb[cnt + *exponent];
 276       *exponent = 0;
 277     }
 278   else
 279     base = _tens_in_limb[cnt];
 280
 281   if (*nsize == 0)
 282     {
 283       n[0] = low;
 284       *nsize = 1;
 285     }
 286   else
 287     {
 288       mp_limb_t cy;
 289       cy = __mpn_mul_1 (n, n, *nsize, base);
 290       cy += __mpn_add_1 (n, n, *nsize, low);
 291       if (cy != 0)
 292         n[(*nsize)++] = cy;
 293     }
 294   return str;
 295 }
 296
 297
 298 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
 299    with the COUNT most significant bits of LIMB.
 300
 301    Tege doesn't like this function so I have to write it here myself. :)
 302    --drepper */
 303 static inline void
 304 __mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count,
 305                 mp_limb_t limb)
 306 {
 307   if (count == BITS_PER_MP_LIMB)
 308     {
 309       /* Optimize the case of shifting by exactly a word:
 310          just copy words, with no actual bit-shifting.  */
 311       mp_size_t i;
 312       for (i = size - 1; i > 0; --i)
 313         ptr[i] = ptr[i - 1];
 314       ptr[0] = limb;
 315     }
 316   else
 317     {
 318       (void) __mpn_lshift (ptr, ptr, size, count);
 319       ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
 320     }
 321 }
 322
 323
 324 #define INTERNAL(x) INTERNAL1(x)
 325 #define INTERNAL1(x) __##x##_internal
 326
 327 /* This file defines a function to check for correct grouping.  */
 328 #include "grouping.h"
 329
 330
 331 /* Return a floating point number with the value of the given string NPTR.
 332    Set *ENDPTR to the character after the last used one.  If the number is
 333    smaller than the smallest representable number, set `errno' to ERANGE and
 334    return 0.0.  If the number is too big to be represented, set `errno' to
 335    ERANGE and return HUGE_VAL with the approriate sign.  */
 336 FLOAT
 337 INTERNAL (STRTOF) (nptr, endptr, group)
 338      const STRING_TYPE *nptr;
 339      STRING_TYPE **endptr;
 340      int group;
 341 {
 342   int negative;                 /* The sign of the number.  */
 343   MPN_VAR (num);                /* MP representation of the number.  */
 344   int exponent;                 /* Exponent of the number.  */
 345
 346   /* When we have to compute fractional digits we form a fraction with a
 347      second multi-precision number (and we sometimes need a second for
 348      temporary results).  */
 349   MPN_VAR (den);
 350
 351   /* Representation for the return value.  */
 352   mp_limb_t retval[RETURN_LIMB_SIZE];
 353   /* Number of bits currently in result value.  */
 354   int bits;
 355
 356   /* Running pointer after the last character processed in the string.  */
 357   const STRING_TYPE *cp, *tp;
 358   /* Start of significant part of the number.  */
 359   const STRING_TYPE *startp, *start_of_digits;
 360   /* Points at the character following the integer and fractional digits.  */
 361   const STRING_TYPE *expp;
 362   /* Total number of digit and number of digits in integer part.  */
 363   int dig_no, int_no, lead_zero;
 364   /* Contains the last character read.  */
 365   CHAR_TYPE c;
 366
 367 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
 368    there.  So define it ourselves if it remains undefined.  */
 369 #ifndef _WINT_T
 370   typedef unsigned int wint_t;
 371 #endif
 372   /* The radix character of the current locale.  */
 373   wint_t decimal;
 374   /* The thousands character of the current locale.  */
 375   wint_t thousands;
 376   /* The numeric grouping specification of the current locale,
 377      in the format described in <locale.h>.  */
 378   const char *grouping;
 379
 380   if (group)
 381     {
 382       grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
 383       if (*grouping <= 0 || *grouping == CHAR_MAX)
 384         grouping = NULL;
 385       else
 386         {
 387           /* Figure out the thousands separator character.  */
 388           if (mbtowc ((wchar_t *) &thousands,
 389                       _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP),
 390                       strlen (_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP))) <= 0)
 391             thousands = (wint_t) *_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
 392           if (thousands == L'\0')
 393             grouping = NULL;
 394         }
 395     }
 396   else
 397     {
 398       grouping = NULL;
 399       thousands = L'\0';
 400     }
 401
 402   /* Find the locale's decimal point character.  */
 403   if (mbtowc ((wchar_t *) &decimal, _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT),
 404               strlen (_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT))) <= 0)
 405     decimal = (wint_t) *_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
 406
 407
 408   /* Prepare number representation.  */
 409   exponent = 0;
 410   negative = 0;
 411   bits = 0;
 412
 413   /* Parse string to get maximal legal prefix.  We need the number of
 414      characters of the integer part, the fractional part and the exponent.  */
 415   cp = nptr - 1;
 416   /* Ignore leading white space.  */
 417   do
 418     c = *++cp;
 419   while (ISSPACE (c));
 420
 421   /* Get sign of the result.  */
 422   if (c == L_('-'))
 423     {
 424       negative = 1;
 425       c = *++cp;
 426     }
 427   else if (c == L_('+'))
 428     c = *++cp;
 429
 430   /* Return 0.0 if no legal string is found.
 431      No character is used even if a sign was found.  */
 432   if ((c < L_('0') || c > L_('9'))
 433       && (c != decimal || cp[1] < L_('0') || cp[1] > L_('9')))
 434     RETURN (0.0, nptr);
 435
 436   /* Record the start of the digits, in case we will check their grouping.  */
 437   start_of_digits = startp = cp;
 438
 439   /* Ignore leading zeroes.  This helps us to avoid useless computations.  */
 440   while (c == L_('0') || (thousands != L'\0' && c == thousands))
 441     c = *++cp;
 442
 443   /* If no other digit but a '0' is found the result is 0.0.
 444      Return current read pointer.  */
 445   if ((c < L_('0') || c > L_('9')) && c != decimal)
 446     {
 447       tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping);
 448       /* If TP is at the start of the digits, there was no correctly
 449          grouped prefix of the string; so no number found.  */
 450       RETURN (0.0, tp == start_of_digits ? nptr : tp);
 451     }
 452
 453   /* Remember first significant digit and read following characters until the
 454      decimal point, exponent character or any non-FP number character.  */
 455   startp = cp;
 456   dig_no = 0;
 457   while (dig_no < NDIG ||
 458          /* If parsing grouping info, keep going past useful digits
 459             so we can check all the grouping separators.  */
 460          grouping)
 461     {
 462       if (c >= L_('0') && c <= L_('9'))
 463         ++dig_no;
 464       else if (thousands == L'\0' || c != thousands)
 465         /* Not a digit or separator: end of the integer part.  */
 466         break;
 467       c = *++cp;
 468     }
 469
 470   if (grouping && dig_no > 0)
 471     {
 472       /* Check the grouping of the digits.  */
 473       tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping);
 474       if (cp != tp)
 475         {
 476           /* Less than the entire string was correctly grouped.  */
 477
 478           if (tp == start_of_digits)
 479             /* No valid group of numbers at all: no valid number.  */
 480             RETURN (0.0, nptr);
 481
 482           if (tp < startp)
 483             /* The number is validly grouped, but consists
 484                only of zeroes.  The whole value is zero.  */
 485             RETURN (0.0, tp);
 486
 487           /* Recompute DIG_NO so we won't read more digits than
 488              are properly grouped.  */
 489           cp = tp;
 490           dig_no = 0;
 491           for (tp = startp; tp < cp; ++tp)
 492             if (*tp >= L_('0') && *tp <= L_('9'))
 493               ++dig_no;
 494
 495           int_no = dig_no;
 496           lead_zero = 0;
 497
 498           goto number_parsed;
 499         }
 500     }
 501
 502   if (dig_no >= NDIG)
 503     /* Too many digits to be representable.  Assigning this to EXPONENT
 504        allows us to read the full number but return HUGE_VAL after parsing.  */
 505     exponent = MAX_10_EXP;
 506
 507   /* We have the number digits in the integer part.  Whether these are all or
 508      any is really a fractional digit will be decided later.  */
 509   int_no = dig_no;
 510   lead_zero = int_no == 0 ? -1 : 0;
 511
 512   /* Read the fractional digits.  A special case are the 'american style'
 513      numbers like `16.' i.e. with decimal but without trailing digits.  */
 514   if (c == decimal)
 515     {
 516       c = *++cp;
 517       while (c >= L_('0') && c <= L_('9'))
 518         {
 519           if (c != L_('0') && lead_zero == -1)
 520             lead_zero = dig_no - int_no;
 521           ++dig_no;
 522           c = *++cp;
 523         }
 524     }
 525
 526   /* Remember start of exponent (if any).  */
 527   expp = cp;
 528
 529   /* Read exponent.  */
 530   if (TOLOWER (c) == L_('e'))
 531     {
 532       int exp_negative = 0;
 533
 534       c = *++cp;
 535       if (c == L_('-'))
 536         {
 537           exp_negative = 1;
 538           c = *++cp;
 539         }
 540       else if (c == L_('+'))
 541         c = *++cp;
 542
 543       if (c >= L_('0') && c <= L_('9'))
 544         {
 545           int exp_limit;
 546
 547           /* Get the exponent limit. */
 548           exp_limit = exp_negative ?
 549                 -MIN_10_EXP + MANT_DIG - int_no :
 550                 MAX_10_EXP - int_no + lead_zero;
 551
 552           do
 553             {
 554               exponent *= 10;
 555
 556               if (exponent > exp_limit)
 557                 /* The exponent is too large/small to represent a valid
 558                    number.  */
 559                 {
 560                   FLOAT retval;
 561
 562                   /* Overflow or underflow.  */
 563                   errno = ERANGE;
 564                   retval = (exp_negative ? 0.0 :
 565                             negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
 566
 567                   /* Accept all following digits as part of the exponent.  */
 568                   do
 569                     ++cp;
 570                   while (*cp >= L_('0') && *cp <= L_('9'));
 571
 572                   RETURN (retval, cp);
 573                   /* NOTREACHED */
 574                 }
 575
 576               exponent += c - L_('0');
 577               c = *++cp;
 578             }
 579           while (c >= L_('0') && c <= L_('9'));
 580
 581           if (exp_negative)
 582             exponent = -exponent;
 583         }
 584       else
 585         cp = expp;
 586     }
 587
 588   /* We don't want to have to work with trailing zeroes after the radix.  */
 589   if (dig_no > int_no)
 590     {
 591       while (expp[-1] == L_('0'))
 592         {
 593           --expp;
 594           --dig_no;
 595         }
 596       assert (dig_no >= int_no);
 597     }
 598
 599  number_parsed:
 600
 601   /* The whole string is parsed.  Store the address of the next character.  */
 602   if (endptr)
 603     *endptr = (STRING_TYPE *) cp;
 604
 605   if (dig_no == 0)
 606     return 0.0;
 607
 608   if (lead_zero)
 609     {
 610       /* Find the decimal point */
 611       while (*startp != decimal)
 612         ++startp;
 613       startp += lead_zero + 1;
 614       exponent -= lead_zero;
 615       dig_no -= lead_zero;
 616     }
 617
 618   /* Now we have the number of digits in total and the integer digits as well
 619      as the exponent and its sign.  We can decide whether the read digits are
 620      really integer digits or belong to the fractional part; i.e. we normalize
 621      123e-2 to 1.23.  */
 622   {
 623     register int incr = exponent < 0 ? MAX (-int_no, exponent)
 624                                      : MIN (dig_no - int_no, exponent);
 625     int_no += incr;
 626     exponent -= incr;
 627   }
 628
 629   if (int_no + exponent > MAX_10_EXP + 1)
 630     {
 631       errno = ERANGE;
 632       return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
 633     }
 634
 635   if (exponent < MIN_10_EXP - (DIG + 1))
 636     {
 637       errno = ERANGE;
 638       return 0.0;
 639     }
 640
 641   if (int_no > 0)
 642     {
 643       /* Read the integer part as a multi-precision number to NUM.  */
 644       startp = str_to_mpn (startp, int_no, num, &numsize, &exponent);
 645
 646       if (exponent > 0)
 647         {
 648           /* We now multiply the gained number by the given power of ten.  */
 649           mp_limb_t *psrc = num;
 650           mp_limb_t *pdest = den;
 651           int expbit = 1;
 652           const struct mp_power *ttab = &_fpioconst_pow10[0];
 653
 654           do
 655             {
 656               if ((exponent & expbit) != 0)
 657                 {
 658                   mp_limb_t cy;
 659                   exponent ^= expbit;
 660
 661                   /* FIXME: not the whole multiplication has to be done.
 662                      If we have the needed number of bits we only need the
 663                      information whether more non-zero bits follow.  */
 664                   if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
 665                     cy = __mpn_mul (pdest, psrc, numsize,
 666                                     &ttab->array[_FPIO_CONST_OFFSET],
 667                                     ttab->arraysize - _FPIO_CONST_OFFSET);
 668                   else
 669                     cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
 670                                     ttab->arraysize - _FPIO_CONST_OFFSET,
 671                                     psrc, numsize);
 672                   numsize += ttab->arraysize - _FPIO_CONST_OFFSET;
 673                   if (cy == 0)
 674                     --numsize;
 675                   SWAP (psrc, pdest);
 676                 }
 677               expbit <<= 1;
 678               ++ttab;
 679             }
 680           while (exponent != 0);
 681
 682           if (psrc == den)
 683             memcpy (num, den, numsize * sizeof (mp_limb_t));
 684         }
 685
 686       /* Determine how many bits of the result we already have.  */
 687       count_leading_zeros (bits, num[numsize - 1]);
 688       bits = numsize * BITS_PER_MP_LIMB - bits;
 689
 690       /* Now we know the exponent of the number in base two.
 691          Check it against the maximum possible exponent.  */
 692       if (bits > MAX_EXP)
 693         {
 694           errno = ERANGE;
 695           return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
 696         }
 697
 698       /* We have already the first BITS bits of the result.  Together with
 699          the information whether more non-zero bits follow this is enough
 700          to determine the result.  */
 701       if (bits > MANT_DIG)
 702         {
 703           int i;
 704           const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
 705           const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
 706           const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
 707                                                      : least_idx;
 708           const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
 709                                                      : least_bit - 1;
 710
 711           if (least_bit == 0)
 712             memcpy (retval, &num[least_idx],
 713                     RETURN_LIMB_SIZE * sizeof (mp_limb_t));
 714           else
 715             {
 716               for (i = least_idx; i < numsize - 1; ++i)
 717                 retval[i - least_idx] = (num[i] >> least_bit)
 718                                         | (num[i + 1]
 719                                            << (BITS_PER_MP_LIMB - least_bit));
 720               if (i - least_idx < RETURN_LIMB_SIZE)
 721                 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
 722             }
 723
 724           /* Check whether any limb beside the ones in RETVAL are non-zero.  */
 725           for (i = 0; num[i] == 0; ++i)
 726             ;
 727
 728           return round_and_return (retval, bits - 1, negative,
 729                                    num[round_idx], round_bit,
 730                                    int_no < dig_no || i < round_idx);
 731           /* NOTREACHED */
 732         }
 733       else if (dig_no == int_no)
 734         {
 735           const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
 736           const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
 737
 738           if (target_bit == is_bit)
 739             {
 740               memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
 741                       numsize * sizeof (mp_limb_t));
 742               /* FIXME: the following loop can be avoided if we assume a
 743                  maximal MANT_DIG value.  */
 744               MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
 745             }
 746           else if (target_bit > is_bit)
 747             {
 748               (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
 749                                    num, numsize, target_bit - is_bit);
 750               /* FIXME: the following loop can be avoided if we assume a
 751                  maximal MANT_DIG value.  */
 752               MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
 753             }
 754           else
 755             {
 756               mp_limb_t cy;
 757               assert (numsize < RETURN_LIMB_SIZE);
 758
 759               cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
 760                                  num, numsize, is_bit - target_bit);
 761               retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
 762               /* FIXME: the following loop can be avoided if we assume a
 763                  maximal MANT_DIG value.  */
 764               MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
 765             }
 766
 767           return round_and_return (retval, bits - 1, negative, 0, 0, 0);
 768           /* NOTREACHED */
 769         }
 770
 771       /* Store the bits we already have.  */
 772       memcpy (retval, num, numsize * sizeof (mp_limb_t));
 773 #if RETURN_LIMB_SIZE > 1
 774       if (numsize < RETURN_LIMB_SIZE)
 775         retval[numsize] = 0;
 776 #endif
 777     }
 778
 779   /* We have to compute at least some of the fractional digits.  */
 780   {
 781     /* We construct a fraction and the result of the division gives us
 782        the needed digits.  The denominator is 1.0 multiplied by the
 783        exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
 784        123e-6 gives 123 / 1000000.  */
 785
 786     int expbit;
 787     int cnt;
 788     int neg_exp;
 789     int more_bits;
 790     mp_limb_t cy;
 791     mp_limb_t *psrc = den;
 792     mp_limb_t *pdest = num;
 793     const struct mp_power *ttab = &_fpioconst_pow10[0];
 794
 795     assert (dig_no > int_no && exponent <= 0);
 796
 797
 798     /* For the fractional part we need not process too much digits.  One
 799        decimal digits gives us log_2(10) ~ 3.32 bits.  If we now compute
 800                         ceil(BITS / 3) =: N
 801        digits we should have enough bits for the result.  The remaining
 802        decimal digits give us the information that more bits are following.
 803        This can be used while rounding.  (One added as a safety margin.)  */
 804     if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
 805       {
 806         dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
 807         more_bits = 1;
 808       }
 809     else
 810       more_bits = 0;
 811
 812     neg_exp = dig_no - int_no - exponent;
 813
 814     /* Construct the denominator.  */
 815     densize = 0;
 816     expbit = 1;
 817     do
 818       {
 819         if ((neg_exp & expbit) != 0)
 820           {
 821             mp_limb_t cy;
 822             neg_exp ^= expbit;
 823
 824             if (densize == 0)
 825               {
 826                 densize = ttab->arraysize - _FPIO_CONST_OFFSET;
 827                 memcpy (psrc, &ttab->array[_FPIO_CONST_OFFSET],
 828                         densize * sizeof (mp_limb_t));
 829               }
 830             else
 831               {
 832                 cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
 833                                 ttab->arraysize - _FPIO_CONST_OFFSET,
 834                                 psrc, densize);
 835                 densize += ttab->arraysize - _FPIO_CONST_OFFSET;
 836                 if (cy == 0)
 837                   --densize;
 838                 SWAP (psrc, pdest);
 839               }
 840           }
 841         expbit <<= 1;
 842         ++ttab;
 843       }
 844     while (neg_exp != 0);
 845
 846     if (psrc == num)
 847       memcpy (den, num, densize * sizeof (mp_limb_t));
 848
 849     /* Read the fractional digits from the string.  */
 850     (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent);
 851
 852
 853     /* We now have to shift both numbers so that the highest bit in the
 854        denominator is set.  In the same process we copy the numerator to
 855        a high place in the array so that the division constructs the wanted
 856        digits.  This is done by a "quasi fix point" number representation.
 857
 858        num:   ddddddddddd . 0000000000000000000000
 859               |--- m ---|
 860        den:                            ddddddddddd      n >= m
 861                                        |--- n ---|
 862      */
 863
 864     count_leading_zeros (cnt, den[densize - 1]);
 865
 866     (void) __mpn_lshift (den, den, densize, cnt);
 867     cy = __mpn_lshift (num, num, numsize, cnt);
 868     if (cy != 0)
 869       num[numsize++] = cy;
 870
 871     /* Now we are ready for the division.  But it is not necessary to
 872        do a full multi-precision division because we only need a small
 873        number of bits for the result.  So we do not use __mpn_divmod
 874        here but instead do the division here by hand and stop whenever
 875        the needed number of bits is reached.  The code itself comes
 876        from the GNU MP Library by Torbj\"orn Granlund.  */
 877
 878     exponent = bits;
 879
 880     switch (densize)
 881       {
 882       case 1:
 883         {
 884           mp_limb_t d, n, quot;
 885           int used = 0;
 886
 887           n = num[0];
 888           d = den[0];
 889           assert (numsize == 1 && n < d);
 890
 891           do
 892             {
 893               udiv_qrnnd (quot, n, n, 0, d);
 894
 895 #define got_limb                                                              \
 896               if (bits == 0)                                                  \
 897                 {                                                             \
 898                   register int cnt;                                           \
 899                   if (quot == 0)                                              \
 900                     cnt = BITS_PER_MP_LIMB;                                   \
 901                   else                                                        \
 902                     count_leading_zeros (cnt, quot);                          \
 903                   exponent -= cnt;                                            \
 904                   if (BITS_PER_MP_LIMB - cnt > MANT_DIG)                      \
 905                     {                                                         \
 906                       used = MANT_DIG + cnt;                                  \
 907                       retval[0] = quot >> (BITS_PER_MP_LIMB - used);          \
 908                       bits = MANT_DIG + 1;                                    \
 909                     }                                                         \
 910                   else                                                        \
 911                     {                                                         \
 912                       /* Note that we only clear the second element.  */      \
 913                       /* The conditional is determined at compile time.  */   \
 914                       if (RETURN_LIMB_SIZE > 1)                               \
 915                         retval[1] = 0;                                        \
 916                       retval[0] = quot;                                       \
 917                       bits = -cnt;                                            \
 918                     }                                                         \
 919                 }                                                             \
 920               else if (bits + BITS_PER_MP_LIMB <= MANT_DIG)                   \
 921                 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB,   \
 922                                 quot);                                        \
 923               else                                                            \
 924                 {                                                             \
 925                   used = MANT_DIG - bits;                                     \
 926                   if (used > 0)                                               \
 927                     __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot);    \
 928                 }                                                             \
 929               bits += BITS_PER_MP_LIMB
 930
 931               got_limb;
 932             }
 933           while (bits <= MANT_DIG);
 934
 935           return round_and_return (retval, exponent - 1, negative,
 936                                    quot, BITS_PER_MP_LIMB - 1 - used,
 937                                    more_bits || n != 0);
 938         }
 939       case 2:
 940         {
 941           mp_limb_t d0, d1, n0, n1;
 942           mp_limb_t quot = 0;
 943           int used = 0;
 944
 945           d0 = den[0];
 946           d1 = den[1];
 947
 948           if (numsize < densize)
 949             {
 950               if (num[0] >= d1)
 951                 {
 952                   /* The numerator of the number occupies fewer bits than
 953                      the denominator but the one limb is bigger than the
 954                      high limb of the numerator.  */
 955                   n1 = 0;
 956                   n0 = num[0];
 957                 }
 958               else
 959                 {
 960                   if (bits <= 0)
 961                     exponent -= BITS_PER_MP_LIMB;
 962                   else
 963                     {
 964                       if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
 965                         __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
 966                                         BITS_PER_MP_LIMB, 0);
 967                       else
 968                         {
 969                           used = MANT_DIG - bits;
 970                           if (used > 0)
 971                             __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
 972                         }
 973                       bits += BITS_PER_MP_LIMB;
 974                     }
 975                   n1 = num[0];
 976                   n0 = 0;
 977                 }
 978             }
 979           else
 980             {
 981               n1 = num[1];
 982               n0 = num[0];
 983             }
 984
 985           while (bits <= MANT_DIG)
 986             {
 987               mp_limb_t r;
 988
 989               if (n1 == d1)
 990                 {
 991                   /* QUOT should be either 111..111 or 111..110.  We need
 992                      special treatment of this rare case as normal division
 993                      would give overflow.  */
 994                   quot = ~(mp_limb_t) 0;
 995
 996                   r = n0 + d1;
 997                   if (r < d1)   /* Carry in the addition?  */
 998                     {
 999                       add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
1000                       goto have_quot;
1001                     }
1002                   n1 = d0 - (d0 != 0);
1003                   n0 = -d0;
1004                 }
1005               else
1006                 {
1007                   udiv_qrnnd (quot, r, n1, n0, d1);
1008                   umul_ppmm (n1, n0, d0, quot);
1009                 }
1010
1011             q_test:
1012               if (n1 > r || (n1 == r && n0 > 0))
1013                 {
1014                   /* The estimated QUOT was too large.  */
1015                   --quot;
1016
1017                   sub_ddmmss (n1, n0, n1, n0, 0, d0);
1018                   r += d1;
1019                   if (r >= d1)  /* If not carry, test QUOT again.  */
1020                     goto q_test;
1021                 }
1022               sub_ddmmss (n1, n0, r, 0, n1, n0);
1023
1024             have_quot:
1025               got_limb;
1026             }
1027
1028           return round_and_return (retval, exponent - 1, negative,
1029                                    quot, BITS_PER_MP_LIMB - 1 - used,
1030                                    more_bits || n1 != 0 || n0 != 0);
1031         }
1032       default:
1033         {
1034           int i;
1035           mp_limb_t cy, dX, d1, n0, n1;
1036           mp_limb_t quot = 0;
1037           int used = 0;
1038
1039           dX = den[densize - 1];
1040           d1 = den[densize - 2];
1041
1042           /* The division does not work if the upper limb of the two-limb
1043              numerator is greater than the denominator.  */
1044           if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1045             num[numsize++] = 0;
1046
1047           if (numsize < densize)
1048             {
1049               mp_size_t empty = densize - numsize;
1050
1051               if (bits <= 0)
1052                 {
1053                   register int i;
1054                   for (i = numsize; i > 0; --i)
1055                     num[i + empty] = num[i - 1];
1056                   MPN_ZERO (num, empty + 1);
1057                   exponent -= empty * BITS_PER_MP_LIMB;
1058                 }
1059               else
1060                 {
1061                   if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1062                     {
1063                       /* We make a difference here because the compiler
1064                          cannot optimize the `else' case that good and
1065                          this reflects all currently used FLOAT types
1066                          and GMP implementations.  */
1067                       register int i;
1068 #if RETURN_LIMB_SIZE <= 2
1069                       assert (empty == 1);
1070                       __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1071                                       BITS_PER_MP_LIMB, 0);
1072 #else
1073                       for (i = RETURN_LIMB_SIZE; i > empty; --i)
1074                         retval[i] = retval[i - empty];
1075 #endif
1076                       retval[1] = 0;
1077                       for (i = numsize; i > 0; --i)
1078                         num[i + empty] = num[i - 1];
1079                       MPN_ZERO (num, empty + 1);
1080                     }
1081                   else
1082                     {
1083                       used = MANT_DIG - bits;
1084                       if (used >= BITS_PER_MP_LIMB)
1085                         {
1086                           register int i;
1087                           (void) __mpn_lshift (&retval[used
1088                                                        / BITS_PER_MP_LIMB],
1089                                                retval, RETURN_LIMB_SIZE,
1090                                                used % BITS_PER_MP_LIMB);
1091                           for (i = used / BITS_PER_MP_LIMB; i >= 0; --i)
1092                             retval[i] = 0;
1093                         }
1094                       else if (used > 0)
1095                         __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1096                     }
1097                   bits += empty * BITS_PER_MP_LIMB;
1098                 }
1099             }
1100           else
1101             {
1102               int i;
1103               assert (numsize == densize);
1104               for (i = numsize; i > 0; --i)
1105                 num[i] = num[i - 1];
1106             }
1107
1108           den[densize] = 0;
1109           n0 = num[densize];
1110
1111           while (bits <= MANT_DIG)
1112             {
1113               if (n0 == dX)
1114                 /* This might over-estimate QUOT, but it's probably not
1115                    worth the extra code here to find out.  */
1116                 quot = ~(mp_limb_t) 0;
1117               else
1118                 {
1119                   mp_limb_t r;
1120
1121                   udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1122                   umul_ppmm (n1, n0, d1, quot);
1123
1124                   while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1125                     {
1126                       --quot;
1127                       r += dX;
1128                       if (r < dX) /* I.e. "carry in previous addition?" */
1129                         break;
1130                       n1 -= n0 < d1;
1131                       n0 -= d1;
1132                     }
1133                 }
1134
1135               /* Possible optimization: We already have (q * n0) and (1 * n1)
1136                  after the calculation of QUOT.  Taking advantage of this, we
1137                  could make this loop make two iterations less.  */
1138
1139               cy = __mpn_submul_1 (num, den, densize + 1, quot);
1140
1141               if (num[densize] != cy)
1142                 {
1143                   cy = __mpn_add_n (num, num, den, densize);
1144                   assert (cy != 0);
1145                   --quot;
1146                 }
1147               n0 = num[densize] = num[densize - 1];
1148               for (i = densize - 1; i > 0; --i)
1149                 num[i] = num[i - 1];
1150
1151               got_limb;
1152             }
1153
1154           for (i = densize; num[i] == 0 && i >= 0; --i)
1155             ;
1156           return round_and_return (retval, exponent - 1, negative,
1157                                    quot, BITS_PER_MP_LIMB - 1 - used,
1158                                    more_bits || i >= 0);
1159         }
1160       }
1161   }
1162
1163   /* NOTREACHED */
1164 }
1165 \f
1166 /* External user entry point.  */
1167
1168 FLOAT
1169 STRTOF (nptr, endptr)
1170      const STRING_TYPE *nptr;
1171      STRING_TYPE **endptr;
1172 {
1173   return INTERNAL (STRTOF) (nptr, endptr, 0);
1174 }
1175
1176 #define weak_this(x) weak_symbol(x)
1177 weak_this (STRTOF)