OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [strtod/] [strtod_l.c] - Blame information for rev 740

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn
/* Convert string representing a number to float value, using given locale.
2
   Copyright (C) 1997,1998,2002,2004,2005,2006,2007,2008,2009,2010
3
   Free Software Foundation, Inc.
4
   This file is part of the GNU C Library.
5
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
6
 
7
   The GNU C Library is free software; you can redistribute it and/or
8
   modify it under the terms of the GNU Lesser General Public
9
   License as published by the Free Software Foundation; either
10
   version 2.1 of the 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
   Lesser General Public License for more details.
16
 
17
   You should have received a copy of the GNU Lesser General Public
18
   License along with the GNU C Library; if not, write to the Free
19
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20
   02111-1307 USA.  */
21
 
22
#include <config.h>
23
#include <stdarg.h>
24
#include <string.h>
25
#include <float.h>
26
#include <math.h>
27
#define NDEBUG 1
28
#include <assert.h>
29
#ifdef HAVE_ERRNO_H
30
#include <errno.h>
31
#endif
32
#include "../printf/quadmath-printf.h"
33
#include "../printf/fpioconst.h"
34
 
35
 
36
#undef L_
37
#ifdef USE_WIDE_CHAR
38
# define STRING_TYPE wchar_t
39
# define CHAR_TYPE wint_t
40
# define L_(Ch) L##Ch
41
# define ISSPACE(Ch) __iswspace_l ((Ch), loc)
42
# define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
43
# define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
44
# define TOLOWER(Ch) __towlower_l ((Ch), loc)
45
# define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr)
46
# define STRNCASECMP(S1, S2, N) \
47
  __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
48
# define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
49
#else
50
# define STRING_TYPE char
51
# define CHAR_TYPE char
52
# define L_(Ch) Ch
53
# define ISSPACE(Ch) isspace (Ch)
54
# define ISDIGIT(Ch) isdigit (Ch)
55
# define ISXDIGIT(Ch) isxdigit (Ch)
56
# define TOLOWER(Ch) tolower (Ch)
57
# define TOLOWER_C(Ch) \
58
  ({__typeof(Ch) __tlc = (Ch); \
59
    (__tlc >= 'A' && __tlc <= 'Z') ? __tlc - 'A' + 'a' : __tlc; })
60
# define STRNCASECMP(S1, S2, N) \
61
  __quadmath_strncasecmp_c (S1, S2, N)
62
# ifdef HAVE_STRTOULL
63
#  define STRTOULL(S, E, B) strtoull (S, E, B)
64
# else
65
#  define STRTOULL(S, E, B) strtoul (S, E, B)
66
# endif
67
 
68
static inline int
69
__quadmath_strncasecmp_c (const char *s1, const char *s2, size_t n)
70
{
71
  const unsigned char *p1 = (const unsigned char *) s1;
72
  const unsigned char *p2 = (const unsigned char *) s2;
73
  int result;
74
  if (p1 == p2 || n == 0)
75
    return 0;
76
  while ((result = TOLOWER_C (*p1) - TOLOWER_C (*p2++)) == 0)
77
    if (*p1++ == '\0' || --n == 0)
78
      break;
79
 
80
  return result;
81
}
82
#endif
83
 
84
 
85
/* Constants we need from float.h; select the set for the FLOAT precision.  */
86
#define MANT_DIG        PASTE(FLT,_MANT_DIG)
87
#define DIG             PASTE(FLT,_DIG)
88
#define MAX_EXP         PASTE(FLT,_MAX_EXP)
89
#define MIN_EXP         PASTE(FLT,_MIN_EXP)
90
#define MAX_10_EXP      PASTE(FLT,_MAX_10_EXP)
91
#define MIN_10_EXP      PASTE(FLT,_MIN_10_EXP)
92
 
93
/* Extra macros required to get FLT expanded before the pasting.  */
94
#define PASTE(a,b)      PASTE1(a,b)
95
#define PASTE1(a,b)     a##b
96
 
97
/* Function to construct a floating point number from an MP integer
98
   containing the fraction bits, a base 2 exponent, and a sign flag.  */
99
extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
100
 
101
/* Definitions according to limb size used.  */
102
#if     BITS_PER_MP_LIMB == 32
103
# define MAX_DIG_PER_LIMB       9
104
# define MAX_FAC_PER_LIMB       1000000000UL
105
#elif   BITS_PER_MP_LIMB == 64
106
# define MAX_DIG_PER_LIMB       19
107
# define MAX_FAC_PER_LIMB       10000000000000000000ULL
108
#else
109
# error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
110
#endif
111
 
112
#define _tens_in_limb __quadmath_tens_in_limb
113
extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] attribute_hidden;
114
 
115
#ifndef howmany
116
#define howmany(x,y)            (((x)+((y)-1))/(y))
117
#endif
118
#define SWAP(x, y)              ({ typeof(x) _tmp = x; x = y; y = _tmp; })
119
 
120
#define NDIG                    (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
121
#define HEXNDIG                 ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
122
#define RETURN_LIMB_SIZE                howmany (MANT_DIG, BITS_PER_MP_LIMB)
123
 
124
#define RETURN(val,end)                                                       \
125
    do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end);                 \
126
         return val; } while (0)
127
 
128
/* Maximum size necessary for mpn integers to hold floating point numbers.  */
129
#define MPNSIZE         (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
130
                         + 2)
131
/* Declare an mpn integer variable that big.  */
132
#define MPN_VAR(name)   mp_limb_t name[MPNSIZE]; mp_size_t name##size
133
/* Copy an mpn integer value.  */
134
#define MPN_ASSIGN(dst, src) \
135
        memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
136
 
137
 
138
/* Return a floating point number of the needed type according to the given
139
   multi-precision number after possible rounding.  */
140
static FLOAT
141
round_and_return (mp_limb_t *retval, int exponent, int negative,
142
                  mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
143
{
144
  if (exponent < MIN_EXP - 1)
145
    {
146
      mp_size_t shift = MIN_EXP - 1 - exponent;
147
 
148
      if (shift > MANT_DIG)
149
        {
150
#if defined HAVE_ERRNO_H && defined EDOM
151
          errno = EDOM;
152
#endif
153
          return 0.0;
154
        }
155
 
156
      more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
157
      if (shift == MANT_DIG)
158
        /* This is a special case to handle the very seldom case where
159
           the mantissa will be empty after the shift.  */
160
        {
161
          int i;
162
 
163
          round_limb = retval[RETURN_LIMB_SIZE - 1];
164
          round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
165
          for (i = 0; i < RETURN_LIMB_SIZE; ++i)
166
            more_bits |= retval[i] != 0;
167
          MPN_ZERO (retval, RETURN_LIMB_SIZE);
168
        }
169
      else if (shift >= BITS_PER_MP_LIMB)
170
        {
171
          int i;
172
 
173
          round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
174
          round_bit = (shift - 1) % BITS_PER_MP_LIMB;
175
          for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
176
            more_bits |= retval[i] != 0;
177
          more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
178
                        != 0);
179
 
180
          (void) mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
181
                             RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
182
                             shift % BITS_PER_MP_LIMB);
183
          MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
184
                    shift / BITS_PER_MP_LIMB);
185
        }
