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