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