186
      else if (shift > 0)
187
        {
188
          round_limb = retval[0];
189
          round_bit = shift - 1;
190
          (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
191
        }
192
      /* This is a hook for the m68k long double format, where the
193
         exponent bias is the same for normalized and denormalized
194
         numbers.  */
195
#ifndef DENORM_EXP
196
# define DENORM_EXP (MIN_EXP - 2)
197
#endif
198
      exponent = DENORM_EXP;
199
#if defined HAVE_ERRNO_H && defined ERANGE
200
      errno = ERANGE;
201
#endif
202
    }
203
 
204
  if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
205
      && (more_bits || (retval[0] & 1) != 0
206
          || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
207
    {
208
      mp_limb_t cy = mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
209
 
210
      if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
211
          ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
212
           (retval[RETURN_LIMB_SIZE - 1]
213
            & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
214
        {
215
          ++exponent;
216
          (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
217
          retval[RETURN_LIMB_SIZE - 1]
218
            |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
219
        }
220
      else if (exponent == DENORM_EXP
221
               && (retval[RETURN_LIMB_SIZE - 1]
222
                   & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
223
               != 0)
224
          /* The number was denormalized but now normalized.  */
225
        exponent = MIN_EXP - 1;
226
    }
227
 
228
  if (exponent > MAX_EXP)
229
    return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
230
 
231
  return MPN2FLOAT (retval, exponent, negative);
232
}
233
 
234
 
235
/* Read a multi-precision integer starting at STR with exactly DIGCNT digits
236
   into N.  Return the size of the number limbs in NSIZE at the first
237
   character od the string that is not part of the integer as the function
238
   value.  If the EXPONENT is small enough to be taken as an additional
239
   factor for the resulting number (see code) multiply by it.  */
240
static const STRING_TYPE *
241
str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
242
            int *exponent
243
#ifndef USE_WIDE_CHAR
244
            , const char *decimal, size_t decimal_len, const char *thousands
245
#endif
246
 
247
            )
248
{
249
  /* Number of digits for actual limb.  */
250
  int cnt = 0;
251
  mp_limb_t low = 0;
252
  mp_limb_t start;
253
 
254
  *nsize = 0;
255
  assert (digcnt > 0);
256
  do
257
    {
258
      if (cnt == MAX_DIG_PER_LIMB)
259
        {
260
          if (*nsize == 0)
261
            {
262
              n[0] = low;
263
              *nsize = 1;
264
            }
265
          else
266
            {
267
              mp_limb_t cy;
268
              cy = mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
269
              cy += mpn_add_1 (n, n, *nsize, low);
270
              if (cy != 0)
271
                {
272
                  n[*nsize] = cy;
273
                  ++(*nsize);
274
                }
275
            }
276
          cnt = 0;
277
          low = 0;
278
        }
279
 
280
      /* There might be thousands separators or radix characters in
281
         the string.  But these all can be ignored because we know the
282
         format of the number is correct and we have an exact number
283
         of characters to read.  */
284
#ifdef USE_WIDE_CHAR
285
      if (*str < L'0' || *str > L'9')
286
        ++str;
287
#else
288
      if (*str < '0' || *str > '9')
289
        {
290
          int inner = 0;
291
          if (thousands != NULL && *str == *thousands
292
              && ({ for (inner = 1; thousands[inner] != '\0'; ++inner)
293
                      if (thousands[inner] != str[inner])
294
                        break;
295
                    thousands[inner] == '\0'; }))
296
            str += inner;
297
          else
298
            str += decimal_len;
299
        }
300
#endif
301
      low = low * 10 + *str++ - L_('0');
302
      ++cnt;
303
    }
304
  while (--digcnt > 0);
305
 
306
  if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
307
    {
308
      low *= _tens_in_limb[*exponent];
309
      start = _tens_in_limb[cnt + *exponent];
310
      *exponent = 0;
311
    }
312
  else
313
    start = _tens_in_limb[cnt];
314
 
315
  if (*nsize == 0)
316
    {
317
      n[0] = low;
318
      *nsize = 1;
319
    }
320
  else
321
    {
322
      mp_limb_t cy;
323
      cy = mpn_mul_1 (n, n, *nsize, start);
324
      cy += mpn_add_1 (n, n, *nsize, low);
325
      if (cy != 0)
326
        n[(*nsize)++] = cy;
327
    }
328
 
329
  return str;
330
}
331
 
332
 
333
/* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
334
   with the COUNT most significant bits of LIMB.
335
 
336
   Tege doesn't like this function so I have to write it here myself. :)
337
   --drepper */
338
static inline void
339
__attribute ((always_inline))
340
mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count,
341
              mp_limb_t limb)
342
{
343
  if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB)
344
    {
345
      /* Optimize the case of shifting by exactly a word:
346
         just copy words, with no actual bit-shifting.  */
347
      mp_size_t i;
348
      for (i = size - 1; i > 0; --i)
349
        ptr[i] = ptr[i - 1];
350
      ptr[0] = limb;
351
    }
352
  else
353
    {
354
      (void) mpn_lshift (ptr, ptr, size, count);
355
      ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
356
    }
357
}
358
 
359
 
360
#define INTERNAL(x) INTERNAL1(x)
361
#define INTERNAL1(x) __##x##_internal
362
#ifndef ____STRTOF_INTERNAL
363
# define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
364
#endif
365
 
366
/* This file defines a function to check for correct grouping.  */
367
#include "grouping.h"
368
 
369
 
370
/* Return a floating point number with the value of the given string NPTR.
371
   Set *ENDPTR to the character after the last used one.  If the number is
372
   smaller than the smallest representable number, set `errno' to ERANGE and
373
   return 0.0.  If the number is too big to be represented, set `errno' to
374
   ERANGE and return HUGE_VAL with the appropriate sign.  */
375
FLOAT
376
____STRTOF_INTERNAL (nptr, endptr, group)
377
     const STRING_TYPE *nptr;
378
     STRING_TYPE **endptr;
379
     int group;
