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
;
383 grouping
= _NL_CURRENT (LC_NUMERIC
, GROUPING
);
384 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
388 /* Figure out the thousands separator character. */
389 if (mbtowc ((wchar_t *) &thousands
,
390 _NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
),
391 strlen (_NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
))) <= 0)
392 thousands
= (wint_t) *_NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
);
393 if (thousands
== L
'\0')
403 /* Find the locale's decimal point character. */
404 if (mbtowc ((wchar_t *) &decimal
, _NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
),
405 strlen (_NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
))) <= 0)
406 decimal
= (wint_t) *_NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
);
409 /* Prepare number representation. */
414 /* Parse string to get maximal legal prefix. We need the number of
415 characters of the integer part, the fractional part and the exponent. */
417 /* Ignore leading white space. */
422 /* Get sign of the result. */
428 else if (c
== L_('+'))
431 /* Return 0.0 if no legal string is found.
432 No character is used even if a sign was found. */
433 if ((c
< L_('0') || c
> L_('9'))
434 && ((wint_t) c
!= decimal
|| cp
[1] < L_('0') || cp
[1] > L_('9')))
437 /* Record the start of the digits, in case we will check their grouping. */
438 start_of_digits
= startp
= cp
;
440 /* Ignore leading zeroes. This helps us to avoid useless computations. */
441 while (c
== L_('0') || (thousands
!= L
'\0' && (wint_t) c
== thousands
))
444 /* If no other digit but a '0' is found the result is 0.0.
445 Return current read pointer. */
446 if ((c
< L_('0') || c
> L_('9')) && (wint_t) c
!= decimal
)
448 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
449 /* If TP is at the start of the digits, there was no correctly
450 grouped prefix of the string; so no number found. */
451 RETURN (0.0, tp
== start_of_digits
? nptr
: tp
);
454 /* Remember first significant digit and read following characters until the
455 decimal point, exponent character or any non-FP number character. */
458 while (dig_no
< NDIG
||
459 /* If parsing grouping info, keep going past useful digits
460 so we can check all the grouping separators. */
463 if (c
>= L_('0') && c
<= L_('9'))
465 else if (thousands
== L
'\0' || (wint_t) c
!= thousands
)
466 /* Not a digit or separator: end of the integer part. */
471 if (grouping
&& dig_no
> 0)
473 /* Check the grouping of the digits. */
474 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
477 /* Less than the entire string was correctly grouped. */
479 if (tp
== start_of_digits
)
480 /* No valid group of numbers at all: no valid number. */
484 /* The number is validly grouped, but consists
485 only of zeroes. The whole value is zero. */
488 /* Recompute DIG_NO so we won't read more digits than
489 are properly grouped. */
492 for (tp
= startp
; tp
< cp
; ++tp
)
493 if (*tp
>= L_('0') && *tp
<= L_('9'))
504 /* Too many digits to be representable. Assigning this to EXPONENT
505 allows us to read the full number but return HUGE_VAL after parsing. */
506 exponent
= MAX_10_EXP
;
508 /* We have the number digits in the integer part. Whether these are all or
509 any is really a fractional digit will be decided later. */
511 lead_zero
= int_no
== 0 ? -1 : 0;
513 /* Read the fractional digits. A special case are the 'american style'
514 numbers like `16.' i.e. with decimal but without trailing digits. */
515 if ((wint_t) c
== decimal
)
518 while (c
>= L_('0') && c
<= L_('9'))
520 if (c
!= L_('0') && lead_zero
== -1)
521 lead_zero
= dig_no
- int_no
;
527 /* Remember start of exponent (if any). */
531 if (TOLOWER (c
) == L_('e'))
533 int exp_negative
= 0;
541 else if (c
== L_('+'))
544 if (c
>= L_('0') && c
<= L_('9'))
548 /* Get the exponent limit. */
549 exp_limit
= exp_negative
?
550 -MIN_10_EXP
+ MANT_DIG
- int_no
:
551 MAX_10_EXP
- int_no
+ lead_zero
;
557 if (exponent
> exp_limit
)
558 /* The exponent is too large/small to represent a valid
563 /* Overflow or underflow. */
565 retval
= (exp_negative
? 0.0 :
566 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
568 /* Accept all following digits as part of the exponent. */
571 while (*cp
>= L_('0') && *cp
<= L_('9'));
577 exponent
+= c
- L_('0');
580 while (c
>= L_('0') && c
<= L_('9'));
583 exponent
= -exponent
;
589 /* We don't want to have to work with trailing zeroes after the radix. */
592 while (expp
[-1] == L_('0'))
597 assert (dig_no
>= int_no
);
602 /* The whole string is parsed. Store the address of the next character. */
604 *endptr
= (STRING_TYPE
*) cp
;
611 /* Find the decimal point */
612 while ((wint_t) *startp
!= decimal
)
614 startp
+= lead_zero
+ 1;
615 exponent
-= lead_zero
;
619 /* Now we have the number of digits in total and the integer digits as well
620 as the exponent and its sign. We can decide whether the read digits are
621 really integer digits or belong to the fractional part; i.e. we normalize
624 register int incr
= exponent
< 0 ? MAX (-int_no
, exponent
)
625 : MIN (dig_no
- int_no
, exponent
);
630 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
633 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
636 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
644 /* Read the integer part as a multi-precision number to NUM. */
645 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
);
649 /* We now multiply the gained number by the given power of ten. */
650 mp_limb_t
*psrc
= num
;
651 mp_limb_t
*pdest
= den
;
653 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
657 if ((exponent
& expbit
) != 0)
662 /* FIXME: not the whole multiplication has to be done.
663 If we have the needed number of bits we only need the
664 information whether more non-zero bits follow. */
665 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
666 cy
= __mpn_mul (pdest
, psrc
, numsize
,
667 &ttab
->array
[_FPIO_CONST_OFFSET
],
668 ttab
->arraysize
- _FPIO_CONST_OFFSET
);
670 cy
= __mpn_mul (pdest
, &ttab
->array
[_FPIO_CONST_OFFSET
],
671 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
673 numsize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
681 while (exponent
!= 0);
684 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
687 /* Determine how many bits of the result we already have. */
688 count_leading_zeros (bits
, num
[numsize
- 1]);
689 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
691 /* Now we know the exponent of the number in base two.
692 Check it against the maximum possible exponent. */
696 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
699 /* We have already the first BITS bits of the result. Together with
700 the information whether more non-zero bits follow this is enough
701 to determine the result. */
705 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
706 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
707 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
709 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
713 memcpy (retval
, &num
[least_idx
],
714 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
717 for (i
= least_idx
; i
< numsize
- 1; ++i
)
718 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
720 << (BITS_PER_MP_LIMB
- least_bit
));
721 if (i
- least_idx
< RETURN_LIMB_SIZE
)
722 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
725 /* Check whether any limb beside the ones in RETVAL are non-zero. */
726 for (i
= 0; num
[i
] == 0; ++i
)
729 return round_and_return (retval
, bits
- 1, negative
,
730 num
[round_idx
], round_bit
,
731 int_no
< dig_no
|| i
< round_idx
);
734 else if (dig_no
== int_no
)
736 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
737 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
739 if (target_bit
== is_bit
)
741 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
742 numsize
* sizeof (mp_limb_t
));
743 /* FIXME: the following loop can be avoided if we assume a
744 maximal MANT_DIG value. */
745 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
747 else if (target_bit
> is_bit
)
749 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
750 num
, numsize
, target_bit
- is_bit
);
751 /* FIXME: the following loop can be avoided if we assume a
752 maximal MANT_DIG value. */
753 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
758 assert (numsize
< RETURN_LIMB_SIZE
);
760 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
761 num
, numsize
, is_bit
- target_bit
);
762 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
763 /* FIXME: the following loop can be avoided if we assume a
764 maximal MANT_DIG value. */
765 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
768 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
772 /* Store the bits we already have. */
773 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
774 #if RETURN_LIMB_SIZE > 1
775 if (numsize
< RETURN_LIMB_SIZE
)
780 /* We have to compute at least some of the fractional digits. */
782 /* We construct a fraction and the result of the division gives us
783 the needed digits. The denominator is 1.0 multiplied by the
784 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
785 123e-6 gives 123 / 1000000. */
792 mp_limb_t
*psrc
= den
;
793 mp_limb_t
*pdest
= num
;
794 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
796 assert (dig_no
> int_no
&& exponent
<= 0);
799 /* For the fractional part we need not process too much digits. One
800 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
802 digits we should have enough bits for the result. The remaining
803 decimal digits give us the information that more bits are following.
804 This can be used while rounding. (One added as a safety margin.) */
805 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 1)
807 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 1;
813 neg_exp
= dig_no
- int_no
- exponent
;
815 /* Construct the denominator. */
820 if ((neg_exp
& expbit
) != 0)
827 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
828 memcpy (psrc
, &ttab
->array
[_FPIO_CONST_OFFSET
],
829 densize
* sizeof (mp_limb_t
));
833 cy
= __mpn_mul (pdest
, &ttab
->array
[_FPIO_CONST_OFFSET
],
834 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
836 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
845 while (neg_exp
!= 0);
848 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
850 /* Read the fractional digits from the string. */
851 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
);
854 /* We now have to shift both numbers so that the highest bit in the
855 denominator is set. In the same process we copy the numerator to
856 a high place in the array so that the division constructs the wanted
857 digits. This is done by a "quasi fix point" number representation.
859 num: ddddddddddd . 0000000000000000000000
861 den: ddddddddddd n >= m
865 count_leading_zeros (cnt
, den
[densize
- 1]);
867 (void) __mpn_lshift (den
, den
, densize
, cnt
);
868 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
872 /* Now we are ready for the division. But it is not necessary to
873 do a full multi-precision division because we only need a small
874 number of bits for the result. So we do not use __mpn_divmod
875 here but instead do the division here by hand and stop whenever
876 the needed number of bits is reached. The code itself comes
877 from the GNU MP Library by Torbj\"orn Granlund. */
885 mp_limb_t d
, n
, quot
;
890 assert (numsize
== 1 && n
< d
);
894 udiv_qrnnd (quot
, n
, n
, 0, d
);
901 cnt = BITS_PER_MP_LIMB; \
903 count_leading_zeros (cnt, quot); \
905 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
907 used = MANT_DIG + cnt; \
908 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
909 bits = MANT_DIG + 1; \
913 /* Note that we only clear the second element. */ \
914 /* The conditional is determined at compile time. */ \
915 if (RETURN_LIMB_SIZE > 1) \
921 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
922 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
926 used = MANT_DIG - bits; \
928 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
930 bits += BITS_PER_MP_LIMB
934 while (bits
<= MANT_DIG
);
936 return round_and_return (retval
, exponent
- 1, negative
,
937 quot
, BITS_PER_MP_LIMB
- 1 - used
,
938 more_bits
|| n
!= 0);
942 mp_limb_t d0
, d1
, n0
, n1
;
949 if (numsize
< densize
)
953 /* The numerator of the number occupies fewer bits than
954 the denominator but the one limb is bigger than the
955 high limb of the numerator. */
962 exponent
-= BITS_PER_MP_LIMB
;
965 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
966 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
967 BITS_PER_MP_LIMB
, 0);
970 used
= MANT_DIG
- bits
;
972 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
974 bits
+= BITS_PER_MP_LIMB
;
986 while (bits
<= MANT_DIG
)
992 /* QUOT should be either 111..111 or 111..110. We need
993 special treatment of this rare case as normal division
994 would give overflow. */
995 quot
= ~(mp_limb_t
) 0;
998 if (r
< d1
) /* Carry in the addition? */
1000 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1003 n1
= d0
- (d0
!= 0);
1008 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1009 umul_ppmm (n1
, n0
, d0
, quot
);
1013 if (n1
> r
|| (n1
== r
&& n0
> 0))
1015 /* The estimated QUOT was too large. */
1018 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1020 if (r
>= d1
) /* If not carry, test QUOT again. */
1023 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1029 return round_and_return (retval
, exponent
- 1, negative
,
1030 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1031 more_bits
|| n1
!= 0 || n0
!= 0);
1036 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1040 dX
= den
[densize
- 1];
1041 d1
= den
[densize
- 2];
1043 /* The division does not work if the upper limb of the two-limb
1044 numerator is greater than the denominator. */
1045 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1048 if (numsize
< densize
)
1050 mp_size_t empty
= densize
- numsize
;
1055 for (i
= numsize
; i
> 0; --i
)
1056 num
[i
+ empty
] = num
[i
- 1];
1057 MPN_ZERO (num
, empty
+ 1);
1058 exponent
-= empty
* BITS_PER_MP_LIMB
;
1062 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1064 /* We make a difference here because the compiler
1065 cannot optimize the `else' case that good and
1066 this reflects all currently used FLOAT types
1067 and GMP implementations. */
1069 #if RETURN_LIMB_SIZE <= 2
1070 assert (empty
== 1);
1071 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1072 BITS_PER_MP_LIMB
, 0);
1074 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1075 retval
[i
] = retval
[i
- empty
];
1078 for (i
= numsize
; i
> 0; --i
)
1079 num
[i
+ empty
] = num
[i
- 1];
1080 MPN_ZERO (num
, empty
+ 1);
1084 used
= MANT_DIG
- bits
;
1085 if (used
>= BITS_PER_MP_LIMB
)
1088 (void) __mpn_lshift (&retval
[used
1089 / BITS_PER_MP_LIMB
],
1090 retval
, RETURN_LIMB_SIZE
,
1091 used
% BITS_PER_MP_LIMB
);
1092 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1096 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1098 bits
+= empty
* BITS_PER_MP_LIMB
;
1104 assert (numsize
== densize
);
1105 for (i
= numsize
; i
> 0; --i
)
1106 num
[i
] = num
[i
- 1];
1112 while (bits
<= MANT_DIG
)
1115 /* This might over-estimate QUOT, but it's probably not
1116 worth the extra code here to find out. */
1117 quot
= ~(mp_limb_t
) 0;
1122 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1123 umul_ppmm (n1
, n0
, d1
, quot
);
1125 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1129 if (r
< dX
) /* I.e. "carry in previous addition?" */
1136 /* Possible optimization: We already have (q * n0) and (1 * n1)
1137 after the calculation of QUOT. Taking advantage of this, we
1138 could make this loop make two iterations less. */
1140 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1142 if (num
[densize
] != cy
)
1144 cy
= __mpn_add_n (num
, num
, den
, densize
);
1148 n0
= num
[densize
] = num
[densize
- 1];
1149 for (i
= densize
- 1; i
> 0; --i
)
1150 num
[i
] = num
[i
- 1];
1155 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1157 return round_and_return (retval
, exponent
- 1, negative
,
1158 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1159 more_bits
|| i
>= 0);
1167 /* External user entry point. */
1170 #ifdef weak_function
1173 STRTOF (nptr
, endptr
)
1174 const STRING_TYPE
*nptr
;
1175 STRING_TYPE
**endptr
;
1177 return INTERNAL (STRTOF
) (nptr
, endptr
, 0);