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