380
{
381
  int negative;                 /* The sign of the number.  */
382
  MPN_VAR (num);                /* MP representation of the number.  */
383
  int exponent;                 /* Exponent of the number.  */
384
 
385
  /* Numbers starting `0X' or `0x' have to be processed with base 16.  */
386
  int base = 10;
387
 
388
  /* When we have to compute fractional digits we form a fraction with a
389
     second multi-precision number (and we sometimes need a second for
390
     temporary results).  */
391
  MPN_VAR (den);
392
 
393
  /* Representation for the return value.  */
394
  mp_limb_t retval[RETURN_LIMB_SIZE];
395
  /* Number of bits currently in result value.  */
396
  int bits;
397
 
398
  /* Running pointer after the last character processed in the string.  */
399
  const STRING_TYPE *cp, *tp;
400
  /* Start of significant part of the number.  */
401
  const STRING_TYPE *startp, *start_of_digits;
402
  /* Points at the character following the integer and fractional digits.  */
403
  const STRING_TYPE *expp;
404
  /* Total number of digit and number of digits in integer part.  */
405
  int dig_no, int_no, lead_zero;
406
  /* Contains the last character read.  */
407
  CHAR_TYPE c;
408
 
409
  /* The radix character of the current locale.  */
410
#ifdef USE_WIDE_CHAR
411
  wchar_t decimal;
412
#else
413
  const char *decimal;
414
  size_t decimal_len;
415
#endif
416
  /* The thousands character of the current locale.  */
417
#ifdef USE_WIDE_CHAR
418
  wchar_t thousands = L'\0';
419
#else
420
  const char *thousands = NULL;
421
#endif
422
  /* The numeric grouping specification of the current locale,
423
     in the format described in <locale.h>.  */
424
  const char *grouping;
425
  /* Used in several places.  */
426
  int cnt;
427
 
428
#if defined USE_LOCALECONV && !defined USE_NL_LANGINFO
429
  const struct lconv *lc = localeconv ();
430
#endif
431
 
432
  if (__builtin_expect (group, 0))
433
    {
434
#ifdef USE_NL_LANGINFO
435
      grouping = nl_langinfo (GROUPING);
436
      if (*grouping <= 0 || *grouping == CHAR_MAX)
437
        grouping = NULL;
438
      else
439
        {
440
          /* Figure out the thousands separator character.  */
441
#ifdef USE_WIDE_CHAR
442
          thousands = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC);
443
          if (thousands == L'\0')
444
            grouping = NULL;
445
#else
446
          thousands = nl_langinfo (THOUSANDS_SEP);
447
          if (*thousands == '\0')
448
            {
449
              thousands = NULL;
450
              grouping = NULL;
451
            }
452
#endif
453
        }
454
#elif defined USE_LOCALECONV
455
      grouping = lc->grouping;
456
      if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX)
457
        grouping = NULL;
458
      else
459
        {
460
          /* Figure out the thousands separator character.  */
461
          thousands = lc->thousands_sep;
462
          if (thousands == NULL || *thousands == '\0')
463
            {
464
              thousands = NULL;
465
              grouping = NULL;
466
            }
467
        }
468
#else
469
      grouping = NULL;
470
#endif
471
    }
472
  else
473
    grouping = NULL;
474
 
475
  /* Find the locale's decimal point character.  */
476
#ifdef USE_WIDE_CHAR
477
  decimal = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
478
  assert (decimal != L'\0');
479
# define decimal_len 1
480
#else
481
#ifdef USE_NL_LANGINFO
482
  decimal = nl_langinfo (DECIMAL_POINT);
483
  decimal_len = strlen (decimal);
484
  assert (decimal_len > 0);
485
#elif defined USE_LOCALECONV
486
  decimal = lc->decimal_point;
487
  if (decimal == NULL || *decimal == '\0')
488
    decimal = ".";
489
  decimal_len = strlen (decimal);
490
#else
491
  decimal = ".";
492
  decimal_len = 1;
493
#endif
494
#endif
495
 
496
  /* Prepare number representation.  */
497
  exponent = 0;
498
  negative = 0;
499
  bits = 0;
500
 
501
  /* Parse string to get maximal legal prefix.  We need the number of
502
     characters of the integer part, the fractional part and the exponent.  */
503
  cp = nptr - 1;
504
  /* Ignore leading white space.  */
505
  do
506
    c = *++cp;
507
  while (ISSPACE (c));
508
 
509
  /* Get sign of the result.  */
510
  if (c == L_('-'))
511
    {
512
      negative = 1;
513
      c = *++cp;
514
    }
515
  else if (c == L_('+'))
516
    c = *++cp;
517
 
518
  /* Return 0.0 if no legal string is found.
519
     No character is used even if a sign was found.  */
520
#ifdef USE_WIDE_CHAR
521
  if (c == (wint_t) decimal
522
      && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9')
523
    {
524
      /* We accept it.  This funny construct is here only to indent
525
         the code correctly.  */
526
    }
527
#else
528
  for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
529
    if (cp[cnt] != decimal[cnt])
530
      break;
531
  if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9')
532
    {
533
      /* We accept it.  This funny construct is here only to indent
534
         the code correctly.  */
535
    }
536
#endif
537
  else if (c < L_('0') || c > L_('9'))
538
    {
539
      /* Check for `INF' or `INFINITY'.  */
540
      CHAR_TYPE lowc = TOLOWER_C (c);
541
 
542
      if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0)
543
        {
544
          /* Return +/- infinity.  */
545
          if (endptr != NULL)
546
            *endptr = (STRING_TYPE *)
547
                      (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0
548
                             ? 8 : 3));
549
 
550
          return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
551
        }
552
 
553
      if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
554
        {
555
          /* Return NaN.  */
556
          FLOAT retval = NAN;
557
 
558
          cp += 3;
559
 
560
          /* Match `(n-char-sequence-digit)'.  */
561
          if (*cp == L_('('))
562
            {
563
              const STRING_TYPE *startp = cp;
564
              do
565
                ++cp;
566
              while ((*cp >= L_('0') && *cp <= L_('9'))
567
                     || ({ CHAR_TYPE lo = TOLOWER (*cp);
568
                           lo >= L_('a') && lo <= L_('z'); })
569
                     || *cp == L_('_'));
570
 
571
              if (*cp != L_(')'))
572
                /* The closing brace is missing.  Only match the NAN
573
                   part.  */
574
                cp = startp;
575
              else
576
                {
577
                  /* This is a system-dependent way to specify the
578
                     bitmask used for the NaN.  We expect it to be
579
                     a number which is put in the mantissa of the
580
                     number.  */
581
                  STRING_TYPE *endp;
582
                  unsigned long long int mant;
583
 
584
                  mant = STRTOULL (startp + 1, &endp, 0);
585
                  if (endp == cp)
586
                    SET_MANTISSA (retval, mant);
587
 
588
                  /* Consume the closing brace.  */
589
                  ++cp;
590
                }
591
            }
592
 
593
          if (endptr != NULL)
594
            *endptr = (STRING_TYPE *) cp;
595
 
596
          return retval;
597
        }
