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