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