598
 
599
      /* It is really a text we do not recognize.  */
600
      RETURN (0.0, nptr);
601
    }
602
 
603
  /* First look whether we are faced with a hexadecimal number.  */
604
  if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
605
    {
606
      /* Okay, it is a hexa-decimal number.  Remember this and skip
607
         the characters.  BTW: hexadecimal numbers must not be
608
         grouped.  */
609
      base = 16;
610
      cp += 2;
611
      c = *cp;
612
      grouping = NULL;
613
    }
614
 
615
  /* Record the start of the digits, in case we will check their grouping.  */
616
  start_of_digits = startp = cp;
617
 
618
  /* Ignore leading zeroes.  This helps us to avoid useless computations.  */
619
#ifdef USE_WIDE_CHAR
620
  while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands))
621
    c = *++cp;
622
#else
623
  if (__builtin_expect (thousands == NULL, 1))
624
    while (c == '0')
625
      c = *++cp;
626
  else
627
    {
628
      /* We also have the multibyte thousands string.  */
629
      while (1)
630
        {
631
          if (c != '0')
632
            {
633
              for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
634
                if (thousands[cnt] != cp[cnt])
635
                  break;
636
              if (thousands[cnt] != '\0')
637
                break;
638
              cp += cnt - 1;
639
            }
640
          c = *++cp;
641
        }
642
    }
643
#endif
644
 
645
  /* If no other digit but a '0' is found the result is 0.0.
646
     Return current read pointer.  */
647
  CHAR_TYPE lowc = TOLOWER (c);
648
  if (!((c >= L_('0') && c <= L_('9'))
649
        || (base == 16 && lowc >= L_('a') && lowc <= L_('f'))
650
        || (
651
#ifdef USE_WIDE_CHAR
652
            c == (wint_t) decimal
653
#else
654
            ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
655
                 if (decimal[cnt] != cp[cnt])
656
                   break;
657
               decimal[cnt] == '\0'; })
658
#endif
659
            /* '0x.' alone is not a valid hexadecimal number.
660
               '.' alone is not valid either, but that has been checked
661
               already earlier.  */
662
            && (base != 16
663
                || cp != start_of_digits
664
                || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9'))
665
                || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]);
666
                      lo >= L_('a') && lo <= L_('f'); })))
667
        || (base == 16 && (cp != start_of_digits
668
                           && lowc == L_('p')))
669
        || (base != 16 && lowc == L_('e'))))
670
    {
671
#ifdef USE_WIDE_CHAR
672
      tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
673
                                         grouping);
674
#else
675
      tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
676
                                         grouping);
677
#endif
678
      /* If TP is at the start of the digits, there was no correctly
679
         grouped prefix of the string; so no number found.  */
680
      RETURN (negative ? -0.0 : 0.0,
681
              tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
682
    }
683
 
684
  /* Remember first significant digit and read following characters until the
685
     decimal point, exponent character or any non-FP number character.  */
686
  startp = cp;
687
  dig_no = 0;
688
  while (1)
689
    {
690
      if ((c >= L_('0') && c <= L_('9'))
691
          || (base == 16
692
              && ({ CHAR_TYPE lo = TOLOWER (c);
693
                    lo >= L_('a') && lo <= L_('f'); })))
694
        ++dig_no;
695
      else
696
        {
697
#ifdef USE_WIDE_CHAR
698
          if (__builtin_expect ((wint_t) thousands == L'\0', 1)
699
              || c != (wint_t) thousands)
700
            /* Not a digit or separator: end of the integer part.  */
701
            break;
702
#else
703
          if (__builtin_expect (thousands == NULL, 1))
704
            break;
705
          else
706
            {
707
              for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
708
                if (thousands[cnt] != cp[cnt])
709
                  break;
710
              if (thousands[cnt] != '\0')
711
                break;
712
              cp += cnt - 1;
713
            }
714
#endif
715
        }
716
      c = *++cp;
717
    }
718
 
719
  if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits)
720
    {
721
      /* Check the grouping of the digits.  */
722
#ifdef USE_WIDE_CHAR
723
      tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
724
                                         grouping);
725
#else
726
      tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
727
                                         grouping);
728
#endif
729
      if (cp != tp)
730
        {
731
          /* Less than the entire string was correctly grouped.  */
732
 
733
          if (tp == start_of_digits)
734
            /* No valid group of numbers at all: no valid number.  */
735
            RETURN (0.0, nptr);
736
 
737
          if (tp < startp)
738
            /* The number is validly grouped, but consists
739
               only of zeroes.  The whole value is zero.  */
740
            RETURN (negative ? -0.0 : 0.0, tp);
741
 
742
          /* Recompute DIG_NO so we won't read more digits than
743
             are properly grouped.  */
744
          cp = tp;
745
          dig_no = 0;
746
          for (tp = startp; tp < cp; ++tp)
747
            if (*tp >= L_('0') && *tp <= L_('9'))
748
              ++dig_no;
749
 
750
          int_no = dig_no;
751
          lead_zero = 0;
752
 
753
          goto number_parsed;
754
        }
755
    }
756
 
757
  /* We have the number of digits in the integer part.  Whether these
758
     are all or any is really a fractional digit will be decided
759
     later.  */
760
  int_no = dig_no;
761
  lead_zero = int_no == 0 ? -1 : 0;
762
 
763
  /* Read the fractional digits.  A special case are the 'american
764
     style' numbers like `16.' i.e. with decimal point but without
765
     trailing digits.  */
766
  if (
767
#ifdef USE_WIDE_CHAR
768
      c == (wint_t) decimal
769
#else
770
      ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
771
           if (decimal[cnt] != cp[cnt])
772
             break;
773
         decimal[cnt] == '\0'; })
774
#endif
775
      )
776
    {
777
      cp += decimal_len;
778
      c = *cp;
779
      while ((c >= L_('0') && c <= L_('9')) ||
780
             (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c);
781
                               lo >= L_('a') && lo <= L_('f'); })))
782
        {
783
          if (c != L_('0') && lead_zero == -1)
784
            lead_zero = dig_no - int_no;
785
          ++dig_no;
786
          c = *++cp;
787
        }
788
    }
789
 
790
  /* Remember start of exponent (if any).  */
791
  expp = cp;
792
 
793
  /* Read exponent.  */
794
  lowc = TOLOWER (c);
795
  if ((base == 16 && lowc == L_('p'))
796
      || (base != 16 && lowc == L_('e')))
