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