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