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