Sat Nov 18 16:46:01 1995 Ulrich Drepper <drepper@gnu.ai.mit.edu>
[kopensolaris-gnu/glibc.git] / stdio-common / printf_fp.c
1 /* Floating point output for `printf'.
2 Copyright (C) 1995 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 (fp, info, args)
132      FILE *fp;
133      const struct printf_info *info;
134      const **const args;
135 {
136   /* The floating-point value to output.  */
137   union
138     {
139       double dbl;
140       LONG_DOUBLE ldbl;
141     }
142   fpnum;
143
144   /* Locale-dependent representation of decimal point.  */
145   wchar_t decimal;
146
147   /* Locale-dependent thousands separator and grouping specification.  */
148   wchar_t thousands_sep;
149   const char *grouping;
150
151   /* "NaN" or "Inf" for the special cases.  */
152   CONST char *special = NULL;
153
154   /* We need just a few limbs for the input before shifting to the right
155      position.  */
156   mp_limb fp_input[(LDBL_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB];
157   /* We need to shift the contents of fp_input by this amount of bits.  */
158   int to_shift;
159
160   /* The significant of the floting-point value in question  */
161   MPN_VAR(frac);
162   /* and the exponent.  */
163   int exponent;
164   /* Sign of the exponent.  */
165   int expsign = 0;
166   /* Sign of float number.  */
167   int is_neg = 0;
168
169   /* Scaling factor.  */
170   MPN_VAR(scale);
171
172   /* Temporary bignum value.  */
173   MPN_VAR(tmp);
174
175   /* Digit which is result of last hack_digit() call.  */
176   int digit;
177
178   /* The type of output format that will be used: 'e'/'E' or 'f'.  */
179   int type;
180
181   /* Counter for number of written characters.  */
182   int done = 0;
183
184   /* General helper (carry limb).  */
185   mp_limb cy;
186
187   char hack_digit (void)
188     {
189       mp_limb hi;
190
191       if (expsign != 0 && type == 'f' && exponent-- > 0)
192         hi = 0;
193       else if (scalesize == 0)
194         {
195           hi = frac[fracsize - 1];
196           cy = __mpn_mul_1 (frac, frac, fracsize - 1, 10);
197           frac[fracsize - 1] = cy;
198         }
199       else
200         {
201           if (fracsize < scalesize)
202             hi = 0;
203           else
204             {
205               hi = __mpn_divmod (tmp, frac, fracsize, scale, scalesize);
206               tmp[fracsize - scalesize] = hi;
207               hi = tmp[0];
208
209               fracsize = __mpn_normal_size (frac, scalesize);
210               if (fracsize == 0)
211                 {
212                   /* We're not prepared for an mpn variable with zero
213                      limbs.  */
214                   fracsize = 1;
215                   return '0' + hi;
216                 }
217             }
218
219           cy = __mpn_mul_1 (frac, frac, fracsize, 10);
220           if (cy != 0)
221             frac[fracsize++] = cy;
222         }
223
224       return '0' + hi;
225     }
226
227
228   /* Figure out the decimal point character.  */
229   if (mbtowc (&decimal, _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT),
230               strlen (_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT))) <= 0)
231     decimal = (wchar_t) *_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
232
233
234   if (info->group)
235     {
236       grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
237       if (*grouping <= 0 || *grouping == CHAR_MAX)
238         grouping = NULL;
239       else
240         {
241           /* Figure out the thousands seperator character.  */
242           if (mbtowc (&thousands_sep, _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP),
243                       strlen (_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP))) <= 0)
244             thousands_sep = (wchar_t) *_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
245           if (thousands_sep == L'\0')
246             grouping = NULL;
247         }
248     }
249   else
250     grouping = NULL;
251
252   /* Fetch the argument value.  */
253   if (info->is_long_double && sizeof (long double) > sizeof (double))
254     {
255       fpnum.ldbl = *(const long double *) args[0];
256
257       /* Check for special values: not a number or infinity.  */
258       if (__isnanl (fpnum.ldbl))
259         {
260           special = "NaN";
261           is_neg = 0;
262         }
263       else if (__isinfl (fpnum.ldbl))
264         {
265           special = "Inf";
266           is_neg = fpnum.ldbl < 0;
267         }
268       else
269         {
270           fracsize = __mpn_extract_long_double (fp_input,
271                                                 (sizeof (fp_input) /
272                                                  sizeof (fp_input[0])),
273                                                 &exponent, &is_neg,
274                                                 fpnum.ldbl);
275           to_shift = 1 + fracsize * BITS_PER_MP_LIMB - LDBL_MANT_DIG;
276         }
277     }
278   else
279     {
280       fpnum.dbl = *(const double *) args[0];
281
282       /* Check for special values: not a number or infinity.  */
283       if (__isnan (fpnum.dbl))
284         {
285           special = "NaN";
286           is_neg = 0;
287         }
288       else if (__isinf (fpnum.dbl))
289         {
290           special = "Inf";
291           is_neg = fpnum.dbl < 0;
292         }
293       else
294         {
295           fracsize = __mpn_extract_double (fp_input,
296                                            (sizeof (fp_input)
297                                             / sizeof (fp_input[0])),
298                                            &exponent, &is_neg, fpnum.dbl);
299           to_shift = 1 + fracsize * BITS_PER_MP_LIMB - DBL_MANT_DIG;
300         }
301     }
302
303   if (special)
304     {
305       int width = info->prec > info->width ? info->prec : info->width;
306
307       if (is_neg || info->showsign || info->space)
308         --width;
309       width -= 3;
310
311       if (!info->left && width > 0)
312         PADN (' ', width);
313
314       if (is_neg)
315         outchar ('-');
316       else if (info->showsign)
317         outchar ('+');
318       else if (info->space)
319         outchar (' ');
320
321       PRINT (special, 3);
322
323       if (info->left && width > 0)
324         PADN (' ', width);
325
326       return done;
327     }
328
329
330   /* We need three multiprecision variables.  Now that we have the exponent
331      of the number we can allocate the needed memory.  It would be more
332      efficient to use variables of the fixed maximum size but because this
333      would be really big it could lead to memory problems.  */
334   {
335     mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1)
336                              / BITS_PER_MP_LIMB + 3) * sizeof (mp_limb);
337     frac = (mp_limb *) alloca (bignum_size);
338     tmp = (mp_limb *) alloca (bignum_size);
339     scale = (mp_limb *) alloca (bignum_size);
340   }
341
342   /* We now have to distinguish between numbers with positive and negative
343      exponents because the method used for the one is not applicable/efficient
344      for the other.  */
345   scalesize = 0;
346   if (exponent > 2)
347     {
348       /* |FP| >= 1.0.  */
349       int scaleexpo = 0;
350       int explog = LDBL_MAX_10_EXP_LOG;
351       int exp10 = 0;
352       const struct mp_power *tens = &_fpioconst_pow10[explog + 1];
353       int cnt_h, cnt_l, i;
354
355       if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0)
356         {
357           MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
358                          fp_input, fracsize);
359           fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
360         }
361       else
362         {
363           cy = __mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
364                              fp_input, fracsize,
365                              (exponent + to_shift) % BITS_PER_MP_LIMB);
366           fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
367           if (cy)
368             frac[fracsize++] = cy;
369         }
370       MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB);
371
372       assert (tens > &_fpioconst_pow10[0]);
373       do
374         {
375           --tens;
376
377           /* The number of the product of two binary numbers with n and m
378              bits respectively has m+n or m+n-1 bits.   */
379           if (exponent >= scaleexpo + tens->p_expo - 1)
380             {
381               if (scalesize == 0)
382                 MPN_ASSIGN (tmp, tens->array);
383               else
384                 {
385                   cy = __mpn_mul (tmp, scale, scalesize,
386                                   tens->array + 2, tens->arraysize - 2);
387                   tmpsize = scalesize + tens->arraysize - 2;
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 - 2)
518                 cy = __mpn_mul (tmp, &tens->array[2], tens->arraysize - 2,
519                                 frac, fracsize);
520               else
521                 cy = __mpn_mul (tmp, frac, fracsize,
522                                 &tens->array[2], tens->arraysize - 2);
523               tmpsize = fracsize + tens->arraysize - 2;
524               if (cy == 0)
525                 --tmpsize;
526
527               count_leading_zeros (cnt_h, tmp[tmpsize - 1]);
528               incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB
529                      + BITS_PER_MP_LIMB - 1 - cnt_h;
530
531               assert (incr <= tens->p_expo);
532
533               /* If we increased the exponent by exactly 3 we have to test
534                  for overflow.  This is done by comparing with 10 shifted
535                  to the right position.  */
536               if (incr == exponent + 3)
537                 if (cnt_h <= BITS_PER_MP_LIMB - 4)
538                   {
539                     topval[0] = 0;
540                     topval[1] = 10 << (BITS_PER_MP_LIMB - 4 - cnt_h);
541                   }
542                 else
543                   {
544                     topval[0] = 10 << (BITS_PER_MP_LIMB - 4);
545                     topval[1] = 0;
546                     (void) __mpn_lshift (topval, topval, 2,
547                                          BITS_PER_MP_LIMB - cnt_h);
548                   }
549
550               /* We have to be careful when multiplying the last factor.
551                  If the result is greater than 1.0 be have to test it
552                  against 10.0.  If it is greater or equal to 10.0 the
553                  multiplication was not valid.  This is because we cannot
554                  determine the number of bits in the result in advance.  */
555               if (incr < exponent + 3
556                   || (incr == exponent + 3 &&
557                       (tmp[tmpsize - 1] < topval[1]
558                        || (tmp[tmpsize - 1] == topval[1]
559                            && tmp[tmpsize - 2] < topval[0]))))
560                 {
561                   /* The factor is right.  Adapt binary and decimal
562                      exponents.  */
563                   exponent -= incr;
564                   exp10 |= 1 << explog;
565
566                   /* If this factor yields a number greater or equal to
567                      1.0, we must not shift the non-fractional digits down. */
568                   if (exponent < 0)
569                     cnt_h += -exponent;
570
571                   /* Now we optimize the number representation.  */
572                   for (i = 0; tmp[i] == 0; ++i);
573                   if (cnt_h == BITS_PER_MP_LIMB - 1)
574                     {
575                       MPN_COPY (frac, tmp + i, tmpsize - i);
576                       fracsize = tmpsize - i;
577                     }
578                   else
579                     {
580                       count_trailing_zeros (cnt_l, tmp[i]);
581
582                       /* Now shift the numbers to their optimal position.  */
583                       if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l)
584                         {
585                           /* We cannot save any memory.  Just roll the
586                              number so that the leading digit is in a
587                              seperate limb.  */
588
589                           cy = __mpn_lshift (frac, tmp, tmpsize, cnt_h + 1);
590                           fracsize = tmpsize + 1;
591                           frac[fracsize - 1] = cy;
592                         }
593                       else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l)
594                         {
595                           (void) __mpn_rshift (frac, tmp + i, tmpsize - i,
596                                                BITS_PER_MP_LIMB - 1 - cnt_h);
597                           fracsize = tmpsize - i;
598                         }
599                       else
600                         {
601                           /* We can only save the memory of the limbs which
602                              are zero.  The non-zero parts occupy the same
603                              number of limbs.  */
604
605                           (void) __mpn_rshift (frac, tmp + (i - 1),
606                                                tmpsize - (i - 1),
607                                                BITS_PER_MP_LIMB - 1 - cnt_h);
608                           fracsize = tmpsize - (i - 1);
609                         }
610                     }
611                   used_limbs = fracsize - 1;
612                 }
613             }
614           --explog;
615         }
616       while (tens != &_fpioconst_pow10[1] && exponent > 0);
617       /* All factors but 10^-1 are tested now.  */
618       if (exponent > 0)
619         {
620           cy = __mpn_mul_1 (tmp, frac, fracsize, 10);
621           tmpsize = fracsize;
622           assert (cy == 0 || tmp[tmpsize - 1] < 20);
623
624           (void) __mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent));
625           fracsize = tmpsize;
626           exp10 |= 1;
627           assert (frac[fracsize - 1] < 10);
628         }
629       exponent = exp10;
630     }
631   else
632     {
633       /* This is a special case.  We don't need a factor because the
634          numbers are in the range of 0.0 <= fp < 8.0.  We simply
635          shift it to the right place and divide it by 1.0 to get the
636          leading digit.  (Of course this division is not really made.)  */
637       assert (0 <= exponent && exponent < 3 &&
638               exponent + to_shift < BITS_PER_MP_LIMB);
639
640       /* Now shift the input value to its right place.  */
641       cy = __mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift));
642       frac[fracsize++] = cy;
643       exponent = 0;
644     }
645
646   {
647     int width = info->width;
648     char *buffer, *startp, *cp;
649     int chars_needed;
650     int expscale;
651     int intdig_max, intdig_no = 0;
652     int fracdig_min, fracdig_max, fracdig_no = 0;
653     int dig_max;
654     int significant;
655
656     if (tolower (info->spec) == 'e')
657       {
658         type = info->spec;
659         intdig_max = 1;
660         fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
661         chars_needed = 1 + 1 + fracdig_max + 1 + 1 + 4;
662         /*             d   .     ddd         e   +-  ddd  */
663         dig_max = INT_MAX;              /* Unlimited.  */
664         significant = 1;                /* Does not matter here.  */
665       }
666     else if (info->spec == 'f')
667       {
668         type = 'f';
669         fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
670         if (expsign == 0)
671           {
672             intdig_max = exponent + 1;
673             /* This can be really big!  */  /* XXX Maybe malloc if too big? */
674             chars_needed = exponent + 1 + 1 + fracdig_max;
675           }
676         else
677           {
678             intdig_max = 1;
679             chars_needed = 1 + 1 + fracdig_max;
680           }
681         dig_max = INT_MAX;              /* Unlimited.  */
682         significant = 1;                /* Does not matter here.  */
683       }
684     else
685       {
686         dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec);
687         if ((expsign == 0 && exponent >= dig_max)
688             || (expsign != 0 && exponent > 4))
689           {
690             type = isupper (info->spec) ? 'E' : 'e';
691             fracdig_max = dig_max - 1;
692             intdig_max = 1;
693             chars_needed = 1 + 1 + fracdig_max + 1 + 1 + 4;
694           }
695         else
696           {
697             type = 'f';
698             intdig_max = expsign == 0 ? exponent + 1 : 0;
699             fracdig_max = dig_max - intdig_max;
700             /* We need space for the significant digits and perhaps for
701                leading zeros when < 1.0.  Pessimistic guess: dig_max.  */
702             chars_needed = dig_max + dig_max + 1;
703           }
704         fracdig_min = info->alt ? fracdig_max : 0;
705         significant = 0;                /* We count significant digits.  */
706       }
707
708     if (grouping)
709       /* Guess the number of groups we will make, and thus how
710          many spaces we need for separator characters.  */
711       chars_needed += guess_grouping (intdig_max, grouping, thousands_sep);
712
713     /* Allocate buffer for output.  We need two more because while rounding
714        it is possible that we need two more characters in front of all the
715        other output.  */
716     buffer = alloca (2 + chars_needed);
717     cp = startp = buffer + 2;   /* Let room for rounding.  */
718
719     /* Do the real work: put digits in allocated buffer.  */
720     if (expsign == 0 || type != 'f')
721       {
722         assert (expsign == 0 || intdig_max == 1);
723         while (intdig_no < intdig_max)
724           {
725             ++intdig_no;
726             *cp++ = hack_digit ();
727           }
728         significant = 1;
729         if (info->alt
730             || fracdig_min > 0
731             || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0)))
732           *cp++ = decimal;
733       }
734     else
735       {
736         /* |fp| < 1.0 and the selected type is 'f', so put "0."
737            in the buffer.  */
738         *cp++ = '0';
739         --exponent;
740         *cp++ = decimal;
741       }
742
743     /* Generate the needed number of fractional digits.  */
744     while (fracdig_no < fracdig_min
745            || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
746       {
747         ++fracdig_no;
748         *cp = hack_digit ();
749         if (*cp != '0')
750           significant = 1;
751         else if (significant == 0)
752           {
753             ++fracdig_max;
754             if (fracdig_min > 0)
755               ++fracdig_min;
756           }
757         ++cp;
758       }
759
760     /* Do rounding.  */
761     digit = hack_digit ();
762     if (digit > '4')
763       {
764         char *tp = cp;
765
766         if (digit == '5')
767           /* This is the critical case.  */
768           if (fracsize == 1 && frac[0] == 0)
769             /* Rest of the number is zero -> round to even.
770                (IEEE 754-1985 4.1 says this is the default rounding.)  */
771             if ((*(cp - 1) & 1) == 0)
772               goto do_expo;
773
774         if (fracdig_no > 0)
775           {
776             /* Process fractional digits.  Terminate if not rounded or
777                radix character is reached.  */
778             while (*--tp != decimal && *tp == '9')
779               *tp = '0';
780             if (*tp != decimal)
781               /* Round up.  */
782               (*tp)++;
783           }
784
785         if (fracdig_no == 0 || *tp == decimal)
786           {
787             /* Round the integer digits.  */
788             if (*(tp - 1) == decimal)
789               --tp;
790
791             while (--tp >= startp && *tp == '9')
792               *tp = '0';
793
794             if (tp >= startp)
795               /* Round up.  */
796               (*tp)++;
797             else
798               /* It is more citical.  All digits were 9's.  */
799               {
800                 if (type != 'f')
801                   {
802                     *startp = '1';
803                     exponent += expsign == 0 ? 1 : -1;
804                   }
805                 else if (intdig_no == dig_max)
806                   {
807                     /* This is the case where for type %g the number fits
808                        really in the range for %f output but after rounding
809                        the number of digits is too big.  */
810                     *--startp = decimal;
811                     *--startp = '1';
812
813                     if (info->alt || fracdig_no > 0)
814                       {
815                         /* Overwrite the old radix character.  */
816                         startp[intdig_no + 2] = '0';
817                         ++fracdig_no;
818                       }
819
820                     fracdig_no += intdig_no;
821                     intdig_no = 1;
822                     fracdig_max = intdig_max - intdig_no;
823                     ++exponent;
824                     /* Now we must print the exponent.  */
825                     type = isupper (info->spec) ? 'E' : 'e';
826                   }
827                 else
828                   {
829                     /* We can simply add another another digit before the
830                        radix.  */
831                     *--startp = '1';
832                     ++intdig_no;
833                   }
834
835                 /* While rounding the number of digits can change.
836                    If the number now exceeds the limits remove some
837                    fractional digits.  */
838                 if (intdig_no + fracdig_no > dig_max)
839                   {
840                     cp -= intdig_no + fracdig_no - dig_max;
841                     fracdig_no -= intdig_no + fracdig_no - dig_max;
842                   }
843               }
844           }
845       }
846
847   do_expo:
848     /* Now remove unnecessary '0' at the end of the string.  */
849     while (fracdig_no > fracdig_min && *(cp - 1) == '0')
850       {
851         --cp;
852         --fracdig_no;
853       }
854     /* If we eliminate all fractional digits we perhaps also can remove
855        the radix character.  */
856     if (fracdig_no == 0 && !info->alt && *(cp - 1) == decimal)
857       --cp;
858
859     if (grouping)
860       /* Add in separator characters, overwriting the same buffer.  */
861       cp = group_number (startp, cp, intdig_no, grouping, thousands_sep);
862
863     /* Write the exponent if it is needed.  */
864     if (type != 'f')
865       {
866         *cp++ = type;
867         *cp++ = expsign ? '-' : '+';
868
869         /* Find the magnitude of the exponent.  */
870         expscale = 10;
871         while (expscale <= exponent)
872           expscale *= 10;
873
874         if (exponent < 10)
875           /* Exponent always has at least two digits.  */
876           *cp++ = '0';
877         else
878           do
879             {
880               expscale /= 10;
881               *cp++ = '0' + (exponent / expscale);
882               exponent %= expscale;
883             }
884           while (expscale > 10);
885         *cp++ = '0' + exponent;
886       }
887
888     /* Compute number of characters which must be filled with the padding
889        character.  */
890     if (is_neg || info->showsign || info->space)
891       --width;
892     width -= cp - startp;
893
894     if (!info->left && info->pad != '0' && width > 0)
895       PADN (info->pad, width);
896
897     if (is_neg)
898       outchar ('-');
899     else if (info->showsign)
900       outchar ('+');
901     else if (info->space)
902       outchar (' ');
903
904     if (!info->left && info->pad == '0' && width > 0)
905       PADN ('0', width);
906
907     PRINT (startp, cp - startp);
908
909     if (info->left && width > 0)
910       PADN (info->pad, width);
911   }
912   return done;
913 }
914 \f
915 /* Return the number of extra grouping characters that will be inserted
916    into a number with INTDIG_MAX integer digits.  */
917
918 static unsigned int
919 guess_grouping (unsigned int intdig_max, const char *grouping, wchar_t sepchar)
920 {
921   unsigned int groups;
922
923   /* We treat all negative values like CHAR_MAX.  */
924
925   if (*grouping == CHAR_MAX || *grouping <= 0)
926     /* No grouping should be done.  */
927     return 0;
928
929   groups = 0;
930   while (intdig_max > (unsigned int) *grouping)
931     {
932       ++groups;
933       intdig_max -= *grouping++;
934
935       if (*grouping == CHAR_MAX || *grouping < 0)
936         /* No more grouping should be done.  */
937         break;
938       else if (*grouping == 0)
939         {
940           /* Same grouping repeats.  */
941           groups += intdig_max / grouping[-1];
942           break;
943         }
944     }
945
946   return groups;
947 }
948
949 /* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
950    There is guaranteed enough space past BUFEND to extend it.
951    Return the new end of buffer.  */
952
953 static char *
954 group_number (char *buf, char *bufend, unsigned int intdig_no,
955               const char *grouping, wchar_t thousands_sep)
956 {
957   unsigned int groups = guess_grouping (intdig_no, grouping, thousands_sep);
958   char *p;
959
960   if (groups == 0)
961     return bufend;
962
963   /* Move the fractional part down.  */
964   memmove (buf + intdig_no + groups, buf + intdig_no,
965            bufend - (buf + intdig_no));
966
967   p = buf + intdig_no + groups - 1;
968   do
969     {
970       unsigned int len = *grouping++;
971       do
972         *p-- = buf[--intdig_no];
973       while (--len > 0);
974       *p-- = thousands_sep;
975
976       if (*grouping == CHAR_MAX || *grouping < 0)
977         /* No more grouping should be done.  */
978         break;
979       else if (*grouping == 0)
980         /* Same grouping repeats.  */
981         --grouping;
982     } while (intdig_no > (unsigned int) *grouping);
983
984   /* Copy the remaining ungrouped digits.  */
985   do
986     *p-- = buf[--intdig_no];
987   while (p > buf);
988
989   return bufend + groups;
990 }