1 /* Read decimal floating point numbers.
2 Copyright (C) 1995 Free Software Foundation, Inc.
3 Contributed by Ulrich Drepper.
5 This file is part of the GNU C Library.
7 The GNU C Library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Library General Public License as
9 published by the Free Software Foundation; either version 2 of the
10 License, or (at your option) any later version.
12 The GNU C Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 Library General Public License for more details.
17 You should have received a copy of the GNU Library General Public
18 License along with the GNU C Library; see the file COPYING.LIB. If
19 not, write to the Free Software Foundation, Inc., 675 Mass Ave,
20 Cambridge, MA 02139, USA. */
22 /* Configuration part. These macros are defined by `strtold.c' and `strtof.c'
23 to produce the `long double' and `float' versions of the reader. */
28 #define MPN2FLOAT __mpn_construct_double
29 #define FLOAT_HUGE_VAL HUGE_VAL
31 /* End of configuration part. */
36 #include "../locale/localeinfo.h"
41 #include <gmp-mparam.h>
43 #include "fpioconst.h"
49 /* Constants we need from float.h; select the set for the FLOAT precision. */
50 #define MANT_DIG PASTE(FLT,_MANT_DIG)
51 #define DIG PASTE(FLT,_DIG)
52 #define MAX_EXP PASTE(FLT,_MAX_EXP)
53 #define MIN_EXP PASTE(FLT,_MIN_EXP)
54 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
55 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
57 /* Extra macros required to get FLT expanded before the pasting. */
58 #define PASTE(a,b) PASTE1(a,b)
59 #define PASTE1(a,b) a##b
61 /* Function to construct a floating point number from an MP integer
62 containing the fraction bits, a base 2 exponent, and a sign flag. */
63 extern FLOAT
MPN2FLOAT (mp_srcptr mpn
, int exponent
, int negative
);
65 /* Definitions according to limb size used. */
66 #if BITS_PER_MP_LIMB == 32
67 # define MAX_DIG_PER_LIMB 9
68 # define MAX_FAC_PER_LIMB 1000000000L
69 #elif BITS_PER_MP_LIMB == 64
70 # define MAX_DIG_PER_LIMB 19
71 # define MAX_FAC_PER_LIMB 10000000000000000000L
73 # error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
77 /* Local data structure. */
78 static const mp_limb _tens_in_limb
[MAX_DIG_PER_LIMB
+ 1] =
81 1000000, 10000000, 100000000,
83 #if BITS_PER_MP_LIMB > 32
84 , 10000000000, 100000000000,
85 1000000000000, 10000000000000, 100000000000000,
86 1000000000000000, 10000000000000000, 100000000000000000,
87 1000000000000000000, 10000000000000000000
89 #if BITS_PER_MP_LIMB > 64
90 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
95 #define howmany(x,y) (((x)+((y)-1))/(y))
97 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
99 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
100 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
102 #define RETURN(val,end) \
103 do { if (endptr != 0) *endptr = (char *) (end); return val; } while (0)
105 /* Maximum size necessary for mpn integers to hold floating point numbers. */
106 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
108 /* Declare an mpn integer variable that big. */
109 #define MPN_VAR(name) mp_limb name[MPNSIZE]; mp_size_t name##size
110 /* Copy an mpn integer value. */
111 #define MPN_ASSIGN(dst, src) \
112 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb))
115 /* Return a floating point number of the needed type according to the given
116 multi-precision number after possible rounding. */
118 round_and_return (mp_limb
*retval
, int exponent
, int negative
,
119 mp_limb round_limb
, mp_size_t round_bit
, int more_bits
)
121 if (exponent
< MIN_EXP
- 1)
123 mp_size_t shift
= MIN_EXP
- 1 - exponent
;
125 if (shift
> MANT_DIG
)
131 more_bits
|= (round_limb
& ((1 << round_bit
) - 1)) != 0;
132 if (shift
== MANT_DIG
)
133 /* This is a special case to handle the very seldom case where
134 the mantissa will be empty after the shift. */
138 round_limb
= retval
[RETURN_LIMB_SIZE
- 1];
139 round_bit
= BITS_PER_MP_LIMB
- 1;
140 for (i
= 0; i
< RETURN_LIMB_SIZE
; ++i
)
141 more_bits
|= retval
[i
] != 0;
142 MPN_ZERO (retval
, RETURN_LIMB_SIZE
);
144 else if (shift
>= BITS_PER_MP_LIMB
)
148 round_limb
= retval
[(shift
- 1) / BITS_PER_MP_LIMB
];
149 round_bit
= (shift
- 1) % BITS_PER_MP_LIMB
;
150 for (i
= 0; i
< (shift
- 1) / BITS_PER_MP_LIMB
; ++i
)
151 more_bits
|= retval
[i
] != 0;
152 more_bits
|= (round_limb
& ((1 << round_bit
) - 1)) != 0;
154 (void) __mpn_rshift (retval
, &retval
[shift
/ BITS_PER_MP_LIMB
],
155 RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
),
156 shift
% BITS_PER_MP_LIMB
);
157 MPN_ZERO (&retval
[RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
)],
158 shift
/ BITS_PER_MP_LIMB
);
162 round_limb
= retval
[0];
163 round_bit
= shift
- 1;
164 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
166 exponent
= MIN_EXP
- 2;
169 if ((round_limb
& (1 << round_bit
)) != 0
170 && (more_bits
|| (retval
[0] & 1) != 0
171 || (round_limb
& ((1 << round_bit
) - 1)) != 0))
173 mp_limb cy
= __mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
175 if (((MANT_DIG
% BITS_PER_MP_LIMB
) == 0 && cy
) ||
176 ((MANT_DIG
% BITS_PER_MP_LIMB
) != 0 &&
177 (retval
[RETURN_LIMB_SIZE
- 1]
178 & (1 << (MANT_DIG
% BITS_PER_MP_LIMB
))) != 0))
181 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, 1);
182 retval
[RETURN_LIMB_SIZE
- 1] |= 1 << ((MANT_DIG
- 1)
185 else if (exponent
== MIN_EXP
- 2
186 && (retval
[RETURN_LIMB_SIZE
- 1]
187 & (1 << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
))) != 0)
188 /* The number was denormalized but now normalized. */
189 exponent
= MIN_EXP
- 1;
192 if (exponent
> MAX_EXP
)
193 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
195 return MPN2FLOAT (retval
, exponent
, negative
);
199 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
200 into N. Return the size of the number limbs in NSIZE at the first
201 character od the string that is not part of the integer as the function
202 value. If the EXPONENT is small enough to be taken as an additional
203 factor for the resulting number (see code) multiply by it. */
204 static inline const char *
205 str_to_mpn (const char *str
, int digcnt
, mp_limb
*n
, mp_size_t
*nsize
,
208 /* Number of digits for actual limb. */
217 if (cnt
== MAX_DIG_PER_LIMB
)
224 cy
= __mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
225 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
234 /* There might be thousands separators or radix characters in the string.
235 But these all can be ignored because we know the format of the number
236 is correct and we have an exact number of characters to read. */
237 while (!isdigit (*str
))
239 low
= low
* 10 + *str
++ - '0';
242 while (--digcnt
> 0);
244 if (*exponent
> 0 && cnt
+ *exponent
<= MAX_DIG_PER_LIMB
)
246 low
*= _tens_in_limb
[*exponent
];
247 base
= _tens_in_limb
[cnt
+ *exponent
];
251 base
= _tens_in_limb
[cnt
];
261 cy
= __mpn_mul_1 (n
, n
, *nsize
, base
);
262 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
270 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
271 with the COUNT most significant bits of LIMB.
273 Tege doesn't like this function so I have to write it here myself. :)
276 __mpn_lshift_1 (mp_limb
*ptr
, mp_size_t size
, unsigned int count
, mp_limb limb
)
278 if (count
== BITS_PER_MP_LIMB
)
280 /* Optimize the case of shifting by exactly a word:
281 just copy words, with no actual bit-shifting. */
283 for (i
= size
- 1; i
> 0; --i
)
289 (void) __mpn_lshift (ptr
, ptr
, size
, count
);
290 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
295 #define INTERNAL(x) INTERNAL1(x)
296 #define INTERNAL1(x) __##x##_internal
298 /* This file defines a function to check for correct grouping. */
299 #include "grouping.h"
302 /* Return a floating point number with the value of the given string NPTR.
303 Set *ENDPTR to the character after the last used one. If the number is
304 smaller than the smallest representable number, set `errno' to ERANGE and
305 return 0.0. If the number is too big to be represented, set `errno' to
306 ERANGE and return HUGE_VAL with the approriate sign. */
308 INTERNAL (STRTOF
) (nptr
, endptr
, group
)
313 int negative
; /* The sign of the number. */
314 MPN_VAR (num
); /* MP representation of the number. */
315 int exponent
; /* Exponent of the number. */
317 /* When we have to compute fractional digits we form a fraction with a
318 second multi-precision number (and we sometimes need a second for
319 temporary results). */
322 /* Representation for the return value. */
323 mp_limb retval
[RETURN_LIMB_SIZE
];
324 /* Number of bits currently in result value. */
327 /* Running pointer after the last character processed in the string. */
329 /* Start of significant part of the number. */
330 const char *startp
, *start_of_digits
;
331 /* Points at the character following the integer and fractional digits. */
333 /* Total number of digit and number of digits in integer part. */
334 int dig_no
, int_no
, lead_zero
;
335 /* Contains the last character read. */
338 /* The radix character of the current locale. */
340 /* The thousands character of the current locale. */
342 /* The numeric grouping specification of the current locale,
343 in the format described in <locale.h>. */
344 const char *grouping
;
348 grouping
= _NL_CURRENT (LC_NUMERIC
, GROUPING
);
349 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
353 /* Figure out the thousands separator character. */
354 if (mbtowc (&thousands
, _NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
),
355 strlen (_NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
))) <= 0)
356 thousands
= (wchar_t) *_NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
);
357 if (thousands
== L
'\0')
367 /* Find the locale's decimal point character. */
368 if (mbtowc (&decimal
, _NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
),
369 strlen (_NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
))) <= 0)
370 decimal
= (wchar_t) *_NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
);
373 /* Prepare number representation. */
378 /* Parse string to get maximal legal prefix. We need the number of
379 characters of the integer part, the fractional part and the exponent. */
381 /* Ignore leading white space. */
386 /* Get sign of the result. */
395 /* Return 0.0 if no legal string is found.
396 No character is used even if a sign was found. */
397 if (!isdigit (c
) && (c
!= decimal
|| !isdigit (cp
[1])))
400 /* Record the start of the digits, in case we will check their grouping. */
401 start_of_digits
= startp
= cp
;
403 /* Ignore leading zeroes. This helps us to avoid useless computations. */
404 while (c
== '0' || (thousands
!= L
'\0' && c
== thousands
))
407 /* If no other digit but a '0' is found the result is 0.0.
408 Return current read pointer. */
409 if (!isdigit (c
) && c
!= decimal
)
411 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
412 /* If TP is at the start of the digits, there was no correctly
413 grouped prefix of the string; so no number found. */
414 RETURN (0.0, tp
== start_of_digits
? nptr
: tp
);
417 /* Remember first significant digit and read following characters until the
418 decimal point, exponent character or any non-FP number character. */
421 while (dig_no
< NDIG
||
422 /* If parsing grouping info, keep going past useful digits
423 so we can check all the grouping separators. */
428 else if (thousands
== L
'\0' || c
!= thousands
)
429 /* Not a digit or separator: end of the integer part. */
434 if (grouping
&& dig_no
> 0)
436 /* Check the grouping of the digits. */
437 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
440 /* Less than the entire string was correctly grouped. */
442 if (tp
== start_of_digits
)
443 /* No valid group of numbers at all: no valid number. */
447 /* The number is validly grouped, but consists
448 only of zeroes. The whole value is zero. */
451 /* Recompute DIG_NO so we won't read more digits than
452 are properly grouped. */
455 for (tp
= startp
; tp
< cp
; ++tp
)
467 /* Too many digits to be representable. Assigning this to EXPONENT
468 allows us to read the full number but return HUGE_VAL after parsing. */
469 exponent
= MAX_10_EXP
;
471 /* We have the number digits in the integer part. Whether these are all or
472 any is really a fractional digit will be decided later. */
474 lead_zero
= int_no
== 0 ? -1 : 0;
476 /* Read the fractional digits. A special case are the 'american style'
477 numbers like `16.' i.e. with decimal but without trailing digits. */
485 if (c
!= '0' && lead_zero
== -1)
486 lead_zero
= dig_no
- int_no
;
494 /* Remember start of exponent (if any). */
498 if (tolower (c
) == 'e')
500 int exp_negative
= 0;
515 /* Get the exponent limit. */
516 exp_limit
= exp_negative
?
517 -MIN_10_EXP
+ MANT_DIG
- int_no
:
518 MAX_10_EXP
- int_no
+ lead_zero
;
524 if (exponent
> exp_limit
)
525 /* The exponent is too large/small to represent a valid
530 /* Overflow or underflow. */
532 retval
= (exp_negative
? 0.0 :
533 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
535 /* Accept all following digits as part of the exponent. */
538 while (isdigit (*cp
));
553 exponent
= -exponent
;
556 /* We don't want to have to work with trailing zeroes after the radix. */
559 while (expp
[-1] == '0')
564 assert (dig_no
>= int_no
);
569 /* The whole string is parsed. Store the address of the next character. */
571 *endptr
= (char *) cp
;
578 /* Find the decimal point */
579 while (*startp
!= decimal
) startp
++;
580 startp
+= lead_zero
+ 1;
581 exponent
-= lead_zero
;
585 /* Now we have the number of digits in total and the integer digits as well
586 as the exponent and its sign. We can decide whether the read digits are
587 really integer digits or belong to the fractional part; i.e. we normalize
590 register int incr
= exponent
< 0 ? MAX (-int_no
, exponent
)
591 : MIN (dig_no
- int_no
, exponent
);
596 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
599 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
602 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
610 /* Read the integer part as a multi-precision number to NUM. */
611 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
);
615 /* We now multiply the gained number by the given power of ten. */
617 mp_limb
*pdest
= den
;
619 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
623 if ((exponent
& expbit
) != 0)
628 /* FIXME: not the whole multiplication has to be done.
629 If we have the needed number of bits we only need the
630 information whether more non-zero bits follow. */
631 if (numsize
>= ttab
->arraysize
- 2)
632 cy
= __mpn_mul (pdest
, psrc
, numsize
,
633 &ttab
->array
[2], ttab
->arraysize
- 2);
635 cy
= __mpn_mul (pdest
, &ttab
->array
[2],
638 numsize
+= ttab
->arraysize
- 2;
646 while (exponent
!= 0);
649 memcpy (num
, den
, numsize
* sizeof (mp_limb
));
652 /* Determine how many bits of the result we already have. */
653 count_leading_zeros (bits
, num
[numsize
- 1]);
654 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
656 /* Now we know the exponent of the number in base two.
657 Check it against the maximum possible exponent. */
661 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
664 /* We have already the first BITS bits of the result. Together with
665 the information whether more non-zero bits follow this is enough
666 to determine the result. */
670 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
671 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
672 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
674 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
678 memcpy (retval
, &num
[least_idx
],
679 RETURN_LIMB_SIZE
* sizeof (mp_limb
));
682 for (i
= least_idx
; i
< numsize
- 1; ++i
)
683 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
685 << (BITS_PER_MP_LIMB
- least_bit
));
686 if (i
- least_idx
< RETURN_LIMB_SIZE
)
687 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
690 /* Check whether any limb beside the ones in RETVAL are non-zero. */
691 for (i
= 0; num
[i
] == 0; ++i
)
694 return round_and_return (retval
, bits
- 1, negative
,
695 num
[round_idx
], round_bit
,
696 int_no
< dig_no
|| i
< round_idx
);
699 else if (dig_no
== int_no
)
701 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
702 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
704 if (target_bit
== is_bit
)
706 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
707 numsize
* sizeof (mp_limb
));
708 /* FIXME: the following loop can be avoided if we assume a
709 maximal MANT_DIG value. */
710 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
712 else if (target_bit
> is_bit
)
714 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
715 num
, numsize
, target_bit
- is_bit
);
716 /* FIXME: the following loop can be avoided if we assume a
717 maximal MANT_DIG value. */
718 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
723 assert (numsize
< RETURN_LIMB_SIZE
);
725 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
726 num
, numsize
, is_bit
- target_bit
);
727 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
728 /* FIXME: the following loop can be avoided if we assume a
729 maximal MANT_DIG value. */
730 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
733 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
737 /* Store the bits we already have. */
738 memcpy (retval
, num
, numsize
* sizeof (mp_limb
));
739 #if RETURN_LIMB_SIZE > 1
740 if (numsize
< RETURN_LIMB_SIZE
)
745 /* We have to compute at least some of the fractional digits. */
747 /* We construct a fraction and the result of the division gives us
748 the needed digits. The denominator is 1.0 multiplied by the
749 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
750 123e-6 gives 123 / 1000000. */
758 mp_limb
*pdest
= num
;
759 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
761 assert (dig_no
> int_no
&& exponent
<= 0);
764 /* For the fractional part we need not process too much digits. One
765 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
767 digits we should have enough bits for the result. The remaining
768 decimal digits give us the information that more bits are following.
769 This can be used while rounding. (One added as a safety margin.) */
770 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 1)
772 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 1;
778 neg_exp
= dig_no
- int_no
- exponent
;
780 /* Construct the denominator. */
785 if ((neg_exp
& expbit
) != 0)
791 memcpy (psrc
, &ttab
->array
[2],
792 (densize
= ttab
->arraysize
- 2) * sizeof (mp_limb
));
795 cy
= __mpn_mul (pdest
, &ttab
->array
[2], ttab
->arraysize
- 2,
797 densize
+= ttab
->arraysize
- 2;
806 while (neg_exp
!= 0);
809 memcpy (den
, num
, densize
* sizeof (mp_limb
));
811 /* Read the fractional digits from the string. */
812 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
);
815 /* We now have to shift both numbers so that the highest bit in the
816 denominator is set. In the same process we copy the numerator to
817 a high place in the array so that the division constructs the wanted
818 digits. This is done by a "quasi fix point" number representation.
820 num: ddddddddddd . 0000000000000000000000
822 den: ddddddddddd n >= m
826 count_leading_zeros (cnt
, den
[densize
- 1]);
828 (void) __mpn_lshift (den
, den
, densize
, cnt
);
829 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
833 /* Now we are ready for the division. But it is not necessary to
834 do a full multi-precision division because we only need a small
835 number of bits for the result. So we do not use __mpn_divmod
836 here but instead do the division here by hand and stop whenever
837 the needed number of bits is reached. The code itself comes
838 from the GNU MP Library by Torbj\"orn Granlund. */
851 assert (numsize
== 1 && n
< d
);
855 udiv_qrnnd (quot
, n
, n
, 0, d
);
862 cnt = BITS_PER_MP_LIMB; \
864 count_leading_zeros (cnt, quot); \
866 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
868 used = MANT_DIG + cnt; \
869 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
870 bits = MANT_DIG + 1; \
874 /* Note that we only clear the second element. */ \
875 /* The conditional is determined at compile time. */ \
876 if (RETURN_LIMB_SIZE > 1) \
882 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
883 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
887 used = MANT_DIG - bits; \
889 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
891 bits += BITS_PER_MP_LIMB
895 while (bits
<= MANT_DIG
);
897 return round_and_return (retval
, exponent
- 1, negative
,
898 quot
, BITS_PER_MP_LIMB
- 1 - used
,
899 more_bits
|| n
!= 0);
903 mp_limb d0
, d1
, n0
, n1
;
910 if (numsize
< densize
)
914 /* The numerator of the number occupies fewer bits than
915 the denominator but the one limb is bigger than the
916 high limb of the numerator. */
923 exponent
-= BITS_PER_MP_LIMB
;
926 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
927 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
928 BITS_PER_MP_LIMB
, 0);
931 used
= MANT_DIG
- bits
;
933 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
935 bits
+= BITS_PER_MP_LIMB
;
947 while (bits
<= MANT_DIG
)
953 /* QUOT should be either 111..111 or 111..110. We need
954 special treatment of this rare case as normal division
955 would give overflow. */
959 if (r
< d1
) /* Carry in the addition? */
961 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
969 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
970 umul_ppmm (n1
, n0
, d0
, quot
);
974 if (n1
> r
|| (n1
== r
&& n0
> 0))
976 /* The estimated QUOT was too large. */
979 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
981 if (r
>= d1
) /* If not carry, test QUOT again. */
984 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
990 return round_and_return (retval
, exponent
- 1, negative
,
991 quot
, BITS_PER_MP_LIMB
- 1 - used
,
992 more_bits
|| n1
!= 0 || n0
!= 0);
997 mp_limb cy
, dX
, d1
, n0
, n1
;
1001 dX
= den
[densize
- 1];
1002 d1
= den
[densize
- 2];
1004 /* The division does not work if the upper limb of the two-limb
1005 numerator is greater than the denominator. */
1006 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1009 if (numsize
< densize
)
1011 mp_size_t empty
= densize
- numsize
;
1016 for (i
= numsize
; i
> 0; --i
)
1017 num
[i
+ empty
] = num
[i
- 1];
1018 MPN_ZERO (num
, empty
+ 1);
1019 exponent
-= empty
* BITS_PER_MP_LIMB
;
1023 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1025 /* We make a difference here because the compiler
1026 cannot optimize the `else' case that good and
1027 this reflects all currently used FLOAT types
1028 and GMP implementations. */
1030 #if RETURN_LIMB_SIZE <= 2
1031 assert (empty
== 1);
1032 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1033 BITS_PER_MP_LIMB
, 0);
1035 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1036 retval
[i
] = retval
[i
- empty
];
1039 for (i
= numsize
; i
> 0; --i
)
1040 num
[i
+ empty
] = num
[i
- 1];
1041 MPN_ZERO (num
, empty
+ 1);
1045 used
= MANT_DIG
- bits
;
1046 if (used
>= BITS_PER_MP_LIMB
)
1049 (void) __mpn_lshift (&retval
[used
1050 / BITS_PER_MP_LIMB
],
1051 retval
, RETURN_LIMB_SIZE
,
1052 used
% BITS_PER_MP_LIMB
);
1053 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1057 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1059 bits
+= empty
* BITS_PER_MP_LIMB
;
1065 assert (numsize
== densize
);
1066 for (i
= numsize
; i
> 0; --i
)
1067 num
[i
] = num
[i
- 1];
1073 while (bits
<= MANT_DIG
)
1076 /* This might over-estimate QUOT, but it's probably not
1077 worth the extra code here to find out. */
1078 quot
= ~(mp_limb
) 0;
1083 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1084 umul_ppmm (n1
, n0
, d1
, quot
);
1086 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1090 if (r
< dX
) /* I.e. "carry in previous addition?" */
1097 /* Possible optimization: We already have (q * n0) and (1 * n1)
1098 after the calculation of QUOT. Taking advantage of this, we
1099 could make this loop make two iterations less. */
1101 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1103 if (num
[densize
] != cy
)
1105 cy
= __mpn_add_n (num
, num
, den
, densize
);
1109 n0
= num
[densize
] = num
[densize
- 1];
1110 for (i
= densize
- 1; i
> 0; --i
)
1111 num
[i
] = num
[i
- 1];
1116 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1118 return round_and_return (retval
, exponent
- 1, negative
,
1119 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1120 more_bits
|| i
>= 0);
1128 /* External user entry point. */
1131 STRTOF (nptr
, endptr
)
1135 return INTERNAL (STRTOF
) (nptr
, endptr
, 0);
1138 #define weak_this(x) weak_symbol(x)