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