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