797
    {
798
      int exp_negative = 0;
799
 
800
      c = *++cp;
801
      if (c == L_('-'))
802
        {
803
          exp_negative = 1;
804
          c = *++cp;
805
        }
806
      else if (c == L_('+'))
807
        c = *++cp;
808
 
809
      if (c >= L_('0') && c <= L_('9'))
810
        {
811
          int exp_limit;
812
 
813
          /* Get the exponent limit. */
814
          if (base == 16)
815
            exp_limit = (exp_negative ?
816
                         -MIN_EXP + MANT_DIG + 4 * int_no :
817
                         MAX_EXP - 4 * int_no + 4 * lead_zero + 3);
818
          else
819
            exp_limit = (exp_negative ?
820
                         -MIN_10_EXP + MANT_DIG + int_no :
821
                         MAX_10_EXP - int_no + lead_zero + 1);
822
 
823
          do
824
            {
825
              exponent *= 10;
826
              exponent += c - L_('0');
827
 
828
              if (__builtin_expect (exponent > exp_limit, 0))
829
                /* The exponent is too large/small to represent a valid
830
                   number.  */
831
                {
832
                  FLOAT result;
833
 
834
                  /* We have to take care for special situation: a joker
835
                     might have written "0.0e100000" which is in fact
836
                     zero.  */
837
                  if (lead_zero == -1)
838
                    result = negative ? -0.0 : 0.0;
839
                  else
840
                    {
841
                      /* Overflow or underflow.  */
842
#if defined HAVE_ERRNO_H && defined ERANGE
843
                      errno = ERANGE;
844
#endif
845
                      result = (exp_negative ? (negative ? -0.0 : 0.0) :
846
                                negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
847
                    }
848
 
849
                  /* Accept all following digits as part of the exponent.  */
850
                  do
851
                    ++cp;
852
                  while (*cp >= L_('0') && *cp <= L_('9'));
853
 
854
                  RETURN (result, cp);
855
                  /* NOTREACHED */
856
                }
857
 
858
              c = *++cp;
859
            }
860
          while (c >= L_('0') && c <= L_('9'));
861
 
862
          if (exp_negative)
863
            exponent = -exponent;
864
        }
865
      else
866
        cp = expp;
867
    }
868
 
869
  /* We don't want to have to work with trailing zeroes after the radix.  */
870
  if (dig_no > int_no)
871
    {
872
      while (expp[-1] == L_('0'))
873
        {
874
          --expp;
875
          --dig_no;
876
        }
877
      assert (dig_no >= int_no);
878
    }
879
 
880
  if (dig_no == int_no && dig_no > 0 && exponent < 0)
881
    do
882
      {
883
        while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1])))
884
          --expp;
885
 
886
        if (expp[-1] != L_('0'))
887
          break;
888
 
889
        --expp;
890
        --dig_no;
891
        --int_no;
892
        exponent += base == 16 ? 4 : 1;
893
      }
894
    while (dig_no > 0 && exponent < 0);
895
 
896
 number_parsed:
897
 
898
  /* The whole string is parsed.  Store the address of the next character.  */
899
  if (endptr)
900
    *endptr = (STRING_TYPE *) cp;
901
 
902
  if (dig_no == 0)
903
    return negative ? -0.0 : 0.0;
904
 
905
  if (lead_zero)
906
    {
907
      /* Find the decimal point */
908
#ifdef USE_WIDE_CHAR
909
      while (*startp != decimal)
910
        ++startp;
911
#else
912
      while (1)
913
        {
914
          if (*startp == decimal[0])
915
            {
916
              for (cnt = 1; decimal[cnt] != '\0'; ++cnt)
917
                if (decimal[cnt] != startp[cnt])
918
                  break;
919
              if (decimal[cnt] == '\0')
920
                break;
921
            }
922
          ++startp;
923
        }
924
#endif
925
      startp += lead_zero + decimal_len;
926
      exponent -= base == 16 ? 4 * lead_zero : lead_zero;
927
      dig_no -= lead_zero;
928
    }
929
 
930
  /* If the BASE is 16 we can use a simpler algorithm.  */
931
  if (base == 16)
932
    {
933
      static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
934
                                     4, 4, 4, 4, 4, 4, 4, 4 };
935
      int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
936
      int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
937
      mp_limb_t val;
938
 
939
      while (!ISXDIGIT (*startp))
940
        ++startp;
941
      while (*startp == L_('0'))
942
        ++startp;
943
      if (ISDIGIT (*startp))
944
        val = *startp++ - L_('0');
945
      else
946
        val = 10 + TOLOWER (*startp++) - L_('a');
947
      bits = nbits[val];
948
      /* We cannot have a leading zero.  */
949
      assert (bits != 0);
950
 
951
      if (pos + 1 >= 4 || pos + 1 >= bits)
952
        {
953
          /* We don't have to care for wrapping.  This is the normal
954
             case so we add the first clause in the `if' expression as
955
             an optimization.  It is a compile-time constant and so does
956
             not cost anything.  */
957
          retval[idx] = val << (pos - bits + 1);
958
          pos -= bits;
959
        }
960
      else
961
        {
962
          retval[idx--] = val >> (bits - pos - 1);
963
          retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
964
          pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
965
        }
966
 
967
      /* Adjust the exponent for the bits we are shifting in.  */
968
      exponent += bits - 1 + (int_no - 1) * 4;
969
 
970
      while (--dig_no > 0 && idx >= 0)
971
        {
972
          if (!ISXDIGIT (*startp))
973
            startp += decimal_len;
974
          if (ISDIGIT (*startp))
975
            val = *startp++ - L_('0');
976
          else
977
            val = 10 + TOLOWER (*startp++) - L_('a');
978
 
979
          if (pos + 1 >= 4)
980
            {
981
              retval[idx] |= val << (pos - 4 + 1);
982
              pos -= 4;
983
            }
984
          else
985
            {
986
              retval[idx--] |= val >> (4 - pos - 1);
987
              val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
988
              if (idx < 0)
989
                return round_and_return (retval, exponent, negative, val,
990
                                         BITS_PER_MP_LIMB - 1, dig_no > 0);
991
 
992
              retval[idx] = val;
993
              pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
994
            }
995
        }
996
 
997
      /* We ran out of digits.  */
998
      MPN_ZERO (retval, idx);
999
 
1000
      return round_and_return (retval, exponent, negative, 0, 0, 0);
1001
    }
1002
 
1003
  /* Now we have the number of digits in total and the integer digits as well
1004
     as the exponent and its sign.  We can decide whether the read digits are
1005
     really integer digits or belong to the fractional part; i.e. we normalize
1006
     123e-2 to 1.23.  */
1007
  {
1008
    register int incr = (exponent < 0 ? MAX (-int_no, exponent)
1009
                         : MIN (dig_no - int_no, exponent));
1010
    int_no += incr;
1011
    exponent -= incr;
1012
  }
