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