Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/strtod/strtod_l.c @ 68:561a7518be6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update gcc-4.6
author Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
date Sun, 21 Aug 2011 07:07:55 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
67:f6334be47118 68:561a7518be6b
1 /* 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 }