32df06a4ba6234c06eeff7b93e8d5d983c3f6033
[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 0.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     while (fracdig_no < fracdig_min
930            || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
931       {
932         ++fracdig_no;
933         *wcp = hack_digit ();
934         if (*wcp++ != L'0')
935           significant = 1;
936         else if (significant == 0)
937           {
938             ++fracdig_max;
939             if (fracdig_min > 0)
940               ++fracdig_min;
941           }
942       }
943
944     /* Do rounding.  */
945     digit = hack_digit ();
946     if (digit > L'4')
947       {
948         wchar_t *wtp = wcp;
949
950         if (digit == L'5'
951             && ((*(wcp - 1) != decimalwc && (*(wcp - 1) & 1) == 0)
952                 || ((*(wcp - 1) == decimalwc && (*(wcp - 2) & 1) == 0))))
953           {
954             /* This is the critical case.        */
955             if (fracsize == 1 && frac[0] == 0)
956               /* Rest of the number is zero -> round to even.
957                  (IEEE 754-1985 4.1 says this is the default rounding.)  */
958               goto do_expo;
959             else if (scalesize == 0)
960               {
961                 /* Here we have to see whether all limbs are zero since no
962                    normalization happened.  */
963                 size_t lcnt = fracsize;
964                 while (lcnt >= 1 && frac[lcnt - 1] == 0)
965                   --lcnt;
966                 if (lcnt == 0)
967                   /* Rest of the number is zero -> round to even.
968                      (IEEE 754-1985 4.1 says this is the default rounding.)  */
969                   goto do_expo;
970               }
971           }
972
973         if (fracdig_no > 0)
974           {
975             /* Process fractional digits.  Terminate if not rounded or
976                radix character is reached.  */
977             while (*--wtp != decimalwc && *wtp == L'9')
978               *wtp = L'0';
979             if (*wtp != decimalwc)
980               /* Round up.  */
981               (*wtp)++;
982           }
983
984         if (fracdig_no == 0 || *wtp == decimalwc)
985           {
986             /* Round the integer digits.  */
987             if (*(wtp - 1) == decimalwc)
988               --wtp;
989
990             while (--wtp >= wstartp && *wtp == L'9')
991               *wtp = L'0';
992
993             if (wtp >= wstartp)
994               /* Round up.  */
995               (*wtp)++;
996             else
997               /* It is more critical.  All digits were 9's.  */
998               {
999                 if (type != 'f')
1000                   {
1001                     *wstartp = '1';
1002                     exponent += expsign == 0 ? 1 : -1;
1003                   }
1004                 else if (intdig_no == dig_max)
1005                   {
1006                     /* This is the case where for type %g the number fits
1007                        really in the range for %f output but after rounding
1008                        the number of digits is too big.  */
1009                     *--wstartp = decimalwc;
1010                     *--wstartp = L'1';
1011
1012                     if (info->alt || fracdig_no > 0)
1013                       {
1014                         /* Overwrite the old radix character.  */
1015                         wstartp[intdig_no + 2] = L'0';
1016                         ++fracdig_no;
1017                       }
1018
1019                     fracdig_no += intdig_no;
1020                     intdig_no = 1;
1021                     fracdig_max = intdig_max - intdig_no;
1022                     ++exponent;
1023                     /* Now we must print the exponent.  */
1024                     type = isupper (info->spec) ? 'E' : 'e';
1025                   }
1026                 else
1027                   {
1028                     /* We can simply add another another digit before the
1029                        radix.  */
1030                     *--wstartp = L'1';
1031                     ++intdig_no;
1032                   }
1033
1034                 /* While rounding the number of digits can change.
1035                    If the number now exceeds the limits remove some
1036                    fractional digits.  */
1037                 if (intdig_no + fracdig_no > dig_max)
1038                   {
1039                     wcp -= intdig_no + fracdig_no - dig_max;
1040                     fracdig_no -= intdig_no + fracdig_no - dig_max;
1041                   }
1042               }
1043           }
1044       }
1045
1046   do_expo:
1047     /* Now remove unnecessary '0' at the end of the string.  */
1048     while (fracdig_no > fracdig_min && *(wcp - 1) == L'0')
1049       {
1050         --wcp;
1051         --fracdig_no;
1052       }
1053     /* If we eliminate all fractional digits we perhaps also can remove
1054        the radix character.  */
1055     if (fracdig_no == 0 && !info->alt && *(wcp - 1) == decimalwc)
1056       --wcp;
1057
1058     if (grouping)
1059       /* Add in separator characters, overwriting the same buffer.  */
1060       wcp = group_number (wstartp, wcp, intdig_no, grouping, thousands_sepwc,
1061                           ngroups);
1062
1063     /* Write the exponent if it is needed.  */
1064     if (type != 'f')
1065       {
1066         *wcp++ = (wchar_t) type;
1067         *wcp++ = expsign ? L'-' : L'+';
1068
1069         /* Find the magnitude of the exponent.  */
1070         expscale = 10;
1071         while (expscale <= exponent)
1072           expscale *= 10;
1073
1074         if (exponent < 10)
1075           /* Exponent always has at least two digits.  */
1076           *wcp++ = L'0';
1077         else
1078           do
1079             {
1080               expscale /= 10;
1081               *wcp++ = L'0' + (exponent / expscale);
1082               exponent %= expscale;
1083             }
1084           while (expscale > 10);
1085         *wcp++ = L'0' + exponent;
1086       }
1087
1088     /* Compute number of characters which must be filled with the padding
1089        character.  */
1090     if (is_neg || info->showsign || info->space)
1091       --width;
1092     width -= wcp - wstartp;
1093
1094     if (!info->left && info->pad != '0' && width > 0)
1095       PADN (info->pad, width);
1096
1097     if (is_neg)
1098       outchar ('-');
1099     else if (info->showsign)
1100       outchar ('+');
1101     else if (info->space)
1102       outchar (' ');
1103
1104     if (!info->left && info->pad == '0' && width > 0)
1105       PADN ('0', width);
1106
1107     {
1108       char *buffer = NULL;
1109       char *cp = NULL;
1110       char *tmpptr;
1111
1112       if (! wide)
1113         {
1114           /* Create the single byte string.  */
1115           size_t decimal_len;
1116           size_t thousands_sep_len;
1117           wchar_t *copywc;
1118
1119           decimal_len = strlen (decimal);
1120
1121           if (thousands_sep == NULL)
1122             thousands_sep_len = 0;
1123           else
1124             thousands_sep_len = strlen (thousands_sep);
1125
1126           if (__builtin_expect (buffer_malloced, 0))
1127             {
1128               buffer = (char *) malloc (2 + chars_needed + decimal_len
1129                                         + ngroups * thousands_sep_len);
1130               if (buffer == NULL)
1131                 {
1132                   /* Signal an error to the caller.  */
1133                   free (wbuffer);
1134                   return -1;
1135                 }
1136             }
1137           else
1138             buffer = (char *) alloca (2 + chars_needed + decimal_len
1139                                       + ngroups * thousands_sep_len);
1140
1141           /* Now copy the wide character string.  Since the character
1142              (except for the decimal point and thousands separator) must
1143              be coming from the ASCII range we can esily convert the
1144              string without mapping tables.  */
1145           for (cp = buffer, copywc = wstartp; copywc < wcp; ++copywc)
1146             if (*copywc == decimalwc)
1147               cp = (char *) __mempcpy (cp, decimal, decimal_len);
1148             else if (*copywc == thousands_sepwc)
1149               cp = (char *) __mempcpy (cp, thousands_sep, thousands_sep_len);
1150             else
1151               *cp++ = (char) *copywc;
1152         }
1153
1154       tmpptr = buffer;
1155       if (__builtin_expect (info->i18n, 0))
1156         {
1157 #ifdef COMPILE_WPRINTF
1158           wstartp = _i18n_number_rewrite (wstartp, wcp);
1159 #else
1160           tmpptr = _i18n_number_rewrite (tmpptr, cp);
1161 #endif
1162         }
1163
1164       PRINT (tmpptr, wstartp, wide ? wcp - wstartp : cp - tmpptr);
1165
1166       /* Free the memory if necessary.  */
1167       if (__builtin_expect (buffer_malloced, 0))
1168         {
1169           free (buffer);
1170           free (wbuffer);
1171         }
1172     }
1173
1174     if (info->left && width > 0)
1175       PADN (info->pad, width);
1176   }
1177   return done;
1178 }
1179 ldbl_hidden_def (___printf_fp, __printf_fp)
1180 ldbl_strong_alias (___printf_fp, __printf_fp)
1181 \f
1182 /* Return the number of extra grouping characters that will be inserted
1183    into a number with INTDIG_MAX integer digits.  */
1184
1185 unsigned int
1186 __guess_grouping (unsigned int intdig_max, const char *grouping)
1187 {
1188   unsigned int groups;
1189
1190   /* We treat all negative values like CHAR_MAX.  */
1191
1192   if (*grouping == CHAR_MAX || *grouping <= 0)
1193     /* No grouping should be done.  */
1194     return 0;
1195
1196   groups = 0;
1197   while (intdig_max > (unsigned int) *grouping)
1198     {
1199       ++groups;
1200       intdig_max -= *grouping++;
1201
1202       if (*grouping == CHAR_MAX
1203 #if CHAR_MIN < 0
1204           || *grouping < 0
1205 #endif
1206           )
1207         /* No more grouping should be done.  */
1208         break;
1209       else if (*grouping == 0)
1210         {
1211           /* Same grouping repeats.  */
1212           groups += (intdig_max - 1) / grouping[-1];
1213           break;
1214         }
1215     }
1216
1217   return groups;
1218 }
1219
1220 /* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
1221    There is guaranteed enough space past BUFEND to extend it.
1222    Return the new end of buffer.  */
1223
1224 static wchar_t *
1225 internal_function
1226 group_number (wchar_t *buf, wchar_t *bufend, unsigned int intdig_no,
1227               const char *grouping, wchar_t thousands_sep, int ngroups)
1228 {
1229   wchar_t *p;
1230
1231   if (ngroups == 0)
1232     return bufend;
1233
1234   /* Move the fractional part down.  */
1235   __wmemmove (buf + intdig_no + ngroups, buf + intdig_no,
1236               bufend - (buf + intdig_no));
1237
1238   p = buf + intdig_no + ngroups - 1;
1239   do
1240     {
1241       unsigned int len = *grouping++;
1242       do
1243         *p-- = buf[--intdig_no];
1244       while (--len > 0);
1245       *p-- = thousands_sep;
1246
1247       if (*grouping == CHAR_MAX
1248 #if CHAR_MIN < 0
1249           || *grouping < 0
1250 #endif
1251           )
1252         /* No more grouping should be done.  */
1253         break;
1254       else if (*grouping == 0)
1255         /* Same grouping repeats.  */
1256         --grouping;
1257     } while (intdig_no > (unsigned int) *grouping);
1258
1259   /* Copy the remaining ungrouped digits.  */
1260   do
1261     *p-- = buf[--intdig_no];
1262   while (p > buf);
1263
1264   return bufend + ngroups;
1265 }