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