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