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