1013
 
1014
  if (__builtin_expect (int_no + exponent > MAX_10_EXP + 1, 0))
1015
    {
1016
#if defined HAVE_ERRNO_H && defined ERANGE
1017
      errno = ERANGE;
1018
#endif
1019
      return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
1020
    }
1021
 
1022
  if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0))
1023
    {
1024
#if defined HAVE_ERRNO_H && defined ERANGE
1025
      errno = ERANGE;
1026
#endif
1027
      return negative ? -0.0 : 0.0;
1028
    }
1029
 
1030
  if (int_no > 0)
1031
    {
1032
      /* Read the integer part as a multi-precision number to NUM.  */
1033
      startp = str_to_mpn (startp, int_no, num, &numsize, &exponent
1034
#ifndef USE_WIDE_CHAR
1035
                           , decimal, decimal_len, thousands
1036
#endif
1037
                           );
1038
 
1039
      if (exponent > 0)
1040
        {
1041
          /* We now multiply the gained number by the given power of ten.  */
1042
          mp_limb_t *psrc = num;
1043
          mp_limb_t *pdest = den;
1044
          int expbit = 1;
1045
          const struct mp_power *ttab = &_fpioconst_pow10[0];
1046
 
1047
          do
1048
            {
1049
              if ((exponent & expbit) != 0)
1050
                {
1051
                  size_t size = ttab->arraysize - _FPIO_CONST_OFFSET;
1052
                  mp_limb_t cy;
1053
                  exponent ^= expbit;
1054
 
1055
                  /* FIXME: not the whole multiplication has to be
1056
                     done.  If we have the needed number of bits we
1057
                     only need the information whether more non-zero
1058
                     bits follow.  */
1059
                  if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
1060
                    cy = mpn_mul (pdest, psrc, numsize,
1061
                                  &__tens[ttab->arrayoff
1062
                                          + _FPIO_CONST_OFFSET],
1063
                                  size);
1064
                  else
1065
                    cy = mpn_mul (pdest, &__tens[ttab->arrayoff
1066
                                                 + _FPIO_CONST_OFFSET],
1067
                                  size, psrc, numsize);
1068
                  numsize += size;
1069
                  if (cy == 0)
1070
                    --numsize;
1071
                  (void) SWAP (psrc, pdest);
1072
                }
1073
              expbit <<= 1;
1074
              ++ttab;
1075
            }
1076
          while (exponent != 0);
1077
 
1078
          if (psrc == den)
1079
            memcpy (num, den, numsize * sizeof (mp_limb_t));
1080
        }
1081
 
1082
      /* Determine how many bits of the result we already have.  */
1083
      count_leading_zeros (bits, num[numsize - 1]);
1084
      bits = numsize * BITS_PER_MP_LIMB - bits;
1085
 
1086
      /* Now we know the exponent of the number in base two.
1087
         Check it against the maximum possible exponent.  */
1088
      if (__builtin_expect (bits > MAX_EXP, 0))
1089
        {
1090
#if defined HAVE_ERRNO_H && defined ERANGE
1091
          errno = ERANGE;
1092
#endif
1093
          return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
1094
        }
1095
 
1096
      /* We have already the first BITS bits of the result.  Together with
1097
         the information whether more non-zero bits follow this is enough
1098
         to determine the result.  */
1099
      if (bits > MANT_DIG)
1100
        {
1101
          int i;
1102
          const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
1103
          const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
1104
          const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
1105
                                                     : least_idx;
1106
          const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
1107
                                                     : least_bit - 1;
1108
 
1109
          if (least_bit == 0)
1110
            memcpy (retval, &num[least_idx],
1111
                    RETURN_LIMB_SIZE * sizeof (mp_limb_t));
1112
          else
1113
            {
1114
              for (i = least_idx; i < numsize - 1; ++i)
1115
                retval[i - least_idx] = (num[i] >> least_bit)
1116
                                        | (num[i + 1]
1117
                                           << (BITS_PER_MP_LIMB - least_bit));
1118
              if (i - least_idx < RETURN_LIMB_SIZE)
1119
                retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
1120
            }
1121
 
1122
          /* Check whether any limb beside the ones in RETVAL are non-zero.  */
1123
          for (i = 0; num[i] == 0; ++i)
1124
            ;
1125
 
1126
          return round_and_return (retval, bits - 1, negative,
1127
                                   num[round_idx], round_bit,
1128
                                   int_no < dig_no || i < round_idx);
1129
          /* NOTREACHED */
1130
        }
1131
      else if (dig_no == int_no)
1132
        {
1133
          const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
1134
          const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
1135
 
1136
          if (target_bit == is_bit)
1137
            {
1138
              memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
1139
                      numsize * sizeof (mp_limb_t));
1140
              /* FIXME: the following loop can be avoided if we assume a
1141
                 maximal MANT_DIG value.  */
1142
              MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
1143
            }
1144
          else if (target_bit > is_bit)
1145
            {
1146
              (void) mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
1147
                                 num, numsize, target_bit - is_bit);
1148
              /* FIXME: the following loop can be avoided if we assume a
1149
                 maximal MANT_DIG value.  */
1150
              MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
1151
            }
1152
          else
1153
            {
1154
              mp_limb_t cy;
1155
              assert (numsize < RETURN_LIMB_SIZE);
1156
 
1157
              cy = mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
1158
                               num, numsize, is_bit - target_bit);
1159
              retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
1160
              /* FIXME: the following loop can be avoided if we assume a
1161
                 maximal MANT_DIG value.  */
1162
              MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
1163
            }
1164
 
1165
          return round_and_return (retval, bits - 1, negative, 0, 0, 0);
1166
          /* NOTREACHED */
1167
        }
1168
 
1169
      /* Store the bits we already have.  */
1170
      memcpy (retval, num, numsize * sizeof (mp_limb_t));
1171
#if RETURN_LIMB_SIZE > 1
1172
      if (numsize < RETURN_LIMB_SIZE)
1173
# if RETURN_LIMB_SIZE == 2
1174
        retval[numsize] = 0;
1175
# else
1176
        MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize);
1177
# endif
1178
#endif
1179
    }
1180
 
1181
  /* We have to compute at least some of the fractional digits.  */
