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