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