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