update from main archive 961030
[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 extern unsigned int __guess_grouping (unsigned int intdig_max,
125                                       const char *grouping, wchar_t sepchar);
126
127
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 = 0;
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 (info->extra == 0)
233     {
234       if (mbtowc (&decimal, _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT),
235                   strlen (_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT))) <= 0)
236         decimal = (wchar_t) *_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
237     }
238   else
239     {
240       if (mbtowc (&decimal, _NL_CURRENT (LC_MONETARY, MON_DECIMAL_POINT),
241                   strlen (_NL_CURRENT (LC_MONETARY, MON_DECIMAL_POINT))) <= 0)
242         decimal = (wchar_t) *_NL_CURRENT (LC_MONETARY, MON_DECIMAL_POINT);
243     }
244   /* Give default value.  */
245   if (decimal == L'\0')
246     decimal = L'.';
247
248
249   if (info->group)
250     {
251       if (info->extra == 0)
252         grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
253       else
254         grouping = _NL_CURRENT (LC_MONETARY, MON_GROUPING);
255
256       if (*grouping <= 0 || *grouping == CHAR_MAX)
257         grouping = NULL;
258       else
259         {
260           /* Figure out the thousands seperator character.  */
261           if (info->extra == 0)
262             {
263               if (mbtowc (&thousands_sep, _NL_CURRENT (LC_NUMERIC,
264                                                        THOUSANDS_SEP),
265                           strlen (_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP)))
266                   <= 0)
267                 thousands_sep = (wchar_t) *_NL_CURRENT (LC_NUMERIC,
268                                                         THOUSANDS_SEP);
269             }
270           else
271             {
272               if (mbtowc (&thousands_sep, _NL_CURRENT (LC_MONETARY,
273                                                        MON_THOUSANDS_SEP),
274                           strlen (_NL_CURRENT (LC_MONETARY,
275                                                MON_THOUSANDS_SEP))) <= 0)
276                 thousands_sep = (wchar_t) *_NL_CURRENT (LC_MONETARY,
277                                                         MON_THOUSANDS_SEP);
278             }
279
280           if (thousands_sep == L'\0')
281             grouping = NULL;
282         }
283     }
284   else
285     grouping = NULL;
286
287   /* Fetch the argument value.  */
288   if (info->is_long_double && sizeof (long double) > sizeof (double))
289     {
290       fpnum.ldbl = *(const long double *) args[0];
291
292       /* Check for special values: not a number or infinity.  */
293       if (__isnanl (fpnum.ldbl))
294         {
295           special = "NaN";
296           is_neg = 0;
297         }
298       else if (__isinfl (fpnum.ldbl))
299         {
300           special = "Inf";
301           is_neg = fpnum.ldbl < 0;
302         }
303       else
304         {
305           fracsize = __mpn_extract_long_double (fp_input,
306                                                 (sizeof (fp_input) /
307                                                  sizeof (fp_input[0])),
308                                                 &exponent, &is_neg,
309                                                 fpnum.ldbl);
310           to_shift = 1 + fracsize * BITS_PER_MP_LIMB - LDBL_MANT_DIG;
311         }
312     }
313   else
314     {
315       fpnum.dbl = *(const double *) args[0];
316
317       /* Check for special values: not a number or infinity.  */
318       if (__isnan (fpnum.dbl))
319         {
320           special = "NaN";
321           is_neg = 0;
322         }
323       else if (__isinf (fpnum.dbl))
324         {
325           special = "Inf";
326           is_neg = fpnum.dbl < 0;
327         }
328       else
329         {
330           fracsize = __mpn_extract_double (fp_input,
331                                            (sizeof (fp_input)
332                                             / sizeof (fp_input[0])),
333                                            &exponent, &is_neg, fpnum.dbl);
334           to_shift = 1 + fracsize * BITS_PER_MP_LIMB - DBL_MANT_DIG;
335         }
336     }
337
338   if (special)
339     {
340       int width = info->prec > info->width ? info->prec : info->width;
341
342       if (is_neg || info->showsign || info->space)
343         --width;
344       width -= 3;
345
346       if (!info->left && width > 0)
347         PADN (' ', width);
348
349       if (is_neg)
350         outchar ('-');
351       else if (info->showsign)
352         outchar ('+');
353       else if (info->space)
354         outchar (' ');
355
356       PRINT (special, 3);
357
358       if (info->left && width > 0)
359         PADN (' ', width);
360
361       return done;
362     }
363
364
365   /* We need three multiprecision variables.  Now that we have the exponent
366      of the number we can allocate the needed memory.  It would be more
367      efficient to use variables of the fixed maximum size but because this
368      would be really big it could lead to memory problems.  */
369   {
370     mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1)
371                              / BITS_PER_MP_LIMB + 4) * sizeof (mp_limb_t);
372     frac = (mp_limb_t *) alloca (bignum_size);
373     tmp = (mp_limb_t *) alloca (bignum_size);
374     scale = (mp_limb_t *) alloca (bignum_size);
375   }
376
377   /* We now have to distinguish between numbers with positive and negative
378      exponents because the method used for the one is not applicable/efficient
379      for the other.  */
380   scalesize = 0;
381   if (exponent > 2)
382     {
383       /* |FP| >= 8.0.  */
384       int scaleexpo = 0;
385       int explog = LDBL_MAX_10_EXP_LOG;
386       int exp10 = 0;
387       const struct mp_power *tens = &_fpioconst_pow10[explog + 1];
388       int cnt_h, cnt_l, i;
389
390       if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0)
391         {
392           MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
393                          fp_input, fracsize);
394           fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
395         }
396       else
397         {
398           cy = __mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
399                              fp_input, fracsize,
400                              (exponent + to_shift) % BITS_PER_MP_LIMB);
401           fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
402           if (cy)
403             frac[fracsize++] = cy;
404         }
405       MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB);
406
407       assert (tens > &_fpioconst_pow10[0]);
408       do
409         {
410           --tens;
411
412           /* The number of the product of two binary numbers with n and m
413              bits respectively has m+n or m+n-1 bits.   */
414           if (exponent >= scaleexpo + tens->p_expo - 1)
415             {
416               if (scalesize == 0)
417                 MPN_ASSIGN (tmp, tens->array);
418               else
419                 {
420                   cy = __mpn_mul (tmp, scale, scalesize,
421                                   &tens->array[_FPIO_CONST_OFFSET],
422                                   tens->arraysize - _FPIO_CONST_OFFSET);
423                   tmpsize = scalesize + tens->arraysize - _FPIO_CONST_OFFSET;
424                   if (cy == 0)
425                     --tmpsize;
426                 }
427
428               if (MPN_GE (frac, tmp))
429                 {
430                   int cnt;
431                   MPN_ASSIGN (scale, tmp);
432                   count_leading_zeros (cnt, scale[scalesize - 1]);
433                   scaleexpo = (scalesize - 2) * BITS_PER_MP_LIMB - cnt - 1;
434                   exp10 |= 1 << explog;
435                 }
436             }
437           --explog;
438         }
439       while (tens > &_fpioconst_pow10[0]);
440       exponent = exp10;
441
442       /* Optimize number representations.  We want to represent the numbers
443          with the lowest number of bytes possible without losing any
444          bytes. Also the highest bit in the scaling factor has to be set
445          (this is a requirement of the MPN division routines).  */
446       if (scalesize > 0)
447         {
448           /* Determine minimum number of zero bits at the end of
449              both numbers.  */
450           for (i = 0; scale[i] == 0 && frac[i] == 0; i++)
451             ;
452
453           /* Determine number of bits the scaling factor is misplaced.  */
454           count_leading_zeros (cnt_h, scale[scalesize - 1]);
455
456           if (cnt_h == 0)
457             {
458               /* The highest bit of the scaling factor is already set.  So
459                  we only have to remove the trailing empty limbs.  */
460               if (i > 0)
461                 {
462                   MPN_COPY_INCR (scale, scale + i, scalesize - i);
463                   scalesize -= i;
464                   MPN_COPY_INCR (frac, frac + i, fracsize - i);
465                   fracsize -= i;
466                 }
467             }
468           else
469             {
470               if (scale[i] != 0)
471                 {
472                   count_trailing_zeros (cnt_l, scale[i]);
473                   if (frac[i] != 0)
474                     {
475                       int cnt_l2;
476                       count_trailing_zeros (cnt_l2, frac[i]);
477                       if (cnt_l2 < cnt_l)
478                         cnt_l = cnt_l2;
479                     }
480                 }
481               else
482                 count_trailing_zeros (cnt_l, frac[i]);
483
484               /* Now shift the numbers to their optimal position.  */
485               if (i == 0 && BITS_PER_MP_LIMB - cnt_h > cnt_l)
486                 {
487                   /* We cannot save any memory.  So just roll both numbers
488                      so that the scaling factor has its highest bit set.  */
489
490                   (void) __mpn_lshift (scale, scale, scalesize, cnt_h);
491                   cy = __mpn_lshift (frac, frac, fracsize, cnt_h);
492                   if (cy != 0)
493                     frac[fracsize++] = cy;
494                 }
495               else if (BITS_PER_MP_LIMB - cnt_h <= cnt_l)
496                 {
497                   /* We can save memory by removing the trailing zero limbs
498                      and by packing the non-zero limbs which gain another
499                      free one. */
500
501                   (void) __mpn_rshift (scale, scale + i, scalesize - i,
502                                        BITS_PER_MP_LIMB - cnt_h);
503                   scalesize -= i + 1;
504                   (void) __mpn_rshift (frac, frac + i, fracsize - i,
505                                        BITS_PER_MP_LIMB - cnt_h);
506                   fracsize -= frac[fracsize - i - 1] == 0 ? i + 1 : i;
507                 }
508               else
509                 {
510                   /* We can only save the memory of the limbs which are zero.
511                      The non-zero parts occupy the same number of limbs.  */
512
513                   (void) __mpn_rshift (scale, scale + (i - 1),
514                                        scalesize - (i - 1),
515                                        BITS_PER_MP_LIMB - cnt_h);
516                   scalesize -= i;
517                   (void) __mpn_rshift (frac, frac + (i - 1),
518                                        fracsize - (i - 1),
519                                        BITS_PER_MP_LIMB - cnt_h);
520                   fracsize -= frac[fracsize - (i - 1) - 1] == 0 ? i : i - 1;
521                 }
522             }
523         }
524     }
525   else if (exponent < 0)
526     {
527       /* |FP| < 1.0.  */
528       int exp10 = 0;
529       int explog = LDBL_MAX_10_EXP_LOG;
530       const struct mp_power *tens = &_fpioconst_pow10[explog + 1];
531       mp_size_t used_limbs = fracsize - 1;
532
533       /* Now shift the input value to its right place.  */
534       cy = __mpn_lshift (frac, fp_input, fracsize, to_shift);
535       frac[fracsize++] = cy;
536       assert (cy == 1 || (frac[fracsize - 2] == 0 && frac[0] == 0));
537
538       expsign = 1;
539       exponent = -exponent;
540
541       assert (tens != &_fpioconst_pow10[0]);
542       do
543         {
544           --tens;
545
546           if (exponent >= tens->m_expo)
547             {
548               int i, incr, cnt_h, cnt_l;
549               mp_limb_t topval[2];
550
551               /* The __mpn_mul function expects the first argument to be
552                  bigger than the second.  */
553               if (fracsize < tens->arraysize - _FPIO_CONST_OFFSET)
554                 cy = __mpn_mul (tmp, &tens->array[_FPIO_CONST_OFFSET],
555                                 tens->arraysize - _FPIO_CONST_OFFSET,
556                                 frac, fracsize);
557               else
558                 cy = __mpn_mul (tmp, frac, fracsize,
559                                 &tens->array[_FPIO_CONST_OFFSET],
560                                 tens->arraysize - _FPIO_CONST_OFFSET);
561               tmpsize = fracsize + tens->arraysize - _FPIO_CONST_OFFSET;
562               if (cy == 0)
563                 --tmpsize;
564
565               count_leading_zeros (cnt_h, tmp[tmpsize - 1]);
566               incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB
567                      + BITS_PER_MP_LIMB - 1 - cnt_h;
568
569               assert (incr <= tens->p_expo);
570
571               /* If we increased the exponent by exactly 3 we have to test
572                  for overflow.  This is done by comparing with 10 shifted
573                  to the right position.  */
574               if (incr == exponent + 3)
575                 if (cnt_h <= BITS_PER_MP_LIMB - 4)
576                   {
577                     topval[0] = 0;
578                     topval[1]
579                       = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4 - cnt_h);
580                   }
581                 else
582                   {
583                     topval[0] = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4);
584                     topval[1] = 0;
585                     (void) __mpn_lshift (topval, topval, 2,
586                                          BITS_PER_MP_LIMB - cnt_h);
587                   }
588
589               /* We have to be careful when multiplying the last factor.
590                  If the result is greater than 1.0 be have to test it
591                  against 10.0.  If it is greater or equal to 10.0 the
592                  multiplication was not valid.  This is because we cannot
593                  determine the number of bits in the result in advance.  */
594               if (incr < exponent + 3
595                   || (incr == exponent + 3 &&
596                       (tmp[tmpsize - 1] < topval[1]
597                        || (tmp[tmpsize - 1] == topval[1]
598                            && tmp[tmpsize - 2] < topval[0]))))
599                 {
600                   /* The factor is right.  Adapt binary and decimal
601                      exponents.  */
602                   exponent -= incr;
603                   exp10 |= 1 << explog;
604
605                   /* If this factor yields a number greater or equal to
606                      1.0, we must not shift the non-fractional digits down. */
607                   if (exponent < 0)
608                     cnt_h += -exponent;
609
610                   /* Now we optimize the number representation.  */
611                   for (i = 0; tmp[i] == 0; ++i);
612                   if (cnt_h == BITS_PER_MP_LIMB - 1)
613                     {
614                       MPN_COPY (frac, tmp + i, tmpsize - i);
615                       fracsize = tmpsize - i;
616                     }
617                   else
618                     {
619                       count_trailing_zeros (cnt_l, tmp[i]);
620
621                       /* Now shift the numbers to their optimal position.  */
622                       if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l)
623                         {
624                           /* We cannot save any memory.  Just roll the
625                              number so that the leading digit is in a
626                              seperate limb.  */
627
628                           cy = __mpn_lshift (frac, tmp, tmpsize, cnt_h + 1);
629                           fracsize = tmpsize + 1;
630                           frac[fracsize - 1] = cy;
631                         }
632                       else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l)
633                         {
634                           (void) __mpn_rshift (frac, tmp + i, tmpsize - i,
635                                                BITS_PER_MP_LIMB - 1 - cnt_h);
636                           fracsize = tmpsize - i;
637                         }
638                       else
639                         {
640                           /* We can only save the memory of the limbs which
641                              are zero.  The non-zero parts occupy the same
642                              number of limbs.  */
643
644                           (void) __mpn_rshift (frac, tmp + (i - 1),
645                                                tmpsize - (i - 1),
646                                                BITS_PER_MP_LIMB - 1 - cnt_h);
647                           fracsize = tmpsize - (i - 1);
648                         }
649                     }
650                   used_limbs = fracsize - 1;
651                 }
652             }
653           --explog;
654         }
655       while (tens != &_fpioconst_pow10[1] && exponent > 0);
656       /* All factors but 10^-1 are tested now.  */
657       if (exponent > 0)
658         {
659           int cnt_l;
660
661           cy = __mpn_mul_1 (tmp, frac, fracsize, 10);
662           tmpsize = fracsize;
663           assert (cy == 0 || tmp[tmpsize - 1] < 20);
664
665           count_trailing_zeros (cnt_l, tmp[0]);
666           if (cnt_l < MIN (4, exponent))
667             {
668               cy = __mpn_lshift (frac, tmp, tmpsize,
669                                  BITS_PER_MP_LIMB - MIN (4, exponent));
670               if (cy != 0)
671                 frac[tmpsize++] = cy;
672             }
673           else
674             (void) __mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent));
675           fracsize = tmpsize;
676           exp10 |= 1;
677           assert (frac[fracsize - 1] < 10);
678         }
679       exponent = exp10;
680     }
681   else
682     {
683       /* This is a special case.  We don't need a factor because the
684          numbers are in the range of 0.0 <= fp < 8.0.  We simply
685          shift it to the right place and divide it by 1.0 to get the
686          leading digit.  (Of course this division is not really made.)  */
687       assert (0 <= exponent && exponent < 3 &&
688               exponent + to_shift < BITS_PER_MP_LIMB);
689
690       /* Now shift the input value to its right place.  */
691       cy = __mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift));
692       frac[fracsize++] = cy;
693       exponent = 0;
694     }
695
696   {
697     int width = info->width;
698     char *buffer, *startp, *cp;
699     int chars_needed;
700     int expscale;
701     int intdig_max, intdig_no = 0;
702     int fracdig_min, fracdig_max, fracdig_no = 0;
703     int dig_max;
704     int significant;
705
706     if (tolower (info->spec) == 'e')
707       {
708         type = info->spec;
709         intdig_max = 1;
710         fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
711         chars_needed = 1 + 1 + fracdig_max + 1 + 1 + 4;
712         /*             d   .     ddd         e   +-  ddd  */
713         dig_max = INT_MAX;              /* Unlimited.  */
714         significant = 1;                /* Does not matter here.  */
715       }
716     else if (info->spec == 'f')
717       {
718         type = 'f';
719         fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
720         if (expsign == 0)
721           {
722             intdig_max = exponent + 1;
723             /* This can be really big!  */  /* XXX Maybe malloc if too big? */
724             chars_needed = exponent + 1 + 1 + fracdig_max;
725           }
726         else
727           {
728             intdig_max = 1;
729             chars_needed = 1 + 1 + fracdig_max;
730           }
731         dig_max = INT_MAX;              /* Unlimited.  */
732         significant = 1;                /* Does not matter here.  */
733       }
734     else
735       {
736         dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec);
737         if ((expsign == 0 && exponent >= dig_max)
738             || (expsign != 0 && exponent > 4))
739           {
740             type = isupper (info->spec) ? 'E' : 'e';
741             fracdig_max = dig_max - 1;
742             intdig_max = 1;
743             chars_needed = 1 + 1 + fracdig_max + 1 + 1 + 4;
744           }
745         else
746           {
747             type = 'f';
748             intdig_max = expsign == 0 ? exponent + 1 : 0;
749             fracdig_max = dig_max - intdig_max;
750             /* We need space for the significant digits and perhaps for
751                leading zeros when < 1.0.  Pessimistic guess: dig_max.  */
752             chars_needed = dig_max + dig_max + 1;
753           }
754         fracdig_min = info->alt ? fracdig_max : 0;
755         significant = 0;                /* We count significant digits.  */
756       }
757
758     if (grouping)
759       /* Guess the number of groups we will make, and thus how
760          many spaces we need for separator characters.  */
761       chars_needed += __guess_grouping (intdig_max, grouping, thousands_sep);
762
763     /* Allocate buffer for output.  We need two more because while rounding
764        it is possible that we need two more characters in front of all the
765        other output.  */
766     buffer = alloca (2 + chars_needed);
767     cp = startp = buffer + 2;   /* Let room for rounding.  */
768
769     /* Do the real work: put digits in allocated buffer.  */
770     if (expsign == 0 || type != 'f')
771       {
772         assert (expsign == 0 || intdig_max == 1);
773         while (intdig_no < intdig_max)
774           {
775             ++intdig_no;
776             *cp++ = hack_digit ();
777           }
778         significant = 1;
779         if (info->alt
780             || fracdig_min > 0
781             || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0)))
782           *cp++ = decimal;
783       }
784     else
785       {
786         /* |fp| < 1.0 and the selected type is 'f', so put "0."
787            in the buffer.  */
788         *cp++ = '0';
789         --exponent;
790         *cp++ = decimal;
791       }
792
793     /* Generate the needed number of fractional digits.  */
794     while (fracdig_no < fracdig_min
795            || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
796       {
797         ++fracdig_no;
798         *cp = hack_digit ();
799         if (*cp != '0')
800           significant = 1;
801         else if (significant == 0)
802           {
803             ++fracdig_max;
804             if (fracdig_min > 0)
805               ++fracdig_min;
806           }
807         ++cp;
808       }
809
810     /* Do rounding.  */
811     digit = hack_digit ();
812     if (digit > '4')
813       {
814         char *tp = cp;
815
816         if (digit == '5')
817           /* This is the critical case.  */
818           if (fracsize == 1 && frac[0] == 0)
819             /* Rest of the number is zero -> round to even.
820                (IEEE 754-1985 4.1 says this is the default rounding.)  */
821             if ((*(cp - 1) & 1) == 0)
822               goto do_expo;
823
824         if (fracdig_no > 0)
825           {
826             /* Process fractional digits.  Terminate if not rounded or
827                radix character is reached.  */
828             while (*--tp != decimal && *tp == '9')
829               *tp = '0';
830             if (*tp != decimal)
831               /* Round up.  */
832               (*tp)++;
833           }
834
835         if (fracdig_no == 0 || *tp == decimal)
836           {
837             /* Round the integer digits.  */
838             if (*(tp - 1) == decimal)
839               --tp;
840
841             while (--tp >= startp && *tp == '9')
842               *tp = '0';
843
844             if (tp >= startp)
845               /* Round up.  */
846               (*tp)++;
847             else
848               /* It is more citical.  All digits were 9's.  */
849               {
850                 if (type != 'f')
851                   {
852                     *startp = '1';
853                     exponent += expsign == 0 ? 1 : -1;
854                   }
855                 else if (intdig_no == dig_max)
856                   {
857                     /* This is the case where for type %g the number fits
858                        really in the range for %f output but after rounding
859                        the number of digits is too big.  */
860                     *--startp = decimal;
861                     *--startp = '1';
862
863                     if (info->alt || fracdig_no > 0)
864                       {
865                         /* Overwrite the old radix character.  */
866                         startp[intdig_no + 2] = '0';
867                         ++fracdig_no;
868                       }
869
870                     fracdig_no += intdig_no;
871                     intdig_no = 1;
872                     fracdig_max = intdig_max - intdig_no;
873                     ++exponent;
874                     /* Now we must print the exponent.  */
875                     type = isupper (info->spec) ? 'E' : 'e';
876                   }
877                 else
878                   {
879                     /* We can simply add another another digit before the
880                        radix.  */
881                     *--startp = '1';
882                     ++intdig_no;
883                   }
884
885                 /* While rounding the number of digits can change.
886                    If the number now exceeds the limits remove some
887                    fractional digits.  */
888                 if (intdig_no + fracdig_no > dig_max)
889                   {
890                     cp -= intdig_no + fracdig_no - dig_max;
891                     fracdig_no -= intdig_no + fracdig_no - dig_max;
892                   }
893               }
894           }
895       }
896
897   do_expo:
898     /* Now remove unnecessary '0' at the end of the string.  */
899     while (fracdig_no > fracdig_min && *(cp - 1) == '0')
900       {
901         --cp;
902         --fracdig_no;
903       }
904     /* If we eliminate all fractional digits we perhaps also can remove
905        the radix character.  */
906     if (fracdig_no == 0 && !info->alt && *(cp - 1) == decimal)
907       --cp;
908
909     if (grouping)
910       /* Add in separator characters, overwriting the same buffer.  */
911       cp = group_number (startp, cp, intdig_no, grouping, thousands_sep);
912
913     /* Write the exponent if it is needed.  */
914     if (type != 'f')
915       {
916         *cp++ = type;
917         *cp++ = expsign ? '-' : '+';
918
919         /* Find the magnitude of the exponent.  */
920         expscale = 10;
921         while (expscale <= exponent)
922           expscale *= 10;
923
924         if (exponent < 10)
925           /* Exponent always has at least two digits.  */
926           *cp++ = '0';
927         else
928           do
929             {
930               expscale /= 10;
931               *cp++ = '0' + (exponent / expscale);
932               exponent %= expscale;
933             }
934           while (expscale > 10);
935         *cp++ = '0' + exponent;
936       }
937
938     /* Compute number of characters which must be filled with the padding
939        character.  */
940     if (is_neg || info->showsign || info->space)
941       --width;
942     width -= cp - startp;
943
944     if (!info->left && info->pad != '0' && width > 0)
945       PADN (info->pad, width);
946
947     if (is_neg)
948       outchar ('-');
949     else if (info->showsign)
950       outchar ('+');
951     else if (info->space)
952       outchar (' ');
953
954     if (!info->left && info->pad == '0' && width > 0)
955       PADN ('0', width);
956
957     PRINT (startp, cp - startp);
958
959     if (info->left && width > 0)
960       PADN (info->pad, width);
961   }
962   return done;
963 }
964 \f
965 /* Return the number of extra grouping characters that will be inserted
966    into a number with INTDIG_MAX integer digits.  */
967
968 unsigned int
969 __guess_grouping (unsigned int intdig_max, const char *grouping,
970                   wchar_t sepchar)
971 {
972   unsigned int groups;
973
974   /* We treat all negative values like CHAR_MAX.  */
975
976   if (*grouping == CHAR_MAX || *grouping <= 0)
977     /* No grouping should be done.  */
978     return 0;
979
980   groups = 0;
981   while (intdig_max > (unsigned int) *grouping)
982     {
983       ++groups;
984       intdig_max -= *grouping++;
985
986       if (*grouping == CHAR_MAX || *grouping < 0)
987         /* No more grouping should be done.  */
988         break;
989       else if (*grouping == 0)
990         {
991           /* Same grouping repeats.  */
992           groups += intdig_max / grouping[-1];
993           break;
994         }
995     }
996
997   return groups;
998 }
999
1000 /* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
1001    There is guaranteed enough space past BUFEND to extend it.
1002    Return the new end of buffer.  */
1003
1004 static char *
1005 group_number (char *buf, char *bufend, unsigned int intdig_no,
1006               const char *grouping, wchar_t thousands_sep)
1007 {
1008   unsigned int groups = __guess_grouping (intdig_no, grouping, thousands_sep);
1009   char *p;
1010
1011   if (groups == 0)
1012     return bufend;
1013
1014   /* Move the fractional part down.  */
1015   memmove (buf + intdig_no + groups, buf + intdig_no,
1016            bufend - (buf + intdig_no));
1017
1018   p = buf + intdig_no + groups - 1;
1019   do
1020     {
1021       unsigned int len = *grouping++;
1022       do
1023         *p-- = buf[--intdig_no];
1024       while (--len > 0);
1025       *p-- = thousands_sep;
1026
1027       if (*grouping == CHAR_MAX || *grouping < 0)
1028         /* No more grouping should be done.  */
1029         break;
1030       else if (*grouping == 0)
1031         /* Same grouping repeats.  */
1032         --grouping;
1033     } while (intdig_no > (unsigned int) *grouping);
1034
1035   /* Copy the remaining ungrouped digits.  */
1036   do
1037     *p-- = buf[--intdig_no];
1038   while (p > buf);
1039
1040   return bufend + groups;
1041 }