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