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