1182
  {
1183
    /* We construct a fraction and the result of the division gives us
1184
       the needed digits.  The denominator is 1.0 multiplied by the
1185
       exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1186
       123e-6 gives 123 / 1000000.  */
1187
 
1188
    int expbit;
1189
    int neg_exp;
1190
    int more_bits;
1191
    mp_limb_t cy;
1192
    mp_limb_t *psrc = den;
1193
    mp_limb_t *pdest = num;
1194
    const struct mp_power *ttab = &_fpioconst_pow10[0];
1195
 
1196
    assert (dig_no > int_no && exponent <= 0);
1197
 
1198
 
1199
    /* For the fractional part we need not process too many digits.  One
1200
       decimal digits gives us log_2(10) ~ 3.32 bits.  If we now compute
1201
                        ceil(BITS / 3) =: N
1202
       digits we should have enough bits for the result.  The remaining
1203
       decimal digits give us the information that more bits are following.
1204
       This can be used while rounding.  (Two added as a safety margin.)  */
1205
    if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 2)
1206
      {
1207
        dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 2;
1208
        more_bits = 1;
1209
      }
1210
    else
1211
      more_bits = 0;
1212
 
1213
    neg_exp = dig_no - int_no - exponent;
1214
 
1215
    /* Construct the denominator.  */
1216
    densize = 0;
1217
    expbit = 1;
1218
    do
1219
      {
1220
        if ((neg_exp & expbit) != 0)
1221
          {
1222
            mp_limb_t cy;
1223
            neg_exp ^= expbit;
1224
 
1225
            if (densize == 0)
1226
              {
1227
                densize = ttab->arraysize - _FPIO_CONST_OFFSET;
1228
                memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET],
1229
                        densize * sizeof (mp_limb_t));
1230
              }
1231
            else
1232
              {
1233
                cy = mpn_mul (pdest, &__tens[ttab->arrayoff
1234
                                             + _FPIO_CONST_OFFSET],
1235
                              ttab->arraysize - _FPIO_CONST_OFFSET,
1236
                              psrc, densize);
1237
                densize += ttab->arraysize - _FPIO_CONST_OFFSET;
1238
                if (cy == 0)
1239
                  --densize;
1240
                (void) SWAP (psrc, pdest);
1241
              }
1242
          }
1243
        expbit <<= 1;
1244
        ++ttab;
1245
      }
1246
    while (neg_exp != 0);
1247
 
1248
    if (psrc == num)
1249
      memcpy (den, num, densize * sizeof (mp_limb_t));
1250
 
1251
    /* Read the fractional digits from the string.  */
1252
    (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent
1253
#ifndef USE_WIDE_CHAR
1254
                       , decimal, decimal_len, thousands
1255
#endif
1256
                       );
1257
 
1258
    /* We now have to shift both numbers so that the highest bit in the
1259
       denominator is set.  In the same process we copy the numerator to
1260
       a high place in the array so that the division constructs the wanted
1261
       digits.  This is done by a "quasi fix point" number representation.
1262
 
1263
       num:   ddddddddddd . 0000000000000000000000
1264
              |--- m ---|
1265
       den:                            ddddddddddd      n >= m
1266
                                       |--- n ---|
1267
     */
1268
 
1269
    count_leading_zeros (cnt, den[densize - 1]);
1270
 
1271
    if (cnt > 0)
1272
      {
1273
        /* Don't call `mpn_shift' with a count of zero since the specification
1274
           does not allow this.  */
1275
        (void) mpn_lshift (den, den, densize, cnt);
1276
        cy = mpn_lshift (num, num, numsize, cnt);
1277
        if (cy != 0)
1278
          num[numsize++] = cy;
1279
      }
1280
 
1281
    /* Now we are ready for the division.  But it is not necessary to
1282
       do a full multi-precision division because we only need a small
1283
       number of bits for the result.  So we do not use mpn_divmod
1284
       here but instead do the division here by hand and stop whenever
1285
       the needed number of bits is reached.  The code itself comes
1286
       from the GNU MP Library by Torbj\"orn Granlund.  */
1287
 
1288
    exponent = bits;
1289
 
1290
    switch (densize)
1291
      {
1292
      case 1:
1293
        {
1294
          mp_limb_t d, n, quot;
1295
          int used = 0;
1296
 
1297
          n = num[0];
1298
          d = den[0];
1299
          assert (numsize == 1 && n < d);
1300
 
1301
          do
1302
            {
1303
              udiv_qrnnd (quot, n, n, 0, d);
1304
 
1305
#define got_limb                                                              \
1306
              if (bits == 0)                                                   \
1307
                {                                                             \
1308
                  register int cnt;                                           \
1309
                  if (quot == 0)                                       \
1310
                    cnt = BITS_PER_MP_LIMB;                                   \
1311
                  else                                                        \
1312
                    count_leading_zeros (cnt, quot);                          \
1313
                  exponent -= cnt;                                            \
1314
                  if (BITS_PER_MP_LIMB - cnt > MANT_DIG)                      \
1315
                    {                                                         \
1316
                      used = MANT_DIG + cnt;                                  \
1317
                      retval[0] = quot >> (BITS_PER_MP_LIMB - used);           \
1318
                      bits = MANT_DIG + 1;                                    \
1319
                    }                                                         \
1320
                  else                                                        \
1321
                    {                                                         \
1322
                      /* Note that we only clear the second element.  */      \
1323
                      /* The conditional is determined at compile time.  */   \
1324
                      if (RETURN_LIMB_SIZE > 1)                               \
1325
                        retval[1] = 0;                                         \
1326
                      retval[0] = quot;                                        \
1327
                      bits = -cnt;                                            \
1328
                    }                                                         \
1329
                }                                                             \
1330
              else if (bits + BITS_PER_MP_LIMB <= MANT_DIG)                   \
1331
                mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB,     \
1332
                              quot);                                          \
1333
              else                                                            \
1334
                {                                                             \
1335
                  used = MANT_DIG - bits;                                     \
1336
                  if (used > 0)                                                \
1337
                    mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot);      \
1338
                }                                                             \
1339
              bits += BITS_PER_MP_LIMB
1340
 
1341
              got_limb;
1342
            }
1343
          while (bits <= MANT_DIG);
1344
 
1345
          return round_and_return (retval, exponent - 1, negative,
1346
                                   quot, BITS_PER_MP_LIMB - 1 - used,
1347
                                   more_bits || n != 0);
1348
        }
1349
      case 2:
1350
        {
1351
          mp_limb_t d0, d1, n0, n1;
1352
          mp_limb_t quot = 0;
1353
          int used = 0;
1354
 
1355
          d0 = den[0];
1356
          d1 = den[1];
1357
 
1358
          if (numsize < densize)
1359
            {
1360
              if (num[0] >= d1)
1361
                {
1362
                  /* The numerator of the number occupies fewer bits than
1363
                     the denominator but the one limb is bigger than the
1364
                     high limb of the numerator.  */
1365
                  n1 = 0;
1366
                  n0 = num[0];
1367
                }
1368
              else
1369
                {
1370
                  if (bits <= 0)
1371
                    exponent -= BITS_PER_MP_LIMB;
1372
                  else
1373
                    {
1374
                      if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
1375
                        mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1376
                                      BITS_PER_MP_LIMB, 0);
1377
                      else
1378
                        {
1379
                          used = MANT_DIG - bits;
1380
                          if (used > 0)
1381
                            mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1382
                        }
1383
                      bits += BITS_PER_MP_LIMB;
1384
                    }
1385
                  n1 = num[0];
1386
                  n0 = 0;
1387
                }
1388
            }
