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
449 && TOLOWER (c
) != L_('e'))
451 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
452 /* If TP is at the start of the digits, there was no correctly
453 grouped prefix of the string; so no number found. */
454 RETURN (0.0, tp
== start_of_digits
? nptr
: tp
);
457 /* Remember first significant digit and read following characters until the
458 decimal point, exponent character or any non-FP number character. */
461 while (dig_no
< NDIG
||
462 /* If parsing grouping info, keep going past useful digits
463 so we can check all the grouping separators. */
466 if (c
>= L_('0') && c
<= L_('9'))
468 else if (thousands
== L
'\0' || (wint_t) c
!= thousands
)
469 /* Not a digit or separator: end of the integer part. */
474 if (grouping
&& dig_no
> 0)
476 /* Check the grouping of the digits. */
477 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
480 /* Less than the entire string was correctly grouped. */
482 if (tp
== start_of_digits
)
483 /* No valid group of numbers at all: no valid number. */
487 /* The number is validly grouped, but consists
488 only of zeroes. The whole value is zero. */
491 /* Recompute DIG_NO so we won't read more digits than
492 are properly grouped. */
495 for (tp
= startp
; tp
< cp
; ++tp
)
496 if (*tp
>= L_('0') && *tp
<= L_('9'))
507 /* Too many digits to be representable. Assigning this to EXPONENT
508 allows us to read the full number but return HUGE_VAL after parsing. */
509 exponent
= MAX_10_EXP
;
511 /* We have the number digits in the integer part. Whether these are all or
512 any is really a fractional digit will be decided later. */
514 lead_zero
= int_no
== 0 ? -1 : 0;
516 /* Read the fractional digits. A special case are the 'american style'
517 numbers like `16.' i.e. with decimal but without trailing digits. */
518 if ((wint_t) c
== decimal
)
521 while (c
>= L_('0') && c
<= L_('9'))
523 if (c
!= L_('0') && lead_zero
== -1)
524 lead_zero
= dig_no
- int_no
;
530 /* Remember start of exponent (if any). */
534 if (TOLOWER (c
) == L_('e'))
536 int exp_negative
= 0;
544 else if (c
== L_('+'))
547 if (c
>= L_('0') && c
<= L_('9'))
551 /* Get the exponent limit. */
552 exp_limit
= exp_negative
?
553 -MIN_10_EXP
+ MANT_DIG
- int_no
:
554 MAX_10_EXP
- int_no
+ lead_zero
;
560 if (exponent
> exp_limit
)
561 /* The exponent is too large/small to represent a valid
566 /* Overflow or underflow. */
567 __set_errno (ERANGE
);
568 retval
= (exp_negative
? 0.0 :
569 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
571 /* Accept all following digits as part of the exponent. */
574 while (*cp
>= L_('0') && *cp
<= L_('9'));
580 exponent
+= c
- L_('0');
583 while (c
>= L_('0') && c
<= L_('9'));
586 exponent
= -exponent
;
592 /* We don't want to have to work with trailing zeroes after the radix. */
595 while (expp
[-1] == L_('0'))
600 assert (dig_no
>= int_no
);
605 /* The whole string is parsed. Store the address of the next character. */
607 *endptr
= (STRING_TYPE
*) cp
;
614 /* Find the decimal point */
615 while ((wint_t) *startp
!= decimal
)
617 startp
+= lead_zero
+ 1;
618 exponent
-= lead_zero
;
622 /* Now we have the number of digits in total and the integer digits as well
623 as the exponent and its sign. We can decide whether the read digits are
624 really integer digits or belong to the fractional part; i.e. we normalize
627 register int incr
= exponent
< 0 ? MAX (-int_no
, exponent
)
628 : MIN (dig_no
- int_no
, exponent
);
633 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
635 __set_errno (ERANGE
);
636 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
639 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
641 __set_errno (ERANGE
);
647 /* Read the integer part as a multi-precision number to NUM. */
648 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
);
652 /* We now multiply the gained number by the given power of ten. */
653 mp_limb_t
*psrc
= num
;
654 mp_limb_t
*pdest
= den
;
656 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
660 if ((exponent
& expbit
) != 0)
665 /* FIXME: not the whole multiplication has to be done.
666 If we have the needed number of bits we only need the
667 information whether more non-zero bits follow. */
668 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
669 cy
= __mpn_mul (pdest
, psrc
, numsize
,
670 &ttab
->array
[_FPIO_CONST_OFFSET
],
671 ttab
->arraysize
- _FPIO_CONST_OFFSET
);
673 cy
= __mpn_mul (pdest
, &ttab
->array
[_FPIO_CONST_OFFSET
],
674 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
676 numsize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
684 while (exponent
!= 0);
687 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
690 /* Determine how many bits of the result we already have. */
691 count_leading_zeros (bits
, num
[numsize
- 1]);
692 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
694 /* Now we know the exponent of the number in base two.
695 Check it against the maximum possible exponent. */
698 __set_errno (ERANGE
);
699 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
702 /* We have already the first BITS bits of the result. Together with
703 the information whether more non-zero bits follow this is enough
704 to determine the result. */
708 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
709 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
710 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
712 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
716 memcpy (retval
, &num
[least_idx
],
717 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
720 for (i
= least_idx
; i
< numsize
- 1; ++i
)
721 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
723 << (BITS_PER_MP_LIMB
- least_bit
));
724 if (i
- least_idx
< RETURN_LIMB_SIZE
)
725 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
728 /* Check whether any limb beside the ones in RETVAL are non-zero. */
729 for (i
= 0; num
[i
] == 0; ++i
)
732 return round_and_return (retval
, bits
- 1, negative
,
733 num
[round_idx
], round_bit
,
734 int_no
< dig_no
|| i
< round_idx
);
737 else if (dig_no
== int_no
)
739 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
740 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
742 if (target_bit
== is_bit
)
744 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
745 numsize
* sizeof (mp_limb_t
));
746 /* FIXME: the following loop can be avoided if we assume a
747 maximal MANT_DIG value. */
748 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
750 else if (target_bit
> is_bit
)
752 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
753 num
, numsize
, target_bit
- is_bit
);
754 /* FIXME: the following loop can be avoided if we assume a
755 maximal MANT_DIG value. */
756 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
761 assert (numsize
< RETURN_LIMB_SIZE
);
763 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
764 num
, numsize
, is_bit
- target_bit
);
765 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
766 /* FIXME: the following loop can be avoided if we assume a
767 maximal MANT_DIG value. */
768 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
771 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
775 /* Store the bits we already have. */
776 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
777 #if RETURN_LIMB_SIZE > 1
778 if (numsize
< RETURN_LIMB_SIZE
)
783 /* We have to compute at least some of the fractional digits. */
785 /* We construct a fraction and the result of the division gives us
786 the needed digits. The denominator is 1.0 multiplied by the
787 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
788 123e-6 gives 123 / 1000000. */
795 mp_limb_t
*psrc
= den
;
796 mp_limb_t
*pdest
= num
;
797 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
799 assert (dig_no
> int_no
&& exponent
<= 0);
802 /* For the fractional part we need not process too much digits. One
803 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
805 digits we should have enough bits for the result. The remaining
806 decimal digits give us the information that more bits are following.
807 This can be used while rounding. (One added as a safety margin.) */
808 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 1)
810 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 1;
816 neg_exp
= dig_no
- int_no
- exponent
;
818 /* Construct the denominator. */
823 if ((neg_exp
& expbit
) != 0)
830 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
831 memcpy (psrc
, &ttab
->array
[_FPIO_CONST_OFFSET
],
832 densize
* sizeof (mp_limb_t
));
836 cy
= __mpn_mul (pdest
, &ttab
->array
[_FPIO_CONST_OFFSET
],
837 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
839 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
848 while (neg_exp
!= 0);
851 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
853 /* Read the fractional digits from the string. */
854 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
);
857 /* We now have to shift both numbers so that the highest bit in the
858 denominator is set. In the same process we copy the numerator to
859 a high place in the array so that the division constructs the wanted
860 digits. This is done by a "quasi fix point" number representation.
862 num: ddddddddddd . 0000000000000000000000
864 den: ddddddddddd n >= m
868 count_leading_zeros (cnt
, den
[densize
- 1]);
870 (void) __mpn_lshift (den
, den
, densize
, cnt
);
871 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
875 /* Now we are ready for the division. But it is not necessary to
876 do a full multi-precision division because we only need a small
877 number of bits for the result. So we do not use __mpn_divmod
878 here but instead do the division here by hand and stop whenever
879 the needed number of bits is reached. The code itself comes
880 from the GNU MP Library by Torbj\"orn Granlund. */
888 mp_limb_t d
, n
, quot
;
893 assert (numsize
== 1 && n
< d
);
897 udiv_qrnnd (quot
, n
, n
, 0, d
);
904 cnt = BITS_PER_MP_LIMB; \
906 count_leading_zeros (cnt, quot); \
908 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
910 used = MANT_DIG + cnt; \
911 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
912 bits = MANT_DIG + 1; \
916 /* Note that we only clear the second element. */ \
917 /* The conditional is determined at compile time. */ \
918 if (RETURN_LIMB_SIZE > 1) \
924 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
925 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
929 used = MANT_DIG - bits; \
931 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
933 bits += BITS_PER_MP_LIMB
937 while (bits
<= MANT_DIG
);
939 return round_and_return (retval
, exponent
- 1, negative
,
940 quot
, BITS_PER_MP_LIMB
- 1 - used
,
941 more_bits
|| n
!= 0);
945 mp_limb_t d0
, d1
, n0
, n1
;
952 if (numsize
< densize
)
956 /* The numerator of the number occupies fewer bits than
957 the denominator but the one limb is bigger than the
958 high limb of the numerator. */
965 exponent
-= BITS_PER_MP_LIMB
;
968 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
969 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
970 BITS_PER_MP_LIMB
, 0);
973 used
= MANT_DIG
- bits
;
975 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
977 bits
+= BITS_PER_MP_LIMB
;
989 while (bits
<= MANT_DIG
)
995 /* QUOT should be either 111..111 or 111..110. We need
996 special treatment of this rare case as normal division
997 would give overflow. */
998 quot
= ~(mp_limb_t
) 0;
1001 if (r
< d1
) /* Carry in the addition? */
1003 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1006 n1
= d0
- (d0
!= 0);
1011 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1012 umul_ppmm (n1
, n0
, d0
, quot
);
1016 if (n1
> r
|| (n1
== r
&& n0
> 0))
1018 /* The estimated QUOT was too large. */
1021 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1023 if (r
>= d1
) /* If not carry, test QUOT again. */
1026 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1032 return round_and_return (retval
, exponent
- 1, negative
,
1033 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1034 more_bits
|| n1
!= 0 || n0
!= 0);
1039 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1043 dX
= den
[densize
- 1];
1044 d1
= den
[densize
- 2];
1046 /* The division does not work if the upper limb of the two-limb
1047 numerator is greater than the denominator. */
1048 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1051 if (numsize
< densize
)
1053 mp_size_t empty
= densize
- numsize
;
1058 for (i
= numsize
; i
> 0; --i
)
1059 num
[i
+ empty
] = num
[i
- 1];
1060 MPN_ZERO (num
, empty
+ 1);
1061 exponent
-= empty
* BITS_PER_MP_LIMB
;
1065 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1067 /* We make a difference here because the compiler
1068 cannot optimize the `else' case that good and
1069 this reflects all currently used FLOAT types
1070 and GMP implementations. */
1072 #if RETURN_LIMB_SIZE <= 2
1073 assert (empty
== 1);
1074 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1075 BITS_PER_MP_LIMB
, 0);
1077 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1078 retval
[i
] = retval
[i
- empty
];
1081 for (i
= numsize
; i
> 0; --i
)
1082 num
[i
+ empty
] = num
[i
- 1];
1083 MPN_ZERO (num
, empty
+ 1);
1087 used
= MANT_DIG
- bits
;
1088 if (used
>= BITS_PER_MP_LIMB
)
1091 (void) __mpn_lshift (&retval
[used
1092 / BITS_PER_MP_LIMB
],
1093 retval
, RETURN_LIMB_SIZE
,
1094 used
% BITS_PER_MP_LIMB
);
1095 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1099 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1101 bits
+= empty
* BITS_PER_MP_LIMB
;
1107 assert (numsize
== densize
);
1108 for (i
= numsize
; i
> 0; --i
)
1109 num
[i
] = num
[i
- 1];
1115 while (bits
<= MANT_DIG
)
1118 /* This might over-estimate QUOT, but it's probably not
1119 worth the extra code here to find out. */
1120 quot
= ~(mp_limb_t
) 0;
1125 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1126 umul_ppmm (n1
, n0
, d1
, quot
);
1128 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1132 if (r
< dX
) /* I.e. "carry in previous addition?" */
1139 /* Possible optimization: We already have (q * n0) and (1 * n1)
1140 after the calculation of QUOT. Taking advantage of this, we
1141 could make this loop make two iterations less. */
1143 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1145 if (num
[densize
] != cy
)
1147 cy
= __mpn_add_n (num
, num
, den
, densize
);
1151 n0
= num
[densize
] = num
[densize
- 1];
1152 for (i
= densize
- 1; i
> 0; --i
)
1153 num
[i
] = num
[i
- 1];
1158 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1160 return round_and_return (retval
, exponent
- 1, negative
,
1161 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1162 more_bits
|| i
>= 0);
1170 /* External user entry point. */
1173 #ifdef weak_function
1176 STRTOF (nptr
, endptr
)
1177 const STRING_TYPE
*nptr
;
1178 STRING_TYPE
**endptr
;
1180 return INTERNAL (STRTOF
) (nptr
, endptr
, 0);