1 /* Read decimal floating point numbers.
2 Copyright (C) 1995, 1996 Free Software Foundation, Inc.
3 Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
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., 59 Temple Place - Suite 330,
20 Boston, MA 02111-1307, USA. */
22 /* Configuration part. These macros are defined by `strtold.c',
23 `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the
24 `long double' and `float' versions of the reader. */
29 # define STRTOF wcstod
31 # define STRTOF strtod
33 # define MPN2FLOAT __mpn_construct_double
34 # define FLOAT_HUGE_VAL HUGE_VAL
40 # define STRING_TYPE wchar_t
41 # define CHAR_TYPE wint_t
43 # define ISSPACE(Ch) iswspace (Ch)
44 # define TOLOWER(Ch) towlower (Ch)
46 # define STRING_TYPE char
47 # define CHAR_TYPE char
49 # define ISSPACE(Ch) isspace (Ch)
50 # define TOLOWER(Ch) tolower (Ch)
52 /* End of configuration part. */
57 #include "../locale/localeinfo.h"
62 /* The gmp headers need some configuration frobs. */
67 #include <gmp-mparam.h>
69 #include "fpioconst.h"
75 /* Constants we need from float.h; select the set for the FLOAT precision. */
76 #define MANT_DIG PASTE(FLT,_MANT_DIG)
77 #define DIG PASTE(FLT,_DIG)
78 #define MAX_EXP PASTE(FLT,_MAX_EXP)
79 #define MIN_EXP PASTE(FLT,_MIN_EXP)
80 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
81 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
83 /* Extra macros required to get FLT expanded before the pasting. */
84 #define PASTE(a,b) PASTE1(a,b)
85 #define PASTE1(a,b) a##b
87 /* Function to construct a floating point number from an MP integer
88 containing the fraction bits, a base 2 exponent, and a sign flag. */
89 extern FLOAT
MPN2FLOAT (mp_srcptr mpn
, int exponent
, int negative
);
91 /* Definitions according to limb size used. */
92 #if BITS_PER_MP_LIMB == 32
93 # define MAX_DIG_PER_LIMB 9
94 # define MAX_FAC_PER_LIMB 1000000000UL
95 #elif BITS_PER_MP_LIMB == 64
96 # define MAX_DIG_PER_LIMB 19
97 # define MAX_FAC_PER_LIMB 10000000000000000000UL
99 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
103 /* Local data structure. */
104 static const mp_limb_t _tens_in_limb
[MAX_DIG_PER_LIMB
+ 1] =
107 1000000, 10000000, 100000000,
109 #if BITS_PER_MP_LIMB > 32
110 , 10000000000, 100000000000,
111 1000000000000, 10000000000000, 100000000000000,
112 1000000000000000, 10000000000000000, 100000000000000000,
113 1000000000000000000, 10000000000000000000U
115 #if BITS_PER_MP_LIMB > 64
116 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
121 #define howmany(x,y) (((x)+((y)-1))/(y))
123 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
125 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
126 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
128 #define RETURN(val,end) \
129 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
130 return val; } while (0)
132 /* Maximum size necessary for mpn integers to hold floating point numbers. */
133 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
135 /* Declare an mpn integer variable that big. */
136 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
137 /* Copy an mpn integer value. */
138 #define MPN_ASSIGN(dst, src) \
139 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
142 /* Return a floating point number of the needed type according to the given
143 multi-precision number after possible rounding. */
145 round_and_return (mp_limb_t
*retval
, int exponent
, int negative
,
146 mp_limb_t round_limb
, mp_size_t round_bit
, int more_bits
)
148 if (exponent
< MIN_EXP
- 1)
150 mp_size_t shift
= MIN_EXP
- 1 - exponent
;
152 if (shift
> MANT_DIG
)
158 more_bits
|= (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0;
159 if (shift
== MANT_DIG
)
160 /* This is a special case to handle the very seldom case where
161 the mantissa will be empty after the shift. */
165 round_limb
= retval
[RETURN_LIMB_SIZE
- 1];
166 round_bit
= BITS_PER_MP_LIMB
- 1;
167 for (i
= 0; i
< RETURN_LIMB_SIZE
; ++i
)
168 more_bits
|= retval
[i
] != 0;
169 MPN_ZERO (retval
, RETURN_LIMB_SIZE
);
171 else if (shift
>= BITS_PER_MP_LIMB
)
175 round_limb
= retval
[(shift
- 1) / BITS_PER_MP_LIMB
];
176 round_bit
= (shift
- 1) % BITS_PER_MP_LIMB
;
177 for (i
= 0; i
< (shift
- 1) / BITS_PER_MP_LIMB
; ++i
)
178 more_bits
|= retval
[i
] != 0;
179 more_bits
|= ((round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1))
182 (void) __mpn_rshift (retval
, &retval
[shift
/ BITS_PER_MP_LIMB
],
183 RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
),
184 shift
% BITS_PER_MP_LIMB
);
185 MPN_ZERO (&retval
[RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
)],
186 shift
/ BITS_PER_MP_LIMB
);
190 round_limb
= retval
[0];
191 round_bit
= shift
- 1;
192 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
194 exponent
= MIN_EXP
- 2;
197 if ((round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0
198 && (more_bits
|| (retval
[0] & 1) != 0
199 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0))
201 mp_limb_t cy
= __mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
203 if (((MANT_DIG
% BITS_PER_MP_LIMB
) == 0 && cy
) ||
204 ((MANT_DIG
% BITS_PER_MP_LIMB
) != 0 &&
205 (retval
[RETURN_LIMB_SIZE
- 1]
206 & (((mp_limb_t
) 1) << (MANT_DIG
% BITS_PER_MP_LIMB
))) != 0))
209 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, 1);
210 retval
[RETURN_LIMB_SIZE
- 1]
211 |= ((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
);
213 else if (exponent
== MIN_EXP
- 2
214 && (retval
[RETURN_LIMB_SIZE
- 1]
215 & (((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
)))
217 /* The number was denormalized but now normalized. */
218 exponent
= MIN_EXP
- 1;
221 if (exponent
> MAX_EXP
)
222 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
224 return MPN2FLOAT (retval
, exponent
, negative
);
228 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
229 into N. Return the size of the number limbs in NSIZE at the first
230 character od the string that is not part of the integer as the function
231 value. If the EXPONENT is small enough to be taken as an additional
232 factor for the resulting number (see code) multiply by it. */
233 static inline const STRING_TYPE
*
234 str_to_mpn (const STRING_TYPE
*str
, int digcnt
, mp_limb_t
*n
, mp_size_t
*nsize
,
237 /* Number of digits for actual limb. */
246 if (cnt
== MAX_DIG_PER_LIMB
)
253 cy
= __mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
254 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
263 /* There might be thousands separators or radix characters in the string.
264 But these all can be ignored because we know the format of the number
265 is correct and we have an exact number of characters to read. */
266 while (*str
< L_('0') || *str
> L_('9'))
268 low
= low
* 10 + *str
++ - L_('0');
271 while (--digcnt
> 0);
273 if (*exponent
> 0 && cnt
+ *exponent
<= MAX_DIG_PER_LIMB
)
275 low
*= _tens_in_limb
[*exponent
];
276 base
= _tens_in_limb
[cnt
+ *exponent
];
280 base
= _tens_in_limb
[cnt
];
290 cy
= __mpn_mul_1 (n
, n
, *nsize
, base
);
291 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
299 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
300 with the COUNT most significant bits of LIMB.
302 Tege doesn't like this function so I have to write it here myself. :)
305 __mpn_lshift_1 (mp_limb_t
*ptr
, mp_size_t size
, unsigned int count
,
308 if (count
== BITS_PER_MP_LIMB
)
310 /* Optimize the case of shifting by exactly a word:
311 just copy words, with no actual bit-shifting. */
313 for (i
= size
- 1; i
> 0; --i
)
319 (void) __mpn_lshift (ptr
, ptr
, size
, count
);
320 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
325 #define INTERNAL(x) INTERNAL1(x)
326 #define INTERNAL1(x) __##x##_internal
328 /* This file defines a function to check for correct grouping. */
329 #include "grouping.h"
332 /* Return a floating point number with the value of the given string NPTR.
333 Set *ENDPTR to the character after the last used one. If the number is
334 smaller than the smallest representable number, set `errno' to ERANGE and
335 return 0.0. If the number is too big to be represented, set `errno' to
336 ERANGE and return HUGE_VAL with the approriate sign. */
338 INTERNAL (STRTOF
) (nptr
, endptr
, group
)
339 const STRING_TYPE
*nptr
;
340 STRING_TYPE
**endptr
;
343 int negative
; /* The sign of the number. */
344 MPN_VAR (num
); /* MP representation of the number. */
345 int exponent
; /* Exponent of the number. */
347 /* When we have to compute fractional digits we form a fraction with a
348 second multi-precision number (and we sometimes need a second for
349 temporary results). */
352 /* Representation for the return value. */
353 mp_limb_t retval
[RETURN_LIMB_SIZE
];
354 /* Number of bits currently in result value. */
357 /* Running pointer after the last character processed in the string. */
358 const STRING_TYPE
*cp
, *tp
;
359 /* Start of significant part of the number. */
360 const STRING_TYPE
*startp
, *start_of_digits
;
361 /* Points at the character following the integer and fractional digits. */
362 const STRING_TYPE
*expp
;
363 /* Total number of digit and number of digits in integer part. */
364 int dig_no
, int_no
, lead_zero
;
365 /* Contains the last character read. */
368 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
369 there. So define it ourselves if it remains undefined. */
371 typedef unsigned int wint_t;
373 /* The radix character of the current locale. */
375 /* The thousands character of the current locale. */
377 /* The numeric grouping specification of the current locale,
378 in the format described in <locale.h>. */
379 const char *grouping
;
381 assert (sizeof (wchar_t) == sizeof (wint_t));
385 grouping
= _NL_CURRENT (LC_NUMERIC
, GROUPING
);
386 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
390 /* Figure out the thousands separator character. */
391 if (mbtowc ((wchar_t *) &thousands
,
392 _NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
),
393 strlen (_NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
))) <= 0)
394 thousands
= (wint_t) *_NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
);
395 if (thousands
== L
'\0')
405 /* Find the locale's decimal point character. */
406 if (mbtowc ((wchar_t *) &decimal
, _NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
),
407 strlen (_NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
))) <= 0)
408 decimal
= (wint_t) *_NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
);
411 /* Prepare number representation. */
416 /* Parse string to get maximal legal prefix. We need the number of
417 characters of the integer part, the fractional part and the exponent. */
419 /* Ignore leading white space. */
424 /* Get sign of the result. */
430 else if (c
== L_('+'))
433 /* Return 0.0 if no legal string is found.
434 No character is used even if a sign was found. */
435 if ((c
< L_('0') || c
> L_('9'))
436 && ((wint_t) c
!= decimal
|| cp
[1] < L_('0') || cp
[1] > L_('9')))
439 /* Record the start of the digits, in case we will check their grouping. */
440 start_of_digits
= startp
= cp
;
442 /* Ignore leading zeroes. This helps us to avoid useless computations. */
443 while (c
== L_('0') || (thousands
!= L
'\0' && (wint_t) c
== thousands
))
446 /* If no other digit but a '0' is found the result is 0.0.
447 Return current read pointer. */
448 if ((c
< L_('0') || c
> L_('9')) && (wint_t) c
!= decimal
)
450 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
451 /* If TP is at the start of the digits, there was no correctly
452 grouped prefix of the string; so no number found. */
453 RETURN (0.0, tp
== start_of_digits
? nptr
: tp
);
456 /* Remember first significant digit and read following characters until the
457 decimal point, exponent character or any non-FP number character. */
460 while (dig_no
< NDIG
||
461 /* If parsing grouping info, keep going past useful digits
462 so we can check all the grouping separators. */
465 if (c
>= L_('0') && c
<= L_('9'))
467 else if (thousands
== L
'\0' || (wint_t) c
!= thousands
)
468 /* Not a digit or separator: end of the integer part. */
473 if (grouping
&& dig_no
> 0)
475 /* Check the grouping of the digits. */
476 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
479 /* Less than the entire string was correctly grouped. */
481 if (tp
== start_of_digits
)
482 /* No valid group of numbers at all: no valid number. */
486 /* The number is validly grouped, but consists
487 only of zeroes. The whole value is zero. */
490 /* Recompute DIG_NO so we won't read more digits than
491 are properly grouped. */
494 for (tp
= startp
; tp
< cp
; ++tp
)
495 if (*tp
>= L_('0') && *tp
<= L_('9'))
506 /* Too many digits to be representable. Assigning this to EXPONENT
507 allows us to read the full number but return HUGE_VAL after parsing. */
508 exponent
= MAX_10_EXP
;
510 /* We have the number digits in the integer part. Whether these are all or
511 any is really a fractional digit will be decided later. */
513 lead_zero
= int_no
== 0 ? -1 : 0;
515 /* Read the fractional digits. A special case are the 'american style'
516 numbers like `16.' i.e. with decimal but without trailing digits. */
517 if ((wint_t) c
== decimal
)
520 while (c
>= L_('0') && c
<= L_('9'))
522 if (c
!= L_('0') && lead_zero
== -1)
523 lead_zero
= dig_no
- int_no
;
529 /* Remember start of exponent (if any). */
533 if (TOLOWER (c
) == L_('e'))
535 int exp_negative
= 0;
543 else if (c
== L_('+'))
546 if (c
>= L_('0') && c
<= L_('9'))
550 /* Get the exponent limit. */
551 exp_limit
= exp_negative
?
552 -MIN_10_EXP
+ MANT_DIG
- int_no
:
553 MAX_10_EXP
- int_no
+ lead_zero
;
559 if (exponent
> exp_limit
)
560 /* The exponent is too large/small to represent a valid
565 /* Overflow or underflow. */
566 __set_errno (ERANGE
);
567 retval
= (exp_negative
? 0.0 :
568 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
570 /* Accept all following digits as part of the exponent. */
573 while (*cp
>= L_('0') && *cp
<= L_('9'));
579 exponent
+= c
- L_('0');
582 while (c
>= L_('0') && c
<= L_('9'));
585 exponent
= -exponent
;
591 /* We don't want to have to work with trailing zeroes after the radix. */
594 while (expp
[-1] == L_('0'))
599 assert (dig_no
>= int_no
);
604 /* The whole string is parsed. Store the address of the next character. */
606 *endptr
= (STRING_TYPE
*) cp
;
613 /* Find the decimal point */
614 while ((wint_t) *startp
!= decimal
)
616 startp
+= lead_zero
+ 1;
617 exponent
-= lead_zero
;
621 /* Now we have the number of digits in total and the integer digits as well
622 as the exponent and its sign. We can decide whether the read digits are
623 really integer digits or belong to the fractional part; i.e. we normalize
626 register int incr
= exponent
< 0 ? MAX (-int_no
, exponent
)
627 : MIN (dig_no
- int_no
, exponent
);
632 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
634 __set_errno (ERANGE
);
635 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
638 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
640 __set_errno (ERANGE
);
646 /* Read the integer part as a multi-precision number to NUM. */
647 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
);
651 /* We now multiply the gained number by the given power of ten. */
652 mp_limb_t
*psrc
= num
;
653 mp_limb_t
*pdest
= den
;
655 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
659 if ((exponent
& expbit
) != 0)
664 /* FIXME: not the whole multiplication has to be done.
665 If we have the needed number of bits we only need the
666 information whether more non-zero bits follow. */
667 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
668 cy
= __mpn_mul (pdest
, psrc
, numsize
,
669 &ttab
->array
[_FPIO_CONST_OFFSET
],
670 ttab
->arraysize
- _FPIO_CONST_OFFSET
);
672 cy
= __mpn_mul (pdest
, &ttab
->array
[_FPIO_CONST_OFFSET
],
673 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
675 numsize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
683 while (exponent
!= 0);
686 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
689 /* Determine how many bits of the result we already have. */
690 count_leading_zeros (bits
, num
[numsize
- 1]);
691 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
693 /* Now we know the exponent of the number in base two.
694 Check it against the maximum possible exponent. */
697 __set_errno (ERANGE
);
698 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
701 /* We have already the first BITS bits of the result. Together with
702 the information whether more non-zero bits follow this is enough
703 to determine the result. */
707 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
708 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
709 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
711 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
715 memcpy (retval
, &num
[least_idx
],
716 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
719 for (i
= least_idx
; i
< numsize
- 1; ++i
)
720 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
722 << (BITS_PER_MP_LIMB
- least_bit
));
723 if (i
- least_idx
< RETURN_LIMB_SIZE
)
724 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
727 /* Check whether any limb beside the ones in RETVAL are non-zero. */
728 for (i
= 0; num
[i
] == 0; ++i
)
731 return round_and_return (retval
, bits
- 1, negative
,
732 num
[round_idx
], round_bit
,
733 int_no
< dig_no
|| i
< round_idx
);
736 else if (dig_no
== int_no
)
738 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
739 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
741 if (target_bit
== is_bit
)
743 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
744 numsize
* sizeof (mp_limb_t
));
745 /* FIXME: the following loop can be avoided if we assume a
746 maximal MANT_DIG value. */
747 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
749 else if (target_bit
> is_bit
)
751 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
752 num
, numsize
, target_bit
- is_bit
);
753 /* FIXME: the following loop can be avoided if we assume a
754 maximal MANT_DIG value. */
755 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
760 assert (numsize
< RETURN_LIMB_SIZE
);
762 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
763 num
, numsize
, is_bit
- target_bit
);
764 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
765 /* FIXME: the following loop can be avoided if we assume a
766 maximal MANT_DIG value. */
767 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
770 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
774 /* Store the bits we already have. */
775 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
776 #if RETURN_LIMB_SIZE > 1
777 if (numsize
< RETURN_LIMB_SIZE
)
782 /* We have to compute at least some of the fractional digits. */
784 /* We construct a fraction and the result of the division gives us
785 the needed digits. The denominator is 1.0 multiplied by the
786 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
787 123e-6 gives 123 / 1000000. */
794 mp_limb_t
*psrc
= den
;
795 mp_limb_t
*pdest
= num
;
796 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
798 assert (dig_no
> int_no
&& exponent
<= 0);
801 /* For the fractional part we need not process too much digits. One
802 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
804 digits we should have enough bits for the result. The remaining
805 decimal digits give us the information that more bits are following.
806 This can be used while rounding. (One added as a safety margin.) */
807 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 1)
809 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 1;
815 neg_exp
= dig_no
- int_no
- exponent
;
817 /* Construct the denominator. */
822 if ((neg_exp
& expbit
) != 0)
829 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
830 memcpy (psrc
, &ttab
->array
[_FPIO_CONST_OFFSET
],
831 densize
* sizeof (mp_limb_t
));
835 cy
= __mpn_mul (pdest
, &ttab
->array
[_FPIO_CONST_OFFSET
],
836 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
838 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
847 while (neg_exp
!= 0);
850 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
852 /* Read the fractional digits from the string. */
853 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
);
856 /* We now have to shift both numbers so that the highest bit in the
857 denominator is set. In the same process we copy the numerator to
858 a high place in the array so that the division constructs the wanted
859 digits. This is done by a "quasi fix point" number representation.
861 num: ddddddddddd . 0000000000000000000000
863 den: ddddddddddd n >= m
867 count_leading_zeros (cnt
, den
[densize
- 1]);
869 (void) __mpn_lshift (den
, den
, densize
, cnt
);
870 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
874 /* Now we are ready for the division. But it is not necessary to
875 do a full multi-precision division because we only need a small
876 number of bits for the result. So we do not use __mpn_divmod
877 here but instead do the division here by hand and stop whenever
878 the needed number of bits is reached. The code itself comes
879 from the GNU MP Library by Torbj\"orn Granlund. */
887 mp_limb_t d
, n
, quot
;
892 assert (numsize
== 1 && n
< d
);
896 udiv_qrnnd (quot
, n
, n
, 0, d
);
903 cnt = BITS_PER_MP_LIMB; \
905 count_leading_zeros (cnt, quot); \
907 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
909 used = MANT_DIG + cnt; \
910 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
911 bits = MANT_DIG + 1; \
915 /* Note that we only clear the second element. */ \
916 /* The conditional is determined at compile time. */ \
917 if (RETURN_LIMB_SIZE > 1) \
923 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
924 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
928 used = MANT_DIG - bits; \
930 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
932 bits += BITS_PER_MP_LIMB
936 while (bits
<= MANT_DIG
);
938 return round_and_return (retval
, exponent
- 1, negative
,
939 quot
, BITS_PER_MP_LIMB
- 1 - used
,
940 more_bits
|| n
!= 0);
944 mp_limb_t d0
, d1
, n0
, n1
;
951 if (numsize
< densize
)
955 /* The numerator of the number occupies fewer bits than
956 the denominator but the one limb is bigger than the
957 high limb of the numerator. */
964 exponent
-= BITS_PER_MP_LIMB
;
967 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
968 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
969 BITS_PER_MP_LIMB
, 0);
972 used
= MANT_DIG
- bits
;
974 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
976 bits
+= BITS_PER_MP_LIMB
;
988 while (bits
<= MANT_DIG
)
994 /* QUOT should be either 111..111 or 111..110. We need
995 special treatment of this rare case as normal division
996 would give overflow. */
997 quot
= ~(mp_limb_t
) 0;
1000 if (r
< d1
) /* Carry in the addition? */
1002 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1005 n1
= d0
- (d0
!= 0);
1010 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1011 umul_ppmm (n1
, n0
, d0
, quot
);
1015 if (n1
> r
|| (n1
== r
&& n0
> 0))
1017 /* The estimated QUOT was too large. */
1020 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1022 if (r
>= d1
) /* If not carry, test QUOT again. */
1025 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1031 return round_and_return (retval
, exponent
- 1, negative
,
1032 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1033 more_bits
|| n1
!= 0 || n0
!= 0);
1038 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1042 dX
= den
[densize
- 1];
1043 d1
= den
[densize
- 2];
1045 /* The division does not work if the upper limb of the two-limb
1046 numerator is greater than the denominator. */
1047 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1050 if (numsize
< densize
)
1052 mp_size_t empty
= densize
- numsize
;
1057 for (i
= numsize
; i
> 0; --i
)
1058 num
[i
+ empty
] = num
[i
- 1];
1059 MPN_ZERO (num
, empty
+ 1);
1060 exponent
-= empty
* BITS_PER_MP_LIMB
;
1064 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1066 /* We make a difference here because the compiler
1067 cannot optimize the `else' case that good and
1068 this reflects all currently used FLOAT types
1069 and GMP implementations. */
1071 #if RETURN_LIMB_SIZE <= 2
1072 assert (empty
== 1);
1073 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1074 BITS_PER_MP_LIMB
, 0);
1076 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1077 retval
[i
] = retval
[i
- empty
];
1080 for (i
= numsize
; i
> 0; --i
)
1081 num
[i
+ empty
] = num
[i
- 1];
1082 MPN_ZERO (num
, empty
+ 1);
1086 used
= MANT_DIG
- bits
;
1087 if (used
>= BITS_PER_MP_LIMB
)
1090 (void) __mpn_lshift (&retval
[used
1091 / BITS_PER_MP_LIMB
],
1092 retval
, RETURN_LIMB_SIZE
,
1093 used
% BITS_PER_MP_LIMB
);
1094 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1098 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1100 bits
+= empty
* BITS_PER_MP_LIMB
;
1106 assert (numsize
== densize
);
1107 for (i
= numsize
; i
> 0; --i
)
1108 num
[i
] = num
[i
- 1];
1114 while (bits
<= MANT_DIG
)
1117 /* This might over-estimate QUOT, but it's probably not
1118 worth the extra code here to find out. */
1119 quot
= ~(mp_limb_t
) 0;
1124 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1125 umul_ppmm (n1
, n0
, d1
, quot
);
1127 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1131 if (r
< dX
) /* I.e. "carry in previous addition?" */
1138 /* Possible optimization: We already have (q * n0) and (1 * n1)
1139 after the calculation of QUOT. Taking advantage of this, we
1140 could make this loop make two iterations less. */
1142 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1144 if (num
[densize
] != cy
)
1146 cy
= __mpn_add_n (num
, num
, den
, densize
);
1150 n0
= num
[densize
] = num
[densize
- 1];
1151 for (i
= densize
- 1; i
> 0; --i
)
1152 num
[i
] = num
[i
- 1];
1157 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1159 return round_and_return (retval
, exponent
- 1, negative
,
1160 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1161 more_bits
|| i
>= 0);
1169 /* External user entry point. */
1172 #ifdef weak_function
1175 STRTOF (nptr
, endptr
)
1176 const STRING_TYPE
*nptr
;
1177 STRING_TYPE
**endptr
;
1179 return INTERNAL (STRTOF
) (nptr
, endptr
, 0);