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