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