1389
          else
1390
            {
1391
              n1 = num[1];
1392
              n0 = num[0];
1393
            }
1394
 
1395
          while (bits <= MANT_DIG)
1396
            {
1397
              mp_limb_t r;
1398
 
1399
              if (n1 == d1)
1400
                {
1401
                  /* QUOT should be either 111..111 or 111..110.  We need
1402
                     special treatment of this rare case as normal division
1403
                     would give overflow.  */
1404
                  quot = ~(mp_limb_t) 0;
1405
 
1406
                  r = n0 + d1;
1407
                  if (r < d1)   /* Carry in the addition?  */
1408
                    {
1409
                      add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
1410
                      goto have_quot;
1411
                    }
1412
                  n1 = d0 - (d0 != 0);
1413
                  n0 = -d0;
1414
                }
1415
              else
1416
                {
1417
                  udiv_qrnnd (quot, r, n1, n0, d1);
1418
                  umul_ppmm (n1, n0, d0, quot);
1419
                }
1420
 
1421
            q_test:
1422
              if (n1 > r || (n1 == r && n0 > 0))
1423
                {
1424
                  /* The estimated QUOT was too large.  */
1425
                  --quot;
1426
 
1427
                  sub_ddmmss (n1, n0, n1, n0, 0, d0);
1428
                  r += d1;
1429
                  if (r >= d1)  /* If not carry, test QUOT again.  */
1430
                    goto q_test;
1431
                }
1432
              sub_ddmmss (n1, n0, r, 0, n1, n0);
1433
 
1434
            have_quot:
1435
              got_limb;
1436
            }
1437
 
1438
          return round_and_return (retval, exponent - 1, negative,
1439
                                   quot, BITS_PER_MP_LIMB - 1 - used,
1440
                                   more_bits || n1 != 0 || n0 != 0);
1441
        }
1442
      default:
1443
        {
1444
          int i;
1445
          mp_limb_t cy, dX, d1, n0, n1;
1446
          mp_limb_t quot = 0;
1447
          int used = 0;
1448
 
1449
          dX = den[densize - 1];
1450
          d1 = den[densize - 2];
1451
 
1452
          /* The division does not work if the upper limb of the two-limb
1453
             numerator is greater than the denominator.  */
1454
          if (mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1455
            num[numsize++] = 0;
1456
 
1457
          if (numsize < densize)
1458
            {
1459
              mp_size_t empty = densize - numsize;
1460
              register int i;
1461
 
1462
              if (bits <= 0)
1463
                exponent -= empty * BITS_PER_MP_LIMB;
1464
              else
1465
                {
1466
                  if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1467
                    {
1468
                      /* We make a difference here because the compiler
1469
                         cannot optimize the `else' case that good and
1470
                         this reflects all currently used FLOAT types
1471
                         and GMP implementations.  */
1472
#if RETURN_LIMB_SIZE <= 2
1473
                      assert (empty == 1);
1474
                      mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1475
                                    BITS_PER_MP_LIMB, 0);
1476
#else
1477
                      for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i)
1478
                        retval[i] = retval[i - empty];
1479
                      while (i >= 0)
1480
                        retval[i--] = 0;
1481
#endif
1482
                    }
1483
                  else
1484
                    {
1485
                      used = MANT_DIG - bits;
1486
                      if (used >= BITS_PER_MP_LIMB)
1487
                        {
1488
                          register int i;
1489
                          (void) mpn_lshift (&retval[used
1490
                                                     / BITS_PER_MP_LIMB],
1491
                                             retval,
1492
                                             (RETURN_LIMB_SIZE
1493
                                              - used / BITS_PER_MP_LIMB),
1494
                                             used % BITS_PER_MP_LIMB);
1495
                          for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i)
1496
                            retval[i] = 0;
1497
                        }
1498
                      else if (used > 0)
1499
                        mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1500
                    }
1501
                  bits += empty * BITS_PER_MP_LIMB;
1502
                }
1503
              for (i = numsize; i > 0; --i)
1504
                num[i + empty] = num[i - 1];
1505
              MPN_ZERO (num, empty + 1);
1506
            }
1507
          else
1508
            {
1509
              int i;
1510
              assert (numsize == densize);
1511
              for (i = numsize; i > 0; --i)
1512
                num[i] = num[i - 1];
1513
            }
1514
 
1515
          den[densize] = 0;
1516
          n0 = num[densize];
1517
 
1518
          while (bits <= MANT_DIG)
1519
            {
1520
              if (n0 == dX)
1521
                /* This might over-estimate QUOT, but it's probably not
1522
                   worth the extra code here to find out.  */
1523
                quot = ~(mp_limb_t) 0;
1524
              else
1525
                {
1526
                  mp_limb_t r;
1527
 
1528
                  udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1529
                  umul_ppmm (n1, n0, d1, quot);
1530
 
1531
                  while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1532
                    {
1533
                      --quot;
1534
                      r += dX;
1535
                      if (r < dX) /* I.e. "carry in previous addition?" */
1536
                        break;
1537
                      n1 -= n0 < d1;
1538
                      n0 -= d1;
1539
                    }
1540
                }
1541
 
1542
              /* Possible optimization: We already have (q * n0) and (1 * n1)
1543
                 after the calculation of QUOT.  Taking advantage of this, we
1544
                 could make this loop make two iterations less.  */
1545
 
1546
              cy = mpn_submul_1 (num, den, densize + 1, quot);
1547
 
1548
              if (num[densize] != cy)
1549
                {
1550
                  cy = mpn_add_n (num, num, den, densize);
1551
                  assert (cy != 0);
1552
                  --quot;
1553
                }
1554
              n0 = num[densize] = num[densize - 1];
1555
              for (i = densize - 1; i > 0; --i)
1556
                num[i] = num[i - 1];
1557
 
1558
              got_limb;
1559
            }
1560
 
1561
          for (i = densize; num[i] == 0 && i >= 0; --i)
1562
            ;
1563
          return round_and_return (retval, exponent - 1, negative,
1564
                                   quot, BITS_PER_MP_LIMB - 1 - used,
1565
                                   more_bits || i >= 0);
1566
        }
1567
      }
1568
  }
1569
 
1570
  /* NOTREACHED */
1571
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.