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