stdio-common/printf_fp.c

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