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 <localeinfo.h>
39 #include "../stdio/gmp.h"
40 #include "../stdio/gmp-impl.h"
41 #include <gmp-mparam.h>
42 #include "../stdio/longlong.h"
43 #include "../stdio/fpioconst.h"
45 /* #define NDEBUG 1 */
49 /* Constants we need from float.h; select the set for the FLOAT precision. */
50 #define MANT_DIG FLT##_MANT_DIG
51 #define MAX_EXP FLT##_MAX_EXP
52 #define MIN_EXP FLT##_MIN_EXP
53 #define MAX_10_EXP FLT##_MAX_10_EXP
54 #define MIN_10_EXP FLT##_MIN_10_EXP
55 #define MAX_10_EXP_LOG FLT##_MAX_10_EXP_LOG
58 /* Function to construct a floating point number from an MP integer
59 containing the fraction bits, a base 2 exponent, and a sign flag. */
60 extern FLOAT
MPN2FLOAT (mp_srcptr mpn
, int exponent
, int negative
);
62 /* Definitions according to limb size used. */
63 #if BITS_PER_MP_LIMB == 32
64 # define MAX_DIG_PER_LIMB 9
65 # define MAX_FAC_PER_LIMB 1000000000L
66 #elif BITS_PER_MP_LIMB == 64
67 # define MAX_DIG_PER_LIMB 19
68 # define MAX_FAC_PER_LIMB 10000000000000000000L
70 # error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
74 /* Local data structure. */
75 static const mp_limb _tens_in_limb
[MAX_DIG_PER_LIMB
] =
78 1000000, 10000000, 100000000
79 #if BITS_PER_MP_LIMB > 32
80 , 1000000000, 10000000000, 100000000000,
81 1000000000000, 10000000000000, 100000000000000,
82 1000000000000000, 10000000000000000, 100000000000000000,
85 #if BITS_PER_MP_LIMB > 64
86 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
91 #define howmany(x,y) (((x)+((y)-1))/(y))
93 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
95 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
96 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
98 #define RETURN(val,end) \
99 do { if (endptr != 0) *endptr = (char *) end; return val; } while (0)
101 /* Maximum size necessary for mpn integers to hold floating point numbers. */
102 #define MPNSIZE (howmany (MAX_EXP + MANT_DIG, BITS_PER_MP_LIMB) + 1)
103 /* Declare an mpn integer variable that big. */
104 #define MPN_VAR(name) mp_limb name[MPNSIZE]; mp_size_t name##size
105 /* Copy an mpn integer value. */
106 #define MPN_ASSIGN(dst, src) \
107 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb))
110 /* Return a floating point number of the needed type according to the given
111 multi-precision number after possible rounding. */
113 round_and_return (mp_limb
*retval
, int exponent
, int negative
,
114 mp_limb round_limb
, mp_size_t round_bit
, int more_bits
)
116 if (exponent
< MIN_EXP
)
118 mp_size_t shift
= MIN_EXP
- 1 - exponent
;
120 if (shift
>= MANT_DIG
)
126 more_bits
|= (round_limb
& ((1 << round_bit
) - 1)) != 0;
127 if (shift
>= BITS_PER_MP_LIMB
)
129 round_limb
= retval
[(shift
- 1) / BITS_PER_MP_LIMB
];
130 round_bit
= (shift
- 1) % BITS_PER_MP_LIMB
;
131 #if RETURN_LIMB_SIZE <= 2
132 assert (RETURN_LIMB_SIZE
== 2);
133 more_bits
|= retval
[0] != 0;
134 retval
[0] = retval
[1];
137 int disp
= shift
/ BITS_PER_MP_LIMB
;
139 while (retval
[i
] == 0 && i
< disp
)
141 more_bits
|= i
< disp
;
142 for (i
= disp
; i
< RETURN_LIMB_SIZE
; ++i
)
143 retval
[i
- disp
] = retval
[i
];
144 MPN_ZERO (&retval
[RETURN_LIMB_SIZE
- disp
], disp
);
146 shift
%= BITS_PER_MP_LIMB
;
150 round_limb
= retval
[0];
151 round_bit
= shift
- 1;
153 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
154 exponent
= MIN_EXP
- 2;
157 if ((round_limb
& (1 << round_bit
)) != 0 &&
158 (more_bits
|| (retval
[0] & 1) != 0 ||
159 (round_limb
& ((1 << round_bit
) - 1)) != 0))
161 mp_limb cy
= __mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
162 if (cy
|| (retval
[RETURN_LIMB_SIZE
- 1]
163 & (1 << (MANT_DIG
% BITS_PER_MP_LIMB
))) != 0)
166 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, 1);
167 retval
[RETURN_LIMB_SIZE
- 1] |= 1 << (MANT_DIG
% BITS_PER_MP_LIMB
);
171 if (exponent
> MAX_EXP
)
172 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
174 return MPN2FLOAT (retval
, exponent
, negative
);
178 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
179 into N. Return the size of the number limbs in NSIZE at the first
180 character od the string that is not part of the integer as the function
181 value. If the EXPONENT is small enough to be taken as an additional
182 factor for the resulting number (see code) multiply by it. */
183 static inline const char *
184 str_to_mpn (const char *str
, int digcnt
, mp_limb
*n
, mp_size_t
*nsize
,
187 /* Number of digits for actual limb. */
196 if (cnt
== MAX_DIG_PER_LIMB
)
203 cy
= __mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
204 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
213 /* There might be thousands separators or radix characters in the string.
214 But these all can be ignored because we know the format of the number
215 is correct and we have an exact number of characters to read. */
216 while (!isdigit (*str
))
218 low
= low
* 10 + *str
++ - '0';
221 while (--digcnt
> 0);
223 if (*exponent
> 0 && cnt
+ *exponent
<= MAX_DIG_PER_LIMB
)
225 low
*= _tens_in_limb
[*exponent
];
226 base
= _tens_in_limb
[cnt
+ *exponent
];
230 base
= _tens_in_limb
[cnt
];
240 cy
= __mpn_mul_1 (n
, n
, *nsize
, base
);
241 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
249 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
250 with the COUNT most significant bits of LIMB.
252 Tege doesn't like this function so I have to write it here myself. :)
255 __mpn_lshift_1 (mp_limb
*ptr
, mp_size_t size
, unsigned int count
, mp_limb limb
)
257 if (count
== BITS_PER_MP_LIMB
)
259 /* Optimize the case of shifting by exactly a word:
260 just copy words, with no actual bit-shifting. */
262 for (i
= size
- 1; i
> 0; --i
)
268 (void) __mpn_lshift (ptr
, ptr
, size
, count
);
269 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
274 /* Return a floating point number with the value of the given string NPTR.
275 Set *ENDPTR to the character after the last used one. If the number is
276 smaller than the smallest representable number, set `errno' to ERANGE and
277 return 0.0. If the number is too big to be represented, set `errno' to
278 ERANGE and return HUGE_VAL with the approriate sign. */
280 STRTOF (nptr
, endptr
)
284 int negative
; /* The sign of the number. */
285 MPN_VAR (num
); /* MP representation of the number. */
286 int exponent
; /* Exponent of the number. */
288 /* When we have to compute fractional digits we form a fraction with a
289 second multi-precision number (and we sometimes need a second for
290 temporary results). */
293 /* Representation for the return value. */
294 mp_limb retval
[RETURN_LIMB_SIZE
];
295 /* Number of bits currently in result value. */
298 /* Running pointer after the last character processed in the string. */
300 /* Start of significant part of the number. */
302 /* Points at the character following the integer and fractional digits. */
304 /* Total number of digit and number of digits in integer part. */
306 /* Contains the last character read. */
309 /* The radix character of the current locale. */
312 /* The thousands character of the current locale. */
314 /* The numeric grouping specification of the current locale,
315 in the format described in <locale.h>. */
316 const char *grouping
;
318 /* Check the grouping of the integer part at [BEGIN,END).
319 Return zero iff a separator is found out of place. */
320 int grouping_ok (const char *begin
, const char *end
)
328 while (*p
!= thousands
&& p
> begin
);
329 if (end
- 1 - p
!= *grouping
++)
330 return 0; /* Wrong number of digits in this group. */
331 end
= p
; /* Correct group; trim it off the end. */
334 --grouping
; /* Same grouping repeats in next iteration. */
335 else if (*grouping
== CHAR_MAX
|| *grouping
< 0)
337 /* No further grouping allowed. */
339 if (*--end
== thousands
)
345 /* Return with no conversion if the grouping of [STARTP,CP) is bad. */
346 #define CHECK_GROUPING if (! grouping_ok (startp, cp)) RETURN (0.0, nptr); else
348 grouping
= _numeric_info
->grouping
; /* Cache the grouping info array. */
349 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
353 /* Figure out the thousands seperator character. */
354 if (mbtowc (&thousands_sep
, _numeric_info
->thousands_sep
,
355 strlen (_numeric_info
->thousands_sep
)) <= 0)
356 thousands
= (wchar_t) *_numeric_info
->thousands_sep
;
357 if (thousands
== L
'\0')
361 #define grouping NULL
362 #define thousands L'\0'
363 #define CHECK_GROUPING ((void) 0)
366 /* Find the locale's decimal point character. */
367 if (mbtowc (&decimal
, _numeric_info
->decimal_point
,
368 strlen (_numeric_info
->decimal_point
)) <= 0)
369 decimal
= (wchar_t) *_numeric_info
->decimal_point
;
372 /* Prepare number representation. */
377 /* Parse string to get maximal legal prefix. We need the number of
378 characters of the interger part, the fractional part and the exponent. */
380 /* Ignore leading white space. */
385 /* Get sign of the result. */
394 /* Return 0.0 if no legal string is found.
395 No character is used even if a sign was found. */
399 /* Record the start of the digits, in case we will check their grouping. */
402 /* Ignore leading zeroes. This helps us to avoid useless computations. */
403 while (c
== '0' || (thousands
!= L
'\0' && c
== thousands
))
408 /* If no other digit but a '0' is found the result is 0.0.
409 Return current read pointer. */
410 if (!isdigit (c
) && c
!= decimal
)
413 /* Remember first significant digit and read following characters until the
414 decimal point, exponent character or any non-FP number character. */
417 while (dig_no
< NDIG
||
418 /* If parsing grouping info, keep going past useful digits
419 so we can check all the grouping separators. */
424 else if (thousands
== L
'\0' || c
!= thousands
)
425 /* Not a digit or separator: end of the integer part. */
433 /* Too many digits to be representable. Assigning this to EXPONENT
434 allows us to read the full number but return HUGE_VAL after parsing. */
435 exponent
= MAX_10_EXP
;
437 /* We have the number digits in the integer part. Whether these are all or
438 any is really a fractional digit will be decided later. */
441 /* Read the fractional digits. */
456 /* Remember start of exponent (if any). */
460 if (tolower (c
) == 'e')
462 int exp_negative
= 0;
477 if ((!exp_negative
&& exponent
* 10 + int_no
> MAX_10_EXP
)
479 && exponent
* 10 + int_no
> -MIN_10_EXP
+ MANT_DIG
))
480 /* The exponent is too large/small to represent a valid
485 /* Overflow or underflow. */
487 retval
= (exp_negative
? 0.0 :
488 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
490 /* Accept all following digits as part of the exponent. */
493 while (isdigit (*cp
));
509 exponent
= -exponent
;
512 /* We don't want to have to work with trailing zeroes after the radix. */
515 while (expp
[-1] == '0')
520 assert (dig_no
>= int_no
);
523 /* The whole string is parsed. Store the address of the next character. */
525 *endptr
= (char *) cp
;
530 /* Now we have the number of digits in total and the integer digits as well
531 as the exponent and its sign. We can decide whether the read digits are
532 really integer digits or belong to the fractional part; i.e. we normalize
535 register int incr
= exponent
< 0 ? MAX (-int_no
, exponent
)
536 : MIN (dig_no
- int_no
, exponent
);
541 if (int_no
+ exponent
> MAX_10_EXP
)
544 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
547 if (int_no
- dig_no
+ exponent
< MIN_10_EXP
- MANT_DIG
)
555 /* Read the integer part as a multi-precision number to NUM. */
556 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
);
560 /* We now multiply the gained number by the given power of ten. */
562 mp_limb
*pdest
= den
;
564 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
566 assert (exponent
< (1 << (MAX_10_EXP_LOG
+ 1)));
569 if ((exponent
& expbit
) != 0)
574 /* FIXME: not the whole multiplication has to be done.
575 If we have the needed number of bits we only need the
576 information whether more non-zero bits follow. */
577 if (numsize
>= ttab
->arraysize
- 2)
578 cy
= __mpn_mul (pdest
, psrc
, numsize
,
579 &ttab
->array
[2], ttab
->arraysize
- 2);
581 cy
= __mpn_mul (pdest
, &ttab
->array
[2],
584 numsize
+= ttab
->arraysize
- 2;
592 while (exponent
!= 0);
595 memcpy (num
, den
, numsize
* sizeof (mp_limb
));
598 /* Determine how many bits of the result we already have. */
599 count_leading_zeros (bits
, num
[numsize
- 1]);
600 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
602 /* We have already the first BITS bits of the result. Together with
603 the information whether more non-zero bits follow this is enough
604 to determine the result. */
607 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
608 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
609 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
611 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
616 memcpy (retval
, &num
[least_idx
],
617 RETURN_LIMB_SIZE
* sizeof (mp_limb
));
619 (void) __mpn_rshift (retval
, &num
[least_idx
],
620 numsize
- least_idx
+ 1, least_bit
);
622 /* Check whether any limb beside the ones in RETVAL are non-zero. */
623 for (i
= 0; num
[i
] == 0; ++i
)
626 return round_and_return (retval
, bits
- 1, negative
,
627 num
[round_idx
], round_bit
,
628 int_no
< dig_no
|| i
< round_idx
);
631 else if (dig_no
== int_no
)
633 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
634 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
636 if (target_bit
== is_bit
)
638 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
639 numsize
* sizeof (mp_limb
));
640 /* FIXME: the following loop can be avoided if we assume a
641 maximal MANT_DIG value. */
642 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
644 else if (target_bit
> is_bit
)
646 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
647 num
, numsize
, target_bit
- is_bit
);
648 /* FIXME: the following loop can be avoided if we assume a
649 maximal MANT_DIG value. */
650 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
655 assert (numsize
< RETURN_LIMB_SIZE
);
657 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
658 num
, numsize
, is_bit
- target_bit
);
659 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
660 /* FIXME: the following loop can be avoided if we assume a
661 maximal MANT_DIG value. */
662 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
665 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
669 /* Store the bits we already have. */
670 memcpy (retval
, num
, numsize
* sizeof (mp_limb
));
671 #if RETURN_LIMB_SIZE > 1
672 if (numsize
< RETURN_LIMB_SIZE
)
677 /* We have to compute at least some of the fractional digits. */
679 /* We construct a fraction and the result of the division gives us
680 the needed digits. The denominator is 1.0 multiplied by the
681 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
682 123e6 gives 123 / 1000000. */
688 mp_limb
*pdest
= num
;
689 int neg_exp
= dig_no
- int_no
- exponent
;
690 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
692 assert (dig_no
> int_no
&& exponent
<= 0);
694 /* Construct the denominator. */
699 if ((neg_exp
& expbit
) != 0)
705 memcpy (psrc
, &ttab
->array
[2],
706 (densize
= ttab
->arraysize
- 2) * sizeof (mp_limb
));
709 cy
= __mpn_mul (pdest
, &ttab
->array
[2], ttab
->arraysize
- 2,
711 densize
+= ttab
->arraysize
- 2;
720 while (neg_exp
!= 0);
723 memcpy (den
, num
, densize
* sizeof (mp_limb
));
725 /* Read the fractional digits from the string. */
726 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
);
729 /* We now have to shift both numbers so that the highest bit in the
730 denominator is set. In the same process we copy the numerator to
731 a high place in the array so that the division constructs the wanted
732 digits. This is done by a "quasi fix point" number representation.
734 num: ddddddddddd . 0000000000000000000000
736 den: ddddddddddd n >= m
740 count_leading_zeros (cnt
, den
[densize
- 1]);
742 (void) __mpn_lshift (den
, den
, densize
, cnt
);
743 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
747 /* Now we are ready for the division. But it is not necessary to
748 do a full multi-precision division because we only need a small
749 number of bits for the result. So we do not use __mpn_divmod
750 here but instead do the division here by hand and stop whenever
751 the needed number of bits is reached. The code itself comes
752 from the GNU MP Library by Torbj\"orn Granlund. */
765 assert (numsize
== 1 && n
< d
);
769 udiv_qrnnd (quot
, n
, n
, 0, d
);
776 cnt = BITS_PER_MP_LIMB; \
778 count_leading_zeros (cnt, quot); \
780 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
782 used = cnt + MANT_DIG; \
783 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
784 bits -= BITS_PER_MP_LIMB - used; \
788 /* Note that we only clear the second element. */ \
794 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
795 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
799 used = MANT_DIG - bits; \
801 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
803 bits += BITS_PER_MP_LIMB
807 while (bits
<= MANT_DIG
);
809 return round_and_return (retval
, exponent
- 1, negative
,
810 quot
, BITS_PER_MP_LIMB
- 1 - used
,
815 mp_limb d0
, d1
, n0
, n1
;
822 if (numsize
< densize
)
825 exponent
-= BITS_PER_MP_LIMB
;
828 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
829 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
830 BITS_PER_MP_LIMB
, 0);
833 used
= MANT_DIG
- bits
;
835 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
837 bits
+= BITS_PER_MP_LIMB
;
848 while (bits
<= MANT_DIG
)
854 /* QUOT should be either 111..111 or 111..110. We need
855 special treatment of this rare case as normal division
856 would give overflow. */
860 if (r
< d1
) /* Carry in the addition? */
862 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
870 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
871 umul_ppmm (n1
, n0
, d0
, quot
);
875 if (n1
> r
|| (n1
== r
&& n0
> 0))
877 /* The estimated QUOT was too large. */
880 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
882 if (r
>= d1
) /* If not carry, test QUOT again. */
885 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
891 return round_and_return (retval
, exponent
- 1, negative
,
892 quot
, BITS_PER_MP_LIMB
- 1 - used
,
898 mp_limb cy
, dX
, d1
, n0
, n1
;
902 dX
= den
[densize
- 1];
903 d1
= den
[densize
- 2];
905 /* The division does not work if the upper limb of the two-limb
906 numerator is greater than the denominator. */
907 if (num
[numsize
- 1] > dX
)
910 if (numsize
< densize
)
912 mp_size_t empty
= densize
- numsize
;
917 for (i
= numsize
; i
> 0; --i
)
918 num
[i
+ empty
] = num
[i
- 1];
919 MPN_ZERO (num
, empty
+ 1);
920 exponent
-= empty
* BITS_PER_MP_LIMB
;
924 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
926 /* We make a difference here because the compiler
927 cannot optimize the `else' case that good and
928 this reflects all currently used FLOAT types
929 and GMP implementations. */
931 #if RETURN_LIMB_SIZE <= 2
933 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
934 BITS_PER_MP_LIMB
, 0);
936 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
937 retval
[i
] = retval
[i
- empty
];
940 for (i
= numsize
; i
> 0; --i
)
941 num
[i
+ empty
] = num
[i
- 1];
942 MPN_ZERO (num
, empty
+ 1);
946 used
= MANT_DIG
- bits
;
947 if (used
>= BITS_PER_MP_LIMB
)
950 (void) __mpn_lshift (&retval
[used
952 retval
, RETURN_LIMB_SIZE
,
953 used
% BITS_PER_MP_LIMB
);
954 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
958 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
960 bits
+= empty
* BITS_PER_MP_LIMB
;
966 assert (numsize
== densize
);
967 for (i
= numsize
; i
> 0; --i
)
974 while (bits
<= MANT_DIG
)
977 /* This might over-estimate QUOT, but it's probably not
978 worth the extra code here to find out. */
984 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
985 umul_ppmm (n1
, n0
, d1
, quot
);
987 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
991 if (r
< dX
) /* I.e. "carry in previous addition?" */
998 /* Possible optimization: We already have (q * n0) and (1 * n1)
999 after the calculation of QUOT. Taking advantage of this, we
1000 could make this loop make two iterations less. */
1002 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1004 if (num
[densize
] != cy
)
1006 cy
= __mpn_add_n (num
, num
, den
, densize
);
1010 n0
= num
[densize
] = num
[densize
- 1];
1011 for (i
= densize
- 1; i
> 0; --i
)
1012 num
[i
] = num
[i
- 1];
1017 for (i
= densize
- 1; num
[i
] != 0 && i
>= 0; --i
)
1019 return round_and_return (retval
, exponent
- 1, negative
,
1020 quot
, BITS_PER_MP_LIMB
- 1 - used
,