1/* Floating point output for `printf'.
2 Copyright (C) 1995-2012 Free Software Foundation, Inc.
3
4 This file is part of the GNU C Library.
5 Written by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
6
7 The GNU C Library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public
9 License as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11
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 Lesser General Public License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
18 License along with the GNU C Library; if not, see
19 <http://www.gnu.org/licenses/>. */
20
21#include <config.h>
22#include <float.h>
23#include <limits.h>
24#include <math.h>
25#include <string.h>
26#include <unistd.h>
27#include <stdlib.h>
28#include <stdbool.h>
29#define NDEBUG
30#include <assert.h>
31#ifdef HAVE_ERRNO_H
32#include <errno.h>
33#endif
34#include <stdio.h>
35#include <stdarg.h>
36#ifdef HAVE_FENV_H
37#include "quadmath-rounding-mode.h"
38#endif
39#include "quadmath-printf.h"
40#include "fpioconst.h"
41
42#ifdef USE_I18N_NUMBER_H
43#include "_i18n_number.h"
44#endif
45
46
47/* Macros for doing the actual output. */
48
49#define outchar(ch) \
50 do \
51 { \
52 register const int outc = (ch); \
53 if (PUTC (outc, fp) == EOF) \
54 { \
55 if (buffer_malloced) \
56 free (wbuffer); \
57 return -1; \
58 } \
59 ++done; \
60 } while (0)
61
62#define PRINT(ptr, wptr, len) \
63 do \
64 { \
65 register size_t outlen = (len); \
66 if (len > 20) \
67 { \
68 if (PUT (fp, wide ? (const char *) wptr : ptr, outlen) != outlen) \
69 { \
70 if (buffer_malloced) \
71 free (wbuffer); \
72 return -1; \
73 } \
74 ptr += outlen; \
75 done += outlen; \
76 } \
77 else \
78 { \
79 if (wide) \
80 while (outlen-- > 0) \
81 outchar (*wptr++); \
82 else \
83 while (outlen-- > 0) \
84 outchar (*ptr++); \
85 } \
86 } while (0)
87
88#define PADN(ch, len) \
89 do \
90 { \
91 if (PAD (fp, ch, len) != len) \
92 { \
93 if (buffer_malloced) \
94 free (wbuffer); \
95 return -1; \
96 } \
97 done += len; \
98 } \
99 while (0)
100
101
102/* We use the GNU MP library to handle large numbers.
103
104 An MP variable occupies a varying number of entries in its array. We keep
105 track of this number for efficiency reasons. Otherwise we would always
106 have to process the whole array. */
107#define MPN_VAR(name) mp_limb_t *name; mp_size_t name##size
108
109#define MPN_ASSIGN(dst,src) \
110 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
111#define MPN_GE(u,v) \
112 (u##size > v##size || (u##size == v##size && mpn_cmp (u, v, u##size) >= 0))
113
114extern mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size,
115 int *expt, int *is_neg,
116 __float128 value) attribute_hidden;
117static unsigned int guess_grouping (unsigned int intdig_max,
118 const char *grouping);
119
120
121static wchar_t *group_number (wchar_t *buf, wchar_t *bufend,
122 unsigned int intdig_no, const char *grouping,
123 wchar_t thousands_sep, int ngroups);
124
125
126int
127__quadmath_printf_fp (struct __quadmath_printf_file *fp,
128 const struct printf_info *info,
129 const void *const *args)
130{
131 /* The floating-point value to output. */
132 __float128 fpnum;
133
134 /* Locale-dependent representation of decimal point. */
135 const char *decimal;
136 wchar_t decimalwc;
137
138 /* Locale-dependent thousands separator and grouping specification. */
139 const char *thousands_sep = NULL;
140 wchar_t thousands_sepwc = L_('\0');
141 const char *grouping;
142
143 /* "NaN" or "Inf" for the special cases. */
144 const char *special = NULL;
145 const wchar_t *wspecial = NULL;
146
147 /* We need just a few limbs for the input before shifting to the right
148 position. */
149 mp_limb_t fp_input[(FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB];
150 /* We need to shift the contents of fp_input by this amount of bits. */
151 int to_shift = 0;
152
153 /* The fraction of the floting-point value in question */
154 MPN_VAR(frac);
155 /* and the exponent. */
156 int exponent;
157 /* Sign of the exponent. */
158 int expsign = 0;
159 /* Sign of float number. */
160 int is_neg = 0;
161
162 /* Scaling factor. */
163 MPN_VAR(scale);
164
165 /* Temporary bignum value. */
166 MPN_VAR(tmp);
167
168 /* Digit which is result of last hack_digit() call. */
169 wchar_t last_digit, next_digit;
170 bool more_bits;
171
172 /* The type of output format that will be used: 'e'/'E' or 'f'. */
173 int type;
174
175 /* Counter for number of written characters. */
176 int done = 0;
177
178 /* General helper (carry limb). */
179 mp_limb_t cy;
180
181 /* Nonzero if this is output on a wide character stream. */
182 int wide = info->wide;
183
184 /* Buffer in which we produce the output. */
185 wchar_t *wbuffer = NULL;
186 /* Flag whether wbuffer is malloc'ed or not. */
187 int buffer_malloced = 0;
188
189 auto wchar_t hack_digit (void);
190
191 wchar_t hack_digit (void)
192 {
193 mp_limb_t hi;
194
195 if (expsign != 0 && type == 'f' && exponent-- > 0)
196 hi = 0;
197 else if (scalesize == 0)
198 {
199 hi = frac[fracsize - 1];
200 frac[fracsize - 1] = mpn_mul_1 (frac, frac, fracsize - 1, 10);
201 }
202 else
203 {
204 if (fracsize < scalesize)
205 hi = 0;
206 else
207 {
208 hi = mpn_divmod (tmp, frac, fracsize, scale, scalesize);
209 tmp[fracsize - scalesize] = hi;
210 hi = tmp[0];
211
212 fracsize = scalesize;
213 while (fracsize != 0 && frac[fracsize - 1] == 0)
214 --fracsize;
215 if (fracsize == 0)
216 {
217 /* We're not prepared for an mpn variable with zero
218 limbs. */
219 fracsize = 1;
220 return L_('0') + hi;
221 }
222 }
223
224 mp_limb_t _cy = mpn_mul_1 (frac, frac, fracsize, 10);
225 if (_cy != 0)
226 frac[fracsize++] = _cy;
227 }
228
229 return L_('0') + hi;
230 }
231
232 /* Figure out the decimal point character. */
233#ifdef USE_NL_LANGINFO
234 if (info->extra == 0)
235 decimal = nl_langinfo (DECIMAL_POINT);
236 else
237 {
238 decimal = nl_langinfo (MON_DECIMAL_POINT);
239 if (*decimal == '\0')
240 decimal = nl_langinfo (DECIMAL_POINT);
241 }
242 /* The decimal point character must never be zero. */
243 assert (*decimal != '\0');
244#elif defined USE_LOCALECONV
245 const struct lconv *lc = localeconv ();
246 if (info->extra == 0)
247 decimal = lc->decimal_point;
248 else
249 {
250 decimal = lc->mon_decimal_point;
251 if (decimal == NULL || *decimal == '\0')
252 decimal = lc->decimal_point;
253 }
254 if (decimal == NULL || *decimal == '\0')
255 decimal = ".";
256#else
257 decimal = ".";
258#endif
259#ifdef USE_NL_LANGINFO_WC
260 if (info->extra == 0)
261 decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
262 else
263 {
264 decimalwc = nl_langinfo_wc (_NL_MONETARY_DECIMAL_POINT_WC);
265 if (decimalwc == L_('\0'))
266 decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
267 }
268 /* The decimal point character must never be zero. */
269 assert (decimalwc != L_('\0'));
270#else
271 decimalwc = L_('.');
272#endif
273
274#if defined USE_NL_LANGINFO && defined USE_NL_LANGINFO_WC
275 if (info->group)
276 {
277 if (info->extra == 0)
278 grouping = nl_langinfo (GROUPING);
279 else
280 grouping = nl_langinfo (MON_GROUPING);
281
282 if (*grouping <= 0 || *grouping == CHAR_MAX)
283 grouping = NULL;
284 else
285 {
286 /* Figure out the thousands separator character. */
287 if (wide)
288 {
289 if (info->extra == 0)
290 thousands_sepwc = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC);
291 else
292 thousands_sepwc = nl_langinfo_wc (_NL_MONETARY_THOUSANDS_SEP_WC);
293
294 if (thousands_sepwc == L_('\0'))
295 grouping = NULL;
296 }
297 else
298 {
299 if (info->extra == 0)
300 thousands_sep = nl_langinfo (THOUSANDS_SEP);
301 else
302 thousands_sep = nl_langinfo (MON_THOUSANDS_SEP);
303 if (*thousands_sep == '\0')
304 grouping = NULL;
305 }
306 }
307 }
308 else
309#elif defined USE_NL_LANGINFO
310 if (info->group && !wide)
311 {
312 if (info->extra == 0)
313 grouping = nl_langinfo (GROUPING);
314 else
315 grouping = nl_langinfo (MON_GROUPING);
316
317 if (*grouping <= 0 || *grouping == CHAR_MAX)
318 grouping = NULL;
319 else
320 {
321 /* Figure out the thousands separator character. */
322 if (info->extra == 0)
323 thousands_sep = nl_langinfo (THOUSANDS_SEP);
324 else
325 thousands_sep = nl_langinfo (MON_THOUSANDS_SEP);
326
327 if (*thousands_sep == '\0')
328 grouping = NULL;
329 }
330 }
331 else
332#elif defined USE_LOCALECONV
333 if (info->group && !wide)
334 {
335 if (info->extra == 0)
336 grouping = lc->grouping;
337 else
338 grouping = lc->mon_grouping;
339
340 if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX)
341 grouping = NULL;
342 else
343 {
344 /* Figure out the thousands separator character. */
345 if (info->extra == 0)
346 thousands_sep = lc->thousands_sep;
347 else
348 thousands_sep = lc->mon_thousands_sep;
349
350 if (thousands_sep == NULL || *thousands_sep == '\0')
351 grouping = NULL;
352 }
353 }
354 else
355#endif
356 grouping = NULL;
357 if (grouping != NULL && !wide)
358 /* If we are printing multibyte characters and there is a
359 multibyte representation for the thousands separator,
360 we must ensure the wide character thousands separator
361 is available, even if it is fake. */
362 thousands_sepwc = (wchar_t) 0xfffffffe;
363
364 /* Fetch the argument value. */
365 {
366 fpnum = **(const __float128 **) args[0];
367
368 /* Check for special values: not a number or infinity. */
369 if (isnanq (fpnum))
370 {
371 ieee854_float128 u = { .value = fpnum };
372 is_neg = u.ieee.negative != 0;
373 if (isupper (info->spec))
374 {
375 special = "NAN";
376 wspecial = L_("NAN");
377 }
378 else
379 {
380 special = "nan";
381 wspecial = L_("nan");
382 }
383 }
384 else if (isinfq (fpnum))
385 {
386 is_neg = fpnum < 0;
387 if (isupper (info->spec))
388 {
389 special = "INF";
390 wspecial = L_("INF");
391 }
392 else
393 {
394 special = "inf";
395 wspecial = L_("inf");
396 }
397 }
398 else
399 {
400 fracsize = mpn_extract_flt128 (fp_input,
401 (sizeof (fp_input) /
402 sizeof (fp_input[0])),
403 &exponent, &is_neg, fpnum);
404 to_shift = 1 + fracsize * BITS_PER_MP_LIMB - FLT128_MANT_DIG;
405 }
406 }
407
408 if (special)
409 {
410 int width = info->width;
411
412 if (is_neg || info->showsign || info->space)
413 --width;
414 width -= 3;
415
416 if (!info->left && width > 0)
417 PADN (' ', width);
418
419 if (is_neg)
420 outchar ('-');
421 else if (info->showsign)
422 outchar ('+');
423 else if (info->space)
424 outchar (' ');
425
426 PRINT (special, wspecial, 3);
427
428 if (info->left && width > 0)
429 PADN (' ', width);
430
431 return done;
432 }
433
434
435 /* We need three multiprecision variables. Now that we have the exponent
436 of the number we can allocate the needed memory. It would be more
437 efficient to use variables of the fixed maximum size but because this
438 would be really big it could lead to memory problems. */
439 {
440 mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1)
441 / BITS_PER_MP_LIMB
442 + (FLT128_MANT_DIG / BITS_PER_MP_LIMB > 2 ? 8 : 4))
443 * sizeof (mp_limb_t);
444 frac = (mp_limb_t *) alloca (bignum_size);
445 tmp = (mp_limb_t *) alloca (bignum_size);
446 scale = (mp_limb_t *) alloca (bignum_size);
447 }
448
449 /* We now have to distinguish between numbers with positive and negative
450 exponents because the method used for the one is not applicable/efficient
451 for the other. */
452 scalesize = 0;
453 if (exponent > 2)
454 {
455 /* |FP| >= 8.0. */
456 int scaleexpo = 0;
457 int explog = FLT128_MAX_10_EXP_LOG;
458 int exp10 = 0;
459 const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
460 int cnt_h, cnt_l, i;
461
462 if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0)
463 {
464 MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
465 fp_input, fracsize);
466 fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
467 }
468 else
469 {
470 cy = mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
471 fp_input, fracsize,
472 (exponent + to_shift) % BITS_PER_MP_LIMB);
473 fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
474 if (cy)
475 frac[fracsize++] = cy;
476 }
477 MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB);
478
479 assert (powers > &_fpioconst_pow10[0]);
480 do
481 {
482 --powers;
483
484 /* The number of the product of two binary numbers with n and m
485 bits respectively has m+n or m+n-1 bits. */
486 if (exponent >= scaleexpo + powers->p_expo - 1)
487 {
488 if (scalesize == 0)
489 {
490 if (FLT128_MANT_DIG > _FPIO_CONST_OFFSET * BITS_PER_MP_LIMB)
491 {
492#define _FPIO_CONST_SHIFT \
493 (((FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) \
494 - _FPIO_CONST_OFFSET)
495 /* 64bit const offset is not enough for
496 IEEE quad long double. */
497 tmpsize = powers->arraysize + _FPIO_CONST_SHIFT;
498 memcpy (tmp + _FPIO_CONST_SHIFT,
499 &__tens[powers->arrayoff],
500 tmpsize * sizeof (mp_limb_t));
501 MPN_ZERO (tmp, _FPIO_CONST_SHIFT);
502 /* Adjust exponent, as scaleexpo will be this much
503 bigger too. */
504 exponent += _FPIO_CONST_SHIFT * BITS_PER_MP_LIMB;
505 }
506 else
507 {
508 tmpsize = powers->arraysize;
509 memcpy (tmp, &__tens[powers->arrayoff],
510 tmpsize * sizeof (mp_limb_t));
511 }
512 }
513 else
514 {
515 cy = mpn_mul (tmp, scale, scalesize,
516 &__tens[powers->arrayoff
517 + _FPIO_CONST_OFFSET],
518 powers->arraysize - _FPIO_CONST_OFFSET);
519 tmpsize = scalesize + powers->arraysize - _FPIO_CONST_OFFSET;
520 if (cy == 0)
521 --tmpsize;
522 }
523
524 if (MPN_GE (frac, tmp))
525 {
526 int cnt;
527 MPN_ASSIGN (scale, tmp);
528 count_leading_zeros (cnt, scale[scalesize - 1]);
529 scaleexpo = (scalesize - 2) * BITS_PER_MP_LIMB - cnt - 1;
530 exp10 |= 1 << explog;
531 }
532 }
533 --explog;
534 }
535 while (powers > &_fpioconst_pow10[0]);
536 exponent = exp10;
537
538 /* Optimize number representations. We want to represent the numbers
539 with the lowest number of bytes possible without losing any
540 bytes. Also the highest bit in the scaling factor has to be set
541 (this is a requirement of the MPN division routines). */
542 if (scalesize > 0)
543 {
544 /* Determine minimum number of zero bits at the end of
545 both numbers. */
546 for (i = 0; scale[i] == 0 && frac[i] == 0; i++)
547 ;
548
549 /* Determine number of bits the scaling factor is misplaced. */
550 count_leading_zeros (cnt_h, scale[scalesize - 1]);
551
552 if (cnt_h == 0)
553 {
554 /* The highest bit of the scaling factor is already set. So
555 we only have to remove the trailing empty limbs. */
556 if (i > 0)
557 {
558 MPN_COPY_INCR (scale, scale + i, scalesize - i);
559 scalesize -= i;
560 MPN_COPY_INCR (frac, frac + i, fracsize - i);
561 fracsize -= i;
562 }
563 }
564 else
565 {
566 if (scale[i] != 0)
567 {
568 count_trailing_zeros (cnt_l, scale[i]);
569 if (frac[i] != 0)
570 {
571 int cnt_l2;
572 count_trailing_zeros (cnt_l2, frac[i]);
573 if (cnt_l2 < cnt_l)
574 cnt_l = cnt_l2;
575 }
576 }
577 else
578 count_trailing_zeros (cnt_l, frac[i]);
579
580 /* Now shift the numbers to their optimal position. */
581 if (i == 0 && BITS_PER_MP_LIMB - cnt_h > cnt_l)
582 {
583 /* We cannot save any memory. So just roll both numbers
584 so that the scaling factor has its highest bit set. */
585
586 (void) mpn_lshift (scale, scale, scalesize, cnt_h);
587 cy = mpn_lshift (frac, frac, fracsize, cnt_h);
588 if (cy != 0)
589 frac[fracsize++] = cy;
590 }
591 else if (BITS_PER_MP_LIMB - cnt_h <= cnt_l)
592 {
593 /* We can save memory by removing the trailing zero limbs
594 and by packing the non-zero limbs which gain another
595 free one. */
596
597 (void) mpn_rshift (scale, scale + i, scalesize - i,
598 BITS_PER_MP_LIMB - cnt_h);
599 scalesize -= i + 1;
600 (void) mpn_rshift (frac, frac + i, fracsize - i,
601 BITS_PER_MP_LIMB - cnt_h);
602 fracsize -= frac[fracsize - i - 1] == 0 ? i + 1 : i;
603 }
604 else
605 {
606 /* We can only save the memory of the limbs which are zero.
607 The non-zero parts occupy the same number of limbs. */
608
609 (void) mpn_rshift (scale, scale + (i - 1),
610 scalesize - (i - 1),
611 BITS_PER_MP_LIMB - cnt_h);
612 scalesize -= i;
613 (void) mpn_rshift (frac, frac + (i - 1),
614 fracsize - (i - 1),
615 BITS_PER_MP_LIMB - cnt_h);
616 fracsize -= frac[fracsize - (i - 1) - 1] == 0 ? i : i - 1;
617 }
618 }
619 }
620 }
621 else if (exponent < 0)
622 {
623 /* |FP| < 1.0. */
624 int exp10 = 0;
625 int explog = FLT128_MAX_10_EXP_LOG;
626 const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
627
628 /* Now shift the input value to its right place. */
629 cy = mpn_lshift (frac, fp_input, fracsize, to_shift);
630 frac[fracsize++] = cy;
631 assert (cy == 1 || (frac[fracsize - 2] == 0 && frac[0] == 0));
632
633 expsign = 1;
634 exponent = -exponent;
635
636 assert (powers != &_fpioconst_pow10[0]);
637 do
638 {
639 --powers;
640
641 if (exponent >= powers->m_expo)
642 {
643 int i, incr, cnt_h, cnt_l;
644 mp_limb_t topval[2];
645
646 /* The mpn_mul function expects the first argument to be
647 bigger than the second. */
648 if (fracsize < powers->arraysize - _FPIO_CONST_OFFSET)
649 cy = mpn_mul (tmp, &__tens[powers->arrayoff
650 + _FPIO_CONST_OFFSET],
651 powers->arraysize - _FPIO_CONST_OFFSET,
652 frac, fracsize);
653 else
654 cy = mpn_mul (tmp, frac, fracsize,
655 &__tens[powers->arrayoff + _FPIO_CONST_OFFSET],
656 powers->arraysize - _FPIO_CONST_OFFSET);
657 tmpsize = fracsize + powers->arraysize - _FPIO_CONST_OFFSET;
658 if (cy == 0)
659 --tmpsize;
660
661 count_leading_zeros (cnt_h, tmp[tmpsize - 1]);
662 incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB
663 + BITS_PER_MP_LIMB - 1 - cnt_h;
664
665 assert (incr <= powers->p_expo);
666
667 /* If we increased the exponent by exactly 3 we have to test
668 for overflow. This is done by comparing with 10 shifted
669 to the right position. */
670 if (incr == exponent + 3)
671 {
672 if (cnt_h <= BITS_PER_MP_LIMB - 4)
673 {
674 topval[0] = 0;
675 topval[1]
676 = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4 - cnt_h);
677 }
678 else
679 {
680 topval[0] = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4);
681 topval[1] = 0;
682 (void) mpn_lshift (topval, topval, 2,
683 BITS_PER_MP_LIMB - cnt_h);
684 }
685 }
686
687 /* We have to be careful when multiplying the last factor.
688 If the result is greater than 1.0 be have to test it
689 against 10.0. If it is greater or equal to 10.0 the
690 multiplication was not valid. This is because we cannot
691 determine the number of bits in the result in advance. */
692 if (incr < exponent + 3
693 || (incr == exponent + 3 &&
694 (tmp[tmpsize - 1] < topval[1]
695 || (tmp[tmpsize - 1] == topval[1]
696 && tmp[tmpsize - 2] < topval[0]))))
697 {
698 /* The factor is right. Adapt binary and decimal
699 exponents. */
700 exponent -= incr;
701 exp10 |= 1 << explog;
702
703 /* If this factor yields a number greater or equal to
704 1.0, we must not shift the non-fractional digits down. */
705 if (exponent < 0)
706 cnt_h += -exponent;
707
708 /* Now we optimize the number representation. */
709 for (i = 0; tmp[i] == 0; ++i);
710 if (cnt_h == BITS_PER_MP_LIMB - 1)
711 {
712 MPN_COPY (frac, tmp + i, tmpsize - i);
713 fracsize = tmpsize - i;
714 }
715 else
716 {
717 count_trailing_zeros (cnt_l, tmp[i]);
718
719 /* Now shift the numbers to their optimal position. */
720 if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l)
721 {
722 /* We cannot save any memory. Just roll the
723 number so that the leading digit is in a
724 separate limb. */
725
726 cy = mpn_lshift (frac, tmp, tmpsize, cnt_h + 1);
727 fracsize = tmpsize + 1;
728 frac[fracsize - 1] = cy;
729 }
730 else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l)
731 {
732 (void) mpn_rshift (frac, tmp + i, tmpsize - i,
733 BITS_PER_MP_LIMB - 1 - cnt_h);
734 fracsize = tmpsize - i;
735 }
736 else
737 {
738 /* We can only save the memory of the limbs which
739 are zero. The non-zero parts occupy the same
740 number of limbs. */
741
742 (void) mpn_rshift (frac, tmp + (i - 1),
743 tmpsize - (i - 1),
744 BITS_PER_MP_LIMB - 1 - cnt_h);
745 fracsize = tmpsize - (i - 1);
746 }
747 }
748 }
749 }
750 --explog;
751 }
752 while (powers != &_fpioconst_pow10[1] && exponent > 0);
753 /* All factors but 10^-1 are tested now. */
754 if (exponent > 0)
755 {
756 int cnt_l;
757
758 cy = mpn_mul_1 (tmp, frac, fracsize, 10);
759 tmpsize = fracsize;
760 assert (cy == 0 || tmp[tmpsize - 1] < 20);
761
762 count_trailing_zeros (cnt_l, tmp[0]);
763 if (cnt_l < MIN (4, exponent))
764 {
765 cy = mpn_lshift (frac, tmp, tmpsize,
766 BITS_PER_MP_LIMB - MIN (4, exponent));
767 if (cy != 0)
768 frac[tmpsize++] = cy;
769 }
770 else
771 (void) mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent));
772 fracsize = tmpsize;
773 exp10 |= 1;
774 assert (frac[fracsize - 1] < 10);
775 }
776 exponent = exp10;
777 }
778 else
779 {
780 /* This is a special case. We don't need a factor because the
781 numbers are in the range of 1.0 <= |fp| < 8.0. We simply
782 shift it to the right place and divide it by 1.0 to get the
783 leading digit. (Of course this division is not really made.) */
784 assert (0 <= exponent && exponent < 3 &&
785 exponent + to_shift < BITS_PER_MP_LIMB);
786
787 /* Now shift the input value to its right place. */
788 cy = mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift));
789 frac[fracsize++] = cy;
790 exponent = 0;
791 }
792
793 {
794 int width = info->width;
795 wchar_t *wstartp, *wcp;
796 size_t chars_needed;
797 int expscale;
798 int intdig_max, intdig_no = 0;
799 int fracdig_min;
800 int fracdig_max;
801 int dig_max;
802 int significant;
803 int ngroups = 0;
804 char spec = tolower (info->spec);
805
806 if (spec == 'e')
807 {
808 type = info->spec;
809 intdig_max = 1;
810 fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
811 chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4;
812 /* d . ddd e +- ddd */
813 dig_max = __INT_MAX__; /* Unlimited. */
814 significant = 1; /* Does not matter here. */
815 }
816 else if (spec == 'f')
817 {
818 type = 'f';
819 fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
820 dig_max = __INT_MAX__; /* Unlimited. */
821 significant = 1; /* Does not matter here. */
822 if (expsign == 0)
823 {
824 intdig_max = exponent + 1;
825 /* This can be really big! */ /* XXX Maybe malloc if too big? */
826 chars_needed = (size_t) exponent + 1 + 1 + (size_t) fracdig_max;
827 }
828 else
829 {
830 intdig_max = 1;
831 chars_needed = 1 + 1 + (size_t) fracdig_max;
832 }
833 }
834 else
835 {
836 dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec);
837 if ((expsign == 0 && exponent >= dig_max)
838 || (expsign != 0 && exponent > 4))
839 {
840 if ('g' - 'G' == 'e' - 'E')
841 type = 'E' + (info->spec - 'G');
842 else
843 type = isupper (info->spec) ? 'E' : 'e';
844 fracdig_max = dig_max - 1;
845 intdig_max = 1;
846 chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4;
847 }
848 else
849 {
850 type = 'f';
851 intdig_max = expsign == 0 ? exponent + 1 : 0;
852 fracdig_max = dig_max - intdig_max;
853 /* We need space for the significant digits and perhaps
854 for leading zeros when < 1.0. The number of leading
855 zeros can be as many as would be required for
856 exponential notation with a negative two-digit
857 exponent, which is 4. */
858 chars_needed = (size_t) dig_max + 1 + 4;
859 }
860 fracdig_min = info->alt ? fracdig_max : 0;
861 significant = 0; /* We count significant digits. */
862 }
863
864 if (grouping)
865 {
866 /* Guess the number of groups we will make, and thus how
867 many spaces we need for separator characters. */
868 ngroups = guess_grouping (intdig_max, grouping);
869 /* Allocate one more character in case rounding increases the
870 number of groups. */
871 chars_needed += ngroups + 1;
872 }
873
874 /* Allocate buffer for output. We need two more because while rounding
875 it is possible that we need two more characters in front of all the
876 other output. If the amount of memory we have to allocate is too
877 large use `malloc' instead of `alloca'. */
878 if (__builtin_expect (chars_needed >= (size_t) -1 / sizeof (wchar_t) - 2
879 || chars_needed < fracdig_max, 0))
880 {
881 /* Some overflow occurred. */
882#if defined HAVE_ERRNO_H && defined ERANGE
883 errno = ERANGE;
884#endif
885 return -1;
886 }
887 size_t wbuffer_to_alloc = (2 + chars_needed) * sizeof (wchar_t);
888 buffer_malloced = wbuffer_to_alloc >= 4096;
889 if (__builtin_expect (buffer_malloced, 0))
890 {
891 wbuffer = (wchar_t *) malloc (wbuffer_to_alloc);
892 if (wbuffer == NULL)
893 /* Signal an error to the caller. */
894 return -1;
895 }
896 else
897 wbuffer = (wchar_t *) alloca (wbuffer_to_alloc);
898 wcp = wstartp = wbuffer + 2; /* Let room for rounding. */
899
900 /* Do the real work: put digits in allocated buffer. */
901 if (expsign == 0 || type != 'f')
902 {
903 assert (expsign == 0 || intdig_max == 1);
904 while (intdig_no < intdig_max)
905 {
906 ++intdig_no;
907 *wcp++ = hack_digit ();
908 }
909 significant = 1;
910 if (info->alt
911 || fracdig_min > 0
912 || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0)))
913 *wcp++ = decimalwc;
914 }
915 else
916 {
917 /* |fp| < 1.0 and the selected type is 'f', so put "0."
918 in the buffer. */
919 *wcp++ = L_('0');
920 --exponent;
921 *wcp++ = decimalwc;
922 }
923
924 /* Generate the needed number of fractional digits. */
925 int fracdig_no = 0;
926 int added_zeros = 0;
927 while (fracdig_no < fracdig_min + added_zeros
928 || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
929 {
930 ++fracdig_no;
931 *wcp = hack_digit ();
932 if (*wcp++ != L_('0'))
933 significant = 1;
934 else if (significant == 0)
935 {
936 ++fracdig_max;
937 if (fracdig_min > 0)
938 ++added_zeros;
939 }
940 }
941
942 /* Do rounding. */
943 last_digit = wcp[-1] != decimalwc ? wcp[-1] : wcp[-2];
944 next_digit =hack_digit ();
945
946 if (next_digit != L_('0') && next_digit != L_('5'))
947 more_bits = true;
948 else if (fracsize == 1 && frac[0] == 0)
949 /* Rest of the number is zero. */
950 more_bits = false;
951 else if (scalesize == 0)
952 {
953 /* Here we have to see whether all limbs are zero since no
954 normalization happened. */
955 size_t lcnt = fracsize;
956 while (lcnt >= 1 && frac[lcnt - 1] == 0)
957 --lcnt;
958 more_bits = lcnt > 0;
959 }
960 else
961 more_bits = true;
962
963#ifdef HAVE_FENV_H
964 int rounding_mode = get_rounding_mode ();
965 if (round_away (is_neg, (last_digit - L_('0')) & 1, next_digit >= L_('5'),
966 more_bits, rounding_mode))
967 {
968 wchar_t *wtp = wcp;
969
970 if (fracdig_no > 0)
971 {
972 /* Process fractional digits. Terminate if not rounded or
973 radix character is reached. */
974 int removed = 0;
975 while (*--wtp != decimalwc && *wtp == L_('9'))
976 {
977 *wtp = L_('0');
978 ++removed;
979 }
980 if (removed == fracdig_min && added_zeros > 0)
981 --added_zeros;
982 if (*wtp != decimalwc)
983 /* Round up. */
984 (*wtp)++;
985 else if (__builtin_expect (spec == 'g' && type == 'f' && info->alt
986 && wtp == wstartp + 1
987 && wstartp[0] == L_('0'),
988 0))
989 /* This is a special case: the rounded number is 1.0,
990 the format is 'g' or 'G', and the alternative format
991 is selected. This means the result must be "1.". */
992 --added_zeros;
993 }
994
995 if (fracdig_no == 0 || *wtp == decimalwc)
996 {
997 /* Round the integer digits. */
998 if (*(wtp - 1) == decimalwc)
999 --wtp;
1000
1001 while (--wtp >= wstartp && *wtp == L_('9'))
1002 *wtp = L_('0');
1003
1004 if (wtp >= wstartp)
1005 /* Round up. */
1006 (*wtp)++;
1007 else
1008 /* It is more critical. All digits were 9's. */
1009 {
1010 if (type != 'f')
1011 {
1012 *wstartp = '1';
1013 exponent += expsign == 0 ? 1 : -1;
1014
1015 /* The above exponent adjustment could lead to 1.0e-00,
1016 e.g. for 0.999999999. Make sure exponent 0 always
1017 uses + sign. */
1018 if (exponent == 0)
1019 expsign = 0;
1020 }
1021 else if (intdig_no == dig_max)
1022 {
1023 /* This is the case where for type %g the number fits
1024 really in the range for %f output but after rounding
1025 the number of digits is too big. */
1026 *--wstartp = decimalwc;
1027 *--wstartp = L_('1');
1028
1029 if (info->alt || fracdig_no > 0)
1030 {
1031 /* Overwrite the old radix character. */
1032 wstartp[intdig_no + 2] = L_('0');
1033 ++fracdig_no;
1034 }
1035
1036 fracdig_no += intdig_no;
1037 intdig_no = 1;
1038 fracdig_max = intdig_max - intdig_no;
1039 ++exponent;
1040 /* Now we must print the exponent. */
1041 type = isupper (info->spec) ? 'E' : 'e';
1042 }
1043 else
1044 {
1045 /* We can simply add another another digit before the
1046 radix. */
1047 *--wstartp = L_('1');
1048 ++intdig_no;
1049 }
1050
1051 /* While rounding the number of digits can change.
1052 If the number now exceeds the limits remove some
1053 fractional digits. */
1054 if (intdig_no + fracdig_no > dig_max)
1055 {
1056 wcp -= intdig_no + fracdig_no - dig_max;
1057 fracdig_no -= intdig_no + fracdig_no - dig_max;
1058 }
1059 }
1060 }
1061 }
1062#endif
1063
1064 /* Now remove unnecessary '0' at the end of the string. */
1065 while (fracdig_no > fracdig_min + added_zeros && *(wcp - 1) == L_('0'))
1066 {
1067 --wcp;
1068 --fracdig_no;
1069 }
1070 /* If we eliminate all fractional digits we perhaps also can remove
1071 the radix character. */
1072 if (fracdig_no == 0 && !info->alt && *(wcp - 1) == decimalwc)
1073 --wcp;
1074
1075 if (grouping)
1076 {
1077 /* Rounding might have changed the number of groups. We allocated
1078 enough memory but we need here the correct number of groups. */
1079 if (intdig_no != intdig_max)
1080 ngroups = guess_grouping (intdig_no, grouping);
1081
1082 /* Add in separator characters, overwriting the same buffer. */
1083 wcp = group_number (wstartp, wcp, intdig_no, grouping, thousands_sepwc,
1084 ngroups);
1085 }
1086
1087 /* Write the exponent if it is needed. */
1088 if (type != 'f')
1089 {
1090 if (__builtin_expect (expsign != 0 && exponent == 4 && spec == 'g', 0))
1091 {
1092 /* This is another special case. The exponent of the number is
1093 really smaller than -4, which requires the 'e'/'E' format.
1094 But after rounding the number has an exponent of -4. */
1095 assert (wcp >= wstartp + 1);
1096 assert (wstartp[0] == L_('1'));
1097 memcpy (wstartp, L_("0.0001"), 6 * sizeof (wchar_t));
1098 wstartp[1] = decimalwc;
1099 if (wcp >= wstartp + 2)
1100 {
1101 size_t cnt;
1102 for (cnt = 0; cnt < wcp - (wstartp + 2); cnt++)
1103 wstartp[6 + cnt] = L_('0');
1104 wcp += 4;
1105 }
1106 else
1107 wcp += 5;
1108 }
1109 else
1110 {
1111 *wcp++ = (wchar_t) type;
1112 *wcp++ = expsign ? L_('-') : L_('+');
1113
1114 /* Find the magnitude of the exponent. */
1115 expscale = 10;
1116 while (expscale <= exponent)
1117 expscale *= 10;
1118
1119 if (exponent < 10)
1120 /* Exponent always has at least two digits. */
1121 *wcp++ = L_('0');
1122 else
1123 do
1124 {
1125 expscale /= 10;
1126 *wcp++ = L_('0') + (exponent / expscale);
1127 exponent %= expscale;
1128 }
1129 while (expscale > 10);
1130 *wcp++ = L_('0') + exponent;
1131 }
1132 }
1133
1134 /* Compute number of characters which must be filled with the padding
1135 character. */
1136 if (is_neg || info->showsign || info->space)
1137 --width;
1138 width -= wcp - wstartp;
1139
1140 if (!info->left && info->pad != '0' && width > 0)
1141 PADN (info->pad, width);
1142
1143 if (is_neg)
1144 outchar ('-');
1145 else if (info->showsign)
1146 outchar ('+');
1147 else if (info->space)
1148 outchar (' ');
1149
1150 if (!info->left && info->pad == '0' && width > 0)
1151 PADN ('0', width);
1152
1153 {
1154 char *buffer = NULL;
1155 char *buffer_end __attribute__((__unused__)) = NULL;
1156 char *cp = NULL;
1157 char *tmpptr;
1158
1159 if (! wide)
1160 {
1161 /* Create the single byte string. */
1162 size_t decimal_len;
1163 size_t thousands_sep_len;
1164 wchar_t *copywc;
1165#ifdef USE_I18N_NUMBER_H
1166 size_t factor = (info->i18n
1167 ? nl_langinfo_wc (_NL_CTYPE_MB_CUR_MAX)
1168 : 1);
1169#else
1170 size_t factor = 1;
1171#endif
1172
1173 decimal_len = strlen (decimal);
1174
1175 if (thousands_sep == NULL)
1176 thousands_sep_len = 0;
1177 else
1178 thousands_sep_len = strlen (thousands_sep);
1179
1180 size_t nbuffer = (2 + chars_needed * factor + decimal_len
1181 + ngroups * thousands_sep_len);
1182 if (__builtin_expect (buffer_malloced, 0))
1183 {
1184 buffer = (char *) malloc (nbuffer);
1185 if (buffer == NULL)
1186 {
1187 /* Signal an error to the caller. */
1188 free (wbuffer);
1189 return -1;
1190 }
1191 }
1192 else
1193 buffer = (char *) alloca (nbuffer);
1194 buffer_end = buffer + nbuffer;
1195
1196 /* Now copy the wide character string. Since the character
1197 (except for the decimal point and thousands separator) must
1198 be coming from the ASCII range we can esily convert the
1199 string without mapping tables. */
1200 for (cp = buffer, copywc = wstartp; copywc < wcp; ++copywc)
1201 if (*copywc == decimalwc)
1202 memcpy (cp, decimal, decimal_len), cp += decimal_len;
1203 else if (*copywc == thousands_sepwc)
1204 memcpy (cp, thousands_sep, thousands_sep_len), cp += thousands_sep_len;
1205 else
1206 *cp++ = (char) *copywc;
1207 }
1208
1209 tmpptr = buffer;
1210#if USE_I18N_NUMBER_H
1211 if (__builtin_expect (info->i18n, 0))
1212 {
1213 tmpptr = _i18n_number_rewrite (tmpptr, cp, buffer_end);
1214 cp = buffer_end;
1215 assert ((uintptr_t) buffer <= (uintptr_t) tmpptr);
1216 assert ((uintptr_t) tmpptr < (uintptr_t) buffer_end);
1217 }
1218#endif
1219
1220 PRINT (tmpptr, wstartp, wide ? wcp - wstartp : cp - tmpptr);
1221
1222 /* Free the memory if necessary. */
1223 if (__builtin_expect (buffer_malloced, 0))
1224 {
1225 free (buffer);
1226 free (wbuffer);
1227 }
1228 }
1229
1230 if (info->left && width > 0)
1231 PADN (info->pad, width);
1232 }
1233 return done;
1234}
1235
1236/* Return the number of extra grouping characters that will be inserted
1237 into a number with INTDIG_MAX integer digits. */
1238
1239static unsigned int
1240guess_grouping (unsigned int intdig_max, const char *grouping)
1241{
1242 unsigned int groups;
1243
1244 /* We treat all negative values like CHAR_MAX. */
1245
1246 if (*grouping == CHAR_MAX || *grouping <= 0)
1247 /* No grouping should be done. */
1248 return 0;
1249
1250 groups = 0;
1251 while (intdig_max > (unsigned int) *grouping)
1252 {
1253 ++groups;
1254 intdig_max -= *grouping++;
1255
1256 if (*grouping == CHAR_MAX
1257#if CHAR_MIN < 0
1258 || *grouping < 0
1259#endif
1260 )
1261 /* No more grouping should be done. */
1262 break;
1263 else if (*grouping == 0)
1264 {
1265 /* Same grouping repeats. */
1266 groups += (intdig_max - 1) / grouping[-1];
1267 break;
1268 }
1269 }
1270
1271 return groups;
1272}
1273
1274/* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
1275 There is guaranteed enough space past BUFEND to extend it.
1276 Return the new end of buffer. */
1277
1278static wchar_t *
1279group_number (wchar_t *buf, wchar_t *bufend, unsigned int intdig_no,
1280 const char *grouping, wchar_t thousands_sep, int ngroups)
1281{
1282 wchar_t *p;
1283
1284 if (ngroups == 0)
1285 return bufend;
1286
1287 /* Move the fractional part down. */
1288 memmove (buf + intdig_no + ngroups, buf + intdig_no,
1289 (bufend - (buf + intdig_no)) * sizeof (wchar_t));
1290
1291 p = buf + intdig_no + ngroups - 1;
1292 do
1293 {
1294 unsigned int len = *grouping++;
1295 do
1296 *p-- = buf[--intdig_no];
1297 while (--len > 0);
1298 *p-- = thousands_sep;
1299
1300 if (*grouping == CHAR_MAX
1301#if CHAR_MIN < 0
1302 || *grouping < 0
1303#endif
1304 )
1305 /* No more grouping should be done. */
1306 break;
1307 else if (*grouping == 0)
1308 /* Same grouping repeats. */
1309 --grouping;
1310 } while (intdig_no > (unsigned int) *grouping);
1311
1312 /* Copy the remaining ungrouped digits. */
1313 do
1314 *p-- = buf[--intdig_no];
1315 while (p > buf);
1316
1317 return bufend + ngroups;
1318}
1319