libquadmath/printf/printf_fp.c

   1 /* Floating point output for `printf'.
   2    Copyright (C) 1995-2003, 2006-2008, 2011 Free Software Foundation, Inc.
   3
   4    This file is part of the GNU C Library.
   5    Written by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
   6
   7    The GNU C Library is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU Lesser General Public
   9    License as published by the Free Software Foundation; either
  10    version 2.1 of the 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    Lesser General Public License for more details.
  16
  17    You should have received a copy of the GNU Lesser General Public
  18    License along with the GNU C Library; if not, write to the Free
  19    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  20    02111-1307 USA.  */
  21
  22 #include <config.h>
  23 #include <float.h>
  24 #include <math.h>
  25 #include <string.h>
  26 #include <unistd.h>
  27 #include <stdlib.h>
  28 #define NDEBUG
  29 #include <assert.h>
  30 #ifdef HAVE_ERRNO_H
  31 #include <errno.h>
  32 #endif
  33 #include <stdio.h>
  34 #include <stdarg.h>
  35 #include "quadmath-printf.h"
  36 #include "fpioconst.h"
  37
  38 #ifdef USE_I18N_NUMBER_H
  39 #include "_i18n_number.h"
  40 #else
  41 #define USE_I18N_NUMBER_H 0
  42 #endif
  43
  44 \f
  45 /* Macros for doing the actual output.  */
  46
  47 #define outchar(ch)                                                           \
  48   do                                                                          \
  49     {                                                                         \
  50       register const int outc = (ch);                                         \
  51       if (PUTC (outc, fp) == EOF)                                             \
  52         {                                                                     \
  53           if (buffer_malloced)                                                \
  54             free (wbuffer);                                                   \
  55           return -1;                                                          \
  56         }                                                                     \
  57       ++done;                                                                 \
  58     } while (0)
  59
  60 #define PRINT(ptr, wptr, len)                                                 \
  61   do                                                                          \
  62     {                                                                         \
  63       register size_t outlen = (len);                                         \
  64       if (len > 20)                                                           \
  65         {                                                                     \
  66           if (PUT (fp, wide ? (const char *) wptr : ptr, outlen) != outlen)   \
  67             {                                                                 \
  68               if (buffer_malloced)                                            \
  69                 free (wbuffer);                                               \
  70               return -1;                                                      \
  71             }                                                                 \
  72           ptr += outlen;                                                      \
  73           done += outlen;                                                     \
  74         }                                                                     \
  75       else                                                                    \
  76         {                                                                     \
  77           if (wide)                                                           \
  78             while (outlen-- > 0)                                              \
  79               outchar (*wptr++);                                              \
  80           else                                                                \
  81             while (outlen-- > 0)                                              \
  82               outchar (*ptr++);                                               \
  83         }                                                                     \
  84     } while (0)
  85
  86 #define PADN(ch, len)                                                         \
  87   do                                                                          \
  88     {                                                                         \
  89       if (PAD (fp, ch, len) != len)                                           \
  90         {                                                                     \
  91           if (buffer_malloced)                                                \
  92             free (wbuffer);                                                   \
  93           return -1;                                                          \
  94         }                                                                     \
  95       done += len;                                                            \
  96     }                                                                         \
  97   while (0)
  98
  99 \f
 100 /* We use the GNU MP library to handle large numbers.
 101
 102    An MP variable occupies a varying number of entries in its array.  We keep
 103    track of this number for efficiency reasons.  Otherwise we would always
 104    have to process the whole array.  */
 105 #define MPN_VAR(name) mp_limb_t *name; mp_size_t name##size
 106
 107 #define MPN_ASSIGN(dst,src)                                                   \
 108   memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
 109 #define MPN_GE(u,v) \
 110   (u##size > v##size || (u##size == v##size && mpn_cmp (u, v, u##size) >= 0))
 111
 112 extern mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size,
 113                                      int *expt, int *is_neg,
 114                                      __float128 value) attribute_hidden;
 115 static unsigned int guess_grouping (unsigned int intdig_max,
 116                                     const char *grouping);
 117
 118
 119 static wchar_t *group_number (wchar_t *buf, wchar_t *bufend,
 120                               unsigned int intdig_no, const char *grouping,
 121                               wchar_t thousands_sep, int ngroups);
 122
 123
 124 int
 125 __quadmath_printf_fp (struct __quadmath_printf_file *fp,
 126                       const struct printf_info *info,
 127                       const void *const *args)
 128 {
 129   /* The floating-point value to output.  */
 130   __float128 fpnum;
 131
 132   /* Locale-dependent representation of decimal point.  */
 133   const char *decimal;
 134   wchar_t decimalwc;
 135
 136   /* Locale-dependent thousands separator and grouping specification.  */
 137   const char *thousands_sep = NULL;
 138   wchar_t thousands_sepwc = L_('\0');
 139   const char *grouping;
 140
 141   /* "NaN" or "Inf" for the special cases.  */
 142   const char *special = NULL;
 143   const wchar_t *wspecial = NULL;
 144
 145   /* We need just a few limbs for the input before shifting to the right
 146      position.  */
 147   mp_limb_t fp_input[(FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB];
 148   /* We need to shift the contents of fp_input by this amount of bits.  */
 149   int to_shift = 0;
 150
 151   /* The fraction of the floting-point value in question  */
 152   MPN_VAR(frac);
 153   /* and the exponent.  */
 154   int exponent;
 155   /* Sign of the exponent.  */
 156   int expsign = 0;
 157   /* Sign of float number.  */
 158   int is_neg = 0;
 159
 160   /* Scaling factor.  */
 161   MPN_VAR(scale);
 162
 163   /* Temporary bignum value.  */
 164   MPN_VAR(tmp);
 165
 166   /* Digit which is result of last hack_digit() call.  */
 167   wchar_t digit;
 168
 169   /* The type of output format that will be used: 'e'/'E' or 'f'.  */
 170   int type;
 171
 172   /* Counter for number of written characters.  */
 173   int done = 0;
 174
 175   /* General helper (carry limb).  */
 176   mp_limb_t cy;
 177
 178   /* Nonzero if this is output on a wide character stream.  */
 179   int wide = info->wide;
 180
 181   /* Buffer in which we produce the output.  */
 182   wchar_t *wbuffer = NULL;
 183   /* Flag whether wbuffer is malloc'ed or not.  */
 184   int buffer_malloced = 0;
 185
 186   auto wchar_t hack_digit (void);
 187
 188   wchar_t hack_digit (void)
 189     {
 190       mp_limb_t hi;
 191
 192       if (expsign != 0 && type == 'f' && exponent-- > 0)
 193         hi = 0;
 194       else if (scalesize == 0)
 195         {
 196           hi = frac[fracsize - 1];
 197           frac[fracsize - 1] = mpn_mul_1 (frac, frac, fracsize - 1, 10);
 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 = scalesize;
 210               while (fracsize != 0 && frac[fracsize - 1] == 0)
 211                 --fracsize;
 212               if (fracsize == 0)
 213                 {
 214                   /* We're not prepared for an mpn variable with zero
 215                      limbs.  */
 216                   fracsize = 1;
 217                   return L_('0') + hi;
 218                 }
 219             }
 220
 221           mp_limb_t _cy = mpn_mul_1 (frac, frac, fracsize, 10);
 222           if (_cy != 0)
 223             frac[fracsize++] = _cy;
 224         }
 225
 226       return L_('0') + hi;
 227     }
 228
 229   /* Figure out the decimal point character.  */
 230 #ifdef USE_LOCALE_SUPPORT
 231   if (info->extra == 0)
 232     {
 233       decimal = nl_langinfo (DECIMAL_POINT);
 234       decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
 235     }
 236   else
 237     {
 238       decimal = nl_langinfo (MON_DECIMAL_POINT);
 239       if (*decimal == '\0')
 240         decimal = nl_langinfo (DECIMAL_POINT);
 241       decimalwc = nl_langinfo_wc (_NL_MONETARY_DECIMAL_POINT_WC);
 242       if (decimalwc == L_('\0'))
 243         decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
 244     }
 245   /* The decimal point character must not be zero.  */
 246   assert (*decimal != '\0');
 247   assert (decimalwc != L_('\0'));
 248 #else
 249   decimal = ".";
 250   decimalwc = L_('.');
 251 #endif
 252
 253 #ifdef USE_LOCALE_SUPPORT
 254   if (info->group)
 255     {
 256       if (info->extra == 0)
 257         grouping = nl_langinfo (GROUPING);
 258       else
 259         grouping = nl_langinfo (MON_GROUPING);
 260
 261       if (*grouping <= 0 || *grouping == CHAR_MAX)
 262         grouping = NULL;
 263       else
 264         {
 265           /* Figure out the thousands separator character.  */
 266           if (wide)
 267             {
 268               if (info->extra == 0)
 269                 thousands_sepwc = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC);
 270               else
 271                 thousands_sepwc = nl_langinfo_wc (_NL_MONETARY_THOUSANDS_SEP_WC);
 272             }
 273           else
 274             {
 275               if (info->extra == 0)
 276                 thousands_sep = nl_langinfo (THOUSANDS_SEP);
 277               else
 278                 thousands_sep = nl_langinfo (MON_THOUSANDS_SEP);
 279             }
 280
 281           if ((wide && thousands_sepwc == L_('\0'))
 282               || (! wide && *thousands_sep == '\0'))
 283             grouping = NULL;
 284           else if (thousands_sepwc == L_('\0'))
 285             /* If we are printing multibyte characters and there is a
 286                multibyte representation for the thousands separator,
 287                we must ensure the wide character thousands separator
 288                is available, even if it is fake.  */
 289             thousands_sepwc = (wchar_t) 0xfffffffe;
 290         }
 291     }
 292   else
 293 #endif
 294     grouping = NULL;
 295
 296   /* Fetch the argument value.  */
 297     {
 298       fpnum = **(const __float128 **) args[0];
 299
 300       /* Check for special values: not a number or infinity.  */
 301       if (isnanq (fpnum))
 302         {
 303           ieee854_float128 u = { .value = fpnum };
 304           is_neg = u.ieee.negative != 0;
 305           if (isupper (info->spec))
 306             {
 307               special = "NAN";
 308               wspecial = L_("NAN");
 309             }
 310             else
 311               {
 312                 special = "nan";
 313                 wspecial = L_("nan");
 314               }
 315         }
 316       else if (isinfq (fpnum))
 317         {
 318           is_neg = fpnum < 0;
 319           if (isupper (info->spec))
 320             {
 321               special = "INF";
 322               wspecial = L_("INF");
 323             }
 324           else
 325             {
 326               special = "inf";
 327               wspecial = L_("inf");
 328             }
 329         }
 330       else
 331         {
 332           fracsize = mpn_extract_flt128 (fp_input,
 333                                          (sizeof (fp_input) /
 334                                           sizeof (fp_input[0])),
 335                                          &exponent, &is_neg, fpnum);
 336           to_shift = 1 + fracsize * BITS_PER_MP_LIMB - FLT128_MANT_DIG;
 337         }
 338     }
 339
 340   if (special)
 341     {
 342       int width = info->width;
 343
 344       if (is_neg || info->showsign || info->space)
 345         --width;
 346       width -= 3;
 347
 348       if (!info->left && width > 0)
 349         PADN (' ', width);
 350
 351       if (is_neg)
 352         outchar ('-');
 353       else if (info->showsign)
 354         outchar ('+');
 355       else if (info->space)
 356         outchar (' ');
 357
 358       PRINT (special, wspecial, 3);
 359
 360       if (info->left && width > 0)
 361         PADN (' ', width);
 362
 363       return done;
 364     }
 365
 366
 367   /* We need three multiprecision variables.  Now that we have the exponent
 368      of the number we can allocate the needed memory.  It would be more
 369      efficient to use variables of the fixed maximum size but because this
 370      would be really big it could lead to memory problems.  */
 371   {
 372     mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1)
 373                              / BITS_PER_MP_LIMB
 374                              + (FLT128_MANT_DIG / BITS_PER_MP_LIMB > 2 ? 8 : 4))
 375                             * sizeof (mp_limb_t);
 376     frac = (mp_limb_t *) alloca (bignum_size);
 377     tmp = (mp_limb_t *) alloca (bignum_size);
 378     scale = (mp_limb_t *) alloca (bignum_size);
 379   }
 380
 381   /* We now have to distinguish between numbers with positive and negative
 382      exponents because the method used for the one is not applicable/efficient
 383      for the other.  */
 384   scalesize = 0;
 385   if (exponent > 2)
 386     {
 387       /* |FP| >= 8.0.  */
 388       int scaleexpo = 0;
 389       int explog = FLT128_MAX_10_EXP_LOG;
 390       int exp10 = 0;
 391       const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
 392       int cnt_h, cnt_l, i;
 393
 394       if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0)
 395         {
 396           MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
 397                          fp_input, fracsize);
 398           fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
 399         }
 400       else
 401         {
 402           cy = mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
 403                            fp_input, fracsize,
 404                            (exponent + to_shift) % BITS_PER_MP_LIMB);
 405           fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
 406           if (cy)
 407             frac[fracsize++] = cy;
 408         }
 409       MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB);
 410
 411       assert (powers > &_fpioconst_pow10[0]);
 412       do
 413         {
 414           --powers;
 415
 416           /* The number of the product of two binary numbers with n and m
 417              bits respectively has m+n or m+n-1 bits.   */
 418           if (exponent >= scaleexpo + powers->p_expo - 1)
 419             {
 420               if (scalesize == 0)
 421                 {
 422                   if (FLT128_MANT_DIG > _FPIO_CONST_OFFSET * BITS_PER_MP_LIMB)
 423                     {
 424 #define _FPIO_CONST_SHIFT \
 425   (((FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) \
 426    - _FPIO_CONST_OFFSET)
 427                       /* 64bit const offset is not enough for
 428                          IEEE quad long double.  */
 429                       tmpsize = powers->arraysize + _FPIO_CONST_SHIFT;
 430                       memcpy (tmp + _FPIO_CONST_SHIFT,
 431                               &__tens[powers->arrayoff],
 432                               tmpsize * sizeof (mp_limb_t));
 433                       MPN_ZERO (tmp, _FPIO_CONST_SHIFT);
 434                       /* Adjust exponent, as scaleexpo will be this much
 435                          bigger too.  */
 436                       exponent += _FPIO_CONST_SHIFT * BITS_PER_MP_LIMB;
 437                     }
 438                   else
 439                     {
 440                       tmpsize = powers->arraysize;
 441                       memcpy (tmp, &__tens[powers->arrayoff],
 442                               tmpsize * sizeof (mp_limb_t));
 443                     }
 444                 }
 445               else
 446                 {
 447                   cy = mpn_mul (tmp, scale, scalesize,
 448                                 &__tens[powers->arrayoff
 449                                         + _FPIO_CONST_OFFSET],
 450                                 powers->arraysize - _FPIO_CONST_OFFSET);
 451                   tmpsize = scalesize + powers->arraysize - _FPIO_CONST_OFFSET;
 452                   if (cy == 0)
 453                     --tmpsize;
 454                 }
 455
 456               if (MPN_GE (frac, tmp))
 457                 {
 458                   int cnt;
 459                   MPN_ASSIGN (scale, tmp);
 460                   count_leading_zeros (cnt, scale[scalesize - 1]);
 461                   scaleexpo = (scalesize - 2) * BITS_PER_MP_LIMB - cnt - 1;
 462                   exp10 |= 1 << explog;
 463                 }
 464             }
 465           --explog;
 466         }
 467       while (powers > &_fpioconst_pow10[0]);
 468       exponent = exp10;
 469
 470       /* Optimize number representations.  We want to represent the numbers
 471          with the lowest number of bytes possible without losing any
 472          bytes. Also the highest bit in the scaling factor has to be set
 473          (this is a requirement of the MPN division routines).  */
 474       if (scalesize > 0)
 475         {
 476           /* Determine minimum number of zero bits at the end of
 477              both numbers.  */
 478           for (i = 0; scale[i] == 0 && frac[i] == 0; i++)
 479             ;
 480
 481           /* Determine number of bits the scaling factor is misplaced.  */
 482           count_leading_zeros (cnt_h, scale[scalesize - 1]);
 483
 484           if (cnt_h == 0)
 485             {
 486               /* The highest bit of the scaling factor is already set.  So
 487                  we only have to remove the trailing empty limbs.  */
 488               if (i > 0)
 489                 {
 490                   MPN_COPY_INCR (scale, scale + i, scalesize - i);
 491                   scalesize -= i;
 492                   MPN_COPY_INCR (frac, frac + i, fracsize - i);
 493                   fracsize -= i;
 494                 }
 495             }
 496           else
 497             {
 498               if (scale[i] != 0)
 499                 {
 500                   count_trailing_zeros (cnt_l, scale[i]);
 501                   if (frac[i] != 0)
 502                     {
 503                       int cnt_l2;
 504                       count_trailing_zeros (cnt_l2, frac[i]);
 505                       if (cnt_l2 < cnt_l)
 506                         cnt_l = cnt_l2;
 507                     }
 508                 }
 509               else
 510                 count_trailing_zeros (cnt_l, frac[i]);
 511
 512               /* Now shift the numbers to their optimal position.  */
 513               if (i == 0 && BITS_PER_MP_LIMB - cnt_h > cnt_l)
 514                 {
 515                   /* We cannot save any memory.  So just roll both numbers
 516                      so that the scaling factor has its highest bit set.  */
 517
 518                   (void) mpn_lshift (scale, scale, scalesize, cnt_h);
 519                   cy = mpn_lshift (frac, frac, fracsize, cnt_h);
 520                   if (cy != 0)
 521                     frac[fracsize++] = cy;
 522                 }
 523               else if (BITS_PER_MP_LIMB - cnt_h <= cnt_l)
 524                 {
 525                   /* We can save memory by removing the trailing zero limbs
 526                      and by packing the non-zero limbs which gain another
 527                      free one. */
 528
 529                   (void) mpn_rshift (scale, scale + i, scalesize - i,
 530                                      BITS_PER_MP_LIMB - cnt_h);
 531                   scalesize -= i + 1;
 532                   (void) mpn_rshift (frac, frac + i, fracsize - i,
 533                                      BITS_PER_MP_LIMB - cnt_h);
 534                   fracsize -= frac[fracsize - i - 1] == 0 ? i + 1 : i;
 535                 }
 536               else
 537                 {
 538                   /* We can only save the memory of the limbs which are zero.
 539                      The non-zero parts occupy the same number of limbs.  */
 540
 541                   (void) mpn_rshift (scale, scale + (i - 1),
 542                                      scalesize - (i - 1),
 543                                      BITS_PER_MP_LIMB - cnt_h);
 544                   scalesize -= i;
 545                   (void) mpn_rshift (frac, frac + (i - 1),
 546                                      fracsize - (i - 1),
 547                                      BITS_PER_MP_LIMB - cnt_h);
 548                   fracsize -= frac[fracsize - (i - 1) - 1] == 0 ? i : i - 1;
 549                 }
 550             }
 551         }
 552     }
 553   else if (exponent < 0)
 554     {
 555       /* |FP| < 1.0.  */
 556       int exp10 = 0;
 557       int explog = FLT128_MAX_10_EXP_LOG;
 558       const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
 559
 560       /* Now shift the input value to its right place.  */
 561       cy = mpn_lshift (frac, fp_input, fracsize, to_shift);
 562       frac[fracsize++] = cy;
 563       assert (cy == 1 || (frac[fracsize - 2] == 0 && frac[0] == 0));
 564
 565       expsign = 1;
 566       exponent = -exponent;
 567
 568       assert (powers != &_fpioconst_pow10[0]);
 569       do
 570         {
 571           --powers;
 572
 573           if (exponent >= powers->m_expo)
 574             {
 575               int i, incr, cnt_h, cnt_l;
 576               mp_limb_t topval[2];
 577
 578               /* The mpn_mul function expects the first argument to be
 579                  bigger than the second.  */
 580               if (fracsize < powers->arraysize - _FPIO_CONST_OFFSET)
 581                 cy = mpn_mul (tmp, &__tens[powers->arrayoff
 582                                            + _FPIO_CONST_OFFSET],
 583                               powers->arraysize - _FPIO_CONST_OFFSET,
 584                               frac, fracsize);
 585               else
 586                 cy = mpn_mul (tmp, frac, fracsize,
 587                               &__tens[powers->arrayoff + _FPIO_CONST_OFFSET],
 588                               powers->arraysize - _FPIO_CONST_OFFSET);
 589               tmpsize = fracsize + powers->arraysize - _FPIO_CONST_OFFSET;
 590               if (cy == 0)
 591                 --tmpsize;
 592
 593               count_leading_zeros (cnt_h, tmp[tmpsize - 1]);
 594               incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB
 595                      + BITS_PER_MP_LIMB - 1 - cnt_h;
 596
 597               assert (incr <= powers->p_expo);
 598
 599               /* If we increased the exponent by exactly 3 we have to test
 600                  for overflow.  This is done by comparing with 10 shifted
 601                  to the right position.  */
 602               if (incr == exponent + 3)
 603                 {
 604                   if (cnt_h <= BITS_PER_MP_LIMB - 4)
 605                     {
 606                       topval[0] = 0;
 607                       topval[1]
 608                         = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4 - cnt_h);
 609                     }
 610                   else
 611                     {
 612                       topval[0] = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4);
 613                       topval[1] = 0;
 614                       (void) mpn_lshift (topval, topval, 2,
 615                                          BITS_PER_MP_LIMB - cnt_h);
 616                     }
 617                 }
 618
 619               /* We have to be careful when multiplying the last factor.
 620                  If the result is greater than 1.0 be have to test it
 621                  against 10.0.  If it is greater or equal to 10.0 the
 622                  multiplication was not valid.  This is because we cannot
 623                  determine the number of bits in the result in advance.  */
 624               if (incr < exponent + 3
 625                   || (incr == exponent + 3 &&
 626                       (tmp[tmpsize - 1] < topval[1]
 627                        || (tmp[tmpsize - 1] == topval[1]
 628                            && tmp[tmpsize - 2] < topval[0]))))
 629                 {
 630                   /* The factor is right.  Adapt binary and decimal
 631                      exponents.  */
 632                   exponent -= incr;
 633                   exp10 |= 1 << explog;
 634
 635                   /* If this factor yields a number greater or equal to
 636                      1.0, we must not shift the non-fractional digits down. */
 637                   if (exponent < 0)
 638                     cnt_h += -exponent;
 639
 640                   /* Now we optimize the number representation.  */
 641                   for (i = 0; tmp[i] == 0; ++i);
 642                   if (cnt_h == BITS_PER_MP_LIMB - 1)
 643                     {
 644                       MPN_COPY (frac, tmp + i, tmpsize - i);
 645                       fracsize = tmpsize - i;
 646                     }
 647                   else
 648                     {
 649                       count_trailing_zeros (cnt_l, tmp[i]);
 650
 651                       /* Now shift the numbers to their optimal position.  */
 652                       if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l)
 653                         {
 654                           /* We cannot save any memory.  Just roll the
 655                              number so that the leading digit is in a
 656                              separate limb.  */
 657
 658                           cy = mpn_lshift (frac, tmp, tmpsize, cnt_h + 1);
 659                           fracsize = tmpsize + 1;
 660                           frac[fracsize - 1] = cy;
 661                         }
 662                       else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l)
 663                         {
 664                           (void) mpn_rshift (frac, tmp + i, tmpsize - i,
 665                                              BITS_PER_MP_LIMB - 1 - cnt_h);
 666                           fracsize = tmpsize - i;
 667                         }
 668                       else
 669                         {
 670                           /* We can only save the memory of the limbs which
 671                              are zero.  The non-zero parts occupy the same
 672                              number of limbs.  */
 673
 674                           (void) mpn_rshift (frac, tmp + (i - 1),
 675                                              tmpsize - (i - 1),
 676                                              BITS_PER_MP_LIMB - 1 - cnt_h);
 677                           fracsize = tmpsize - (i - 1);
 678                         }
 679                     }
 680                 }
 681             }
 682           --explog;
 683         }
 684       while (powers != &_fpioconst_pow10[1] && exponent > 0);
 685       /* All factors but 10^-1 are tested now.  */
 686       if (exponent > 0)
 687         {
 688           int cnt_l;
 689
 690           cy = mpn_mul_1 (tmp, frac, fracsize, 10);
 691           tmpsize = fracsize;
 692           assert (cy == 0 || tmp[tmpsize - 1] < 20);
 693
 694           count_trailing_zeros (cnt_l, tmp[0]);
 695           if (cnt_l < MIN (4, exponent))
 696             {
 697               cy = mpn_lshift (frac, tmp, tmpsize,
 698                                BITS_PER_MP_LIMB - MIN (4, exponent));
 699               if (cy != 0)
 700                 frac[tmpsize++] = cy;
 701             }
 702           else
 703             (void) mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent));
 704           fracsize = tmpsize;
 705           exp10 |= 1;
 706           assert (frac[fracsize - 1] < 10);
 707         }
 708       exponent = exp10;
 709     }
 710   else
 711     {
 712       /* This is a special case.  We don't need a factor because the
 713          numbers are in the range of 1.0 <= |fp| < 8.0.  We simply
 714          shift it to the right place and divide it by 1.0 to get the
 715          leading digit.  (Of course this division is not really made.)  */
 716       assert (0 <= exponent && exponent < 3 &&
 717               exponent + to_shift < BITS_PER_MP_LIMB);
 718
 719       /* Now shift the input value to its right place.  */
 720       cy = mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift));
 721       frac[fracsize++] = cy;
 722       exponent = 0;
 723     }
 724
 725   {
 726     int width = info->width;
 727     wchar_t *wstartp, *wcp;
 728     size_t chars_needed;
 729     int expscale;
 730     int intdig_max, intdig_no = 0;
 731     int fracdig_min;
 732     int fracdig_max;
 733     int dig_max;
 734     int significant;
 735     int ngroups = 0;
 736     char spec = tolower (info->spec);
 737
 738     if (spec == 'e')
 739       {
 740         type = info->spec;
 741         intdig_max = 1;
 742         fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
 743         chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4;
 744         /*             d   .     ddd         e   +-  ddd  */
 745         dig_max = __INT_MAX__;          /* Unlimited.  */
 746         significant = 1;                /* Does not matter here.  */
 747       }
 748     else if (spec == 'f')
 749       {
 750         type = 'f';
 751         fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
 752         dig_max = __INT_MAX__;          /* Unlimited.  */
 753         significant = 1;                /* Does not matter here.  */
 754         if (expsign == 0)
 755           {
 756             intdig_max = exponent + 1;
 757             /* This can be really big!  */  /* XXX Maybe malloc if too big? */
 758             chars_needed = (size_t) exponent + 1 + 1 + (size_t) fracdig_max;
 759           }
 760         else
 761           {
 762             intdig_max = 1;
 763             chars_needed = 1 + 1 + (size_t) fracdig_max;
 764           }
 765       }
 766     else
 767       {
 768         dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec);
 769         if ((expsign == 0 && exponent >= dig_max)
 770             || (expsign != 0 && exponent > 4))
 771           {
 772             if ('g' - 'G' == 'e' - 'E')
 773               type = 'E' + (info->spec - 'G');
 774             else
 775               type = isupper (info->spec) ? 'E' : 'e';
 776             fracdig_max = dig_max - 1;
 777             intdig_max = 1;
 778             chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4;
 779           }
 780         else
 781           {
 782             type = 'f';
 783             intdig_max = expsign == 0 ? exponent + 1 : 0;
 784             fracdig_max = dig_max - intdig_max;
 785             /* We need space for the significant digits and perhaps
 786                for leading zeros when < 1.0.  The number of leading
 787                zeros can be as many as would be required for
 788                exponential notation with a negative two-digit
 789                exponent, which is 4.  */
 790             chars_needed = (size_t) dig_max + 1 + 4;
 791           }
 792         fracdig_min = info->alt ? fracdig_max : 0;
 793         significant = 0;                /* We count significant digits.  */
 794       }
 795
 796     if (grouping)
 797       {
 798         /* Guess the number of groups we will make, and thus how
 799            many spaces we need for separator characters.  */
 800         ngroups = guess_grouping (intdig_max, grouping);
 801         /* Allocate one more character in case rounding increases the
 802            number of groups.  */
 803         chars_needed += ngroups + 1;
 804       }
 805
 806     /* Allocate buffer for output.  We need two more because while rounding
 807        it is possible that we need two more characters in front of all the
 808        other output.  If the amount of memory we have to allocate is too
 809        large use `malloc' instead of `alloca'.  */
 810     if (__builtin_expect (chars_needed >= (size_t) -1 / sizeof (wchar_t) - 2
 811                           || chars_needed < fracdig_max, 0))
 812       {
 813         /* Some overflow occurred.  */
 814 #if defined HAVE_ERRNO_H && defined ERANGE
 815         errno = ERANGE;
 816 #endif
 817         return -1;
 818       }
 819     size_t wbuffer_to_alloc = (2 + chars_needed) * sizeof (wchar_t);
 820     buffer_malloced = wbuffer_to_alloc >= 4096;
 821     if (__builtin_expect (buffer_malloced, 0))
 822       {
 823         wbuffer = (wchar_t *) malloc (wbuffer_to_alloc);
 824         if (wbuffer == NULL)
 825           /* Signal an error to the caller.  */
 826           return -1;
 827       }
 828     else
 829       wbuffer = (wchar_t *) alloca (wbuffer_to_alloc);
 830     wcp = wstartp = wbuffer + 2;        /* Let room for rounding.  */
 831
 832     /* Do the real work: put digits in allocated buffer.  */
 833     if (expsign == 0 || type != 'f')
 834       {
 835         assert (expsign == 0 || intdig_max == 1);
 836         while (intdig_no < intdig_max)
 837           {
 838             ++intdig_no;
 839             *wcp++ = hack_digit ();
 840           }
 841         significant = 1;
 842         if (info->alt
 843             || fracdig_min > 0
 844             || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0)))
 845           *wcp++ = decimalwc;
 846       }
 847     else
 848       {
 849         /* |fp| < 1.0 and the selected type is 'f', so put "0."
 850            in the buffer.  */
 851         *wcp++ = L_('0');
 852         --exponent;
 853         *wcp++ = decimalwc;
 854       }
 855
 856     /* Generate the needed number of fractional digits.  */
 857     int fracdig_no = 0;
 858     int added_zeros = 0;
 859     while (fracdig_no < fracdig_min + added_zeros
 860            || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
 861       {
 862         ++fracdig_no;
 863         *wcp = hack_digit ();
 864         if (*wcp++ != L_('0'))
 865           significant = 1;
 866         else if (significant == 0)
 867           {
 868             ++fracdig_max;
 869             if (fracdig_min > 0)
 870               ++added_zeros;
 871           }
 872       }
 873
 874     /* Do rounding.  */
 875     digit = hack_digit ();
 876     if (digit > L_('4'))
 877       {
 878         wchar_t *wtp = wcp;
 879
 880         if (digit == L_('5')
 881             && ((*(wcp - 1) != decimalwc && (*(wcp - 1) & 1) == 0)
 882                 || ((*(wcp - 1) == decimalwc && (*(wcp - 2) & 1) == 0))))
 883           {
 884             /* This is the critical case.        */
 885             if (fracsize == 1 && frac[0] == 0)
 886               /* Rest of the number is zero -> round to even.
 887                  (IEEE 754-1985 4.1 says this is the default rounding.)  */
 888               goto do_expo;
 889             else if (scalesize == 0)
 890               {
 891                 /* Here we have to see whether all limbs are zero since no
 892                    normalization happened.  */
 893                 size_t lcnt = fracsize;
 894                 while (lcnt >= 1 && frac[lcnt - 1] == 0)
 895                   --lcnt;
 896                 if (lcnt == 0)
 897                   /* Rest of the number is zero -> round to even.
 898                      (IEEE 754-1985 4.1 says this is the default rounding.)  */
 899                   goto do_expo;
 900               }
 901           }
 902
 903         if (fracdig_no > 0)
 904           {
 905             /* Process fractional digits.  Terminate if not rounded or
 906                radix character is reached.  */
 907             int removed = 0;
 908             while (*--wtp != decimalwc && *wtp == L_('9'))
 909               {
 910                 *wtp = L_('0');
 911                 ++removed;
 912               }
 913             if (removed == fracdig_min && added_zeros > 0)
 914               --added_zeros;
 915             if (*wtp != decimalwc)
 916               /* Round up.  */
 917               (*wtp)++;
 918             else if (__builtin_expect (spec == 'g' && type == 'f' && info->alt
 919                                        && wtp == wstartp + 1
 920                                        && wstartp[0] == L_('0'),
 921                                        0))
 922               /* This is a special case: the rounded number is 1.0,
 923                  the format is 'g' or 'G', and the alternative format
 924                  is selected.  This means the result must be "1.".  */
 925               --added_zeros;
 926           }
 927
 928         if (fracdig_no == 0 || *wtp == decimalwc)
 929           {
 930             /* Round the integer digits.  */
 931             if (*(wtp - 1) == decimalwc)
 932               --wtp;
 933
 934             while (--wtp >= wstartp && *wtp == L_('9'))
 935               *wtp = L_('0');
 936
 937             if (wtp >= wstartp)
 938               /* Round up.  */
 939               (*wtp)++;
 940             else
 941               /* It is more critical.  All digits were 9's.  */
 942               {
 943                 if (type != 'f')
 944                   {
 945                     *wstartp = '1';
 946                     exponent += expsign == 0 ? 1 : -1;
 947
 948                     /* The above exponent adjustment could lead to 1.0e-00,
 949                        e.g. for 0.999999999.  Make sure exponent 0 always
 950                        uses + sign.  */
 951                     if (exponent == 0)
 952                       expsign = 0;
 953                   }
 954                 else if (intdig_no == dig_max)
 955                   {
 956                     /* This is the case where for type %g the number fits
 957                        really in the range for %f output but after rounding
 958                        the number of digits is too big.  */
 959                     *--wstartp = decimalwc;
 960                     *--wstartp = L_('1');
 961
 962                     if (info->alt || fracdig_no > 0)
 963                       {
 964                         /* Overwrite the old radix character.  */
 965                         wstartp[intdig_no + 2] = L_('0');
 966                         ++fracdig_no;
 967                       }
 968
 969                     fracdig_no += intdig_no;
 970                     intdig_no = 1;
 971                     fracdig_max = intdig_max - intdig_no;
 972                     ++exponent;
 973                     /* Now we must print the exponent.  */
 974                     type = isupper (info->spec) ? 'E' : 'e';
 975                   }
 976                 else
 977                   {
 978                     /* We can simply add another another digit before the
 979                        radix.  */
 980                     *--wstartp = L_('1');
 981                     ++intdig_no;
 982                   }
 983
 984                 /* While rounding the number of digits can change.
 985                    If the number now exceeds the limits remove some
 986                    fractional digits.  */
 987                 if (intdig_no + fracdig_no > dig_max)
 988                   {
 989                     wcp -= intdig_no + fracdig_no - dig_max;
 990                     fracdig_no -= intdig_no + fracdig_no - dig_max;
 991                   }
 992               }
 993           }
 994       }
 995
 996   do_expo:
 997     /* Now remove unnecessary '0' at the end of the string.  */
 998     while (fracdig_no > fracdig_min + added_zeros && *(wcp - 1) == L_('0'))
 999       {
1000         --wcp;
1001         --fracdig_no;
1002       }
1003     /* If we eliminate all fractional digits we perhaps also can remove
1004        the radix character.  */
1005     if (fracdig_no == 0 && !info->alt && *(wcp - 1) == decimalwc)
1006       --wcp;
1007
1008     if (grouping)
1009       {
1010         /* Rounding might have changed the number of groups.  We allocated
1011            enough memory but we need here the correct number of groups.  */
1012         if (intdig_no != intdig_max)
1013           ngroups = guess_grouping (intdig_no, grouping);
1014
1015         /* Add in separator characters, overwriting the same buffer.  */
1016         wcp = group_number (wstartp, wcp, intdig_no, grouping, thousands_sepwc,
1017                             ngroups);
1018       }
1019
1020     /* Write the exponent if it is needed.  */
1021     if (type != 'f')
1022       {
1023         if (__builtin_expect (expsign != 0 && exponent == 4 && spec == 'g', 0))
1024           {
1025             /* This is another special case.  The exponent of the number is
1026                really smaller than -4, which requires the 'e'/'E' format.
1027                But after rounding the number has an exponent of -4.  */
1028             assert (wcp >= wstartp + 1);
1029             assert (wstartp[0] == L_('1'));
1030             memcpy (wstartp, L_("0.0001"), 6 * sizeof (wchar_t));
1031             wstartp[1] = decimalwc;
1032             if (wcp >= wstartp + 2)
1033               {
1034                 size_t cnt;
1035                 for (cnt = 0; cnt < wcp - (wstartp + 2); cnt++)
1036                   wstartp[6 + cnt] = L_('0');
1037                 wcp += 4;
1038               }
1039             else
1040               wcp += 5;
1041           }
1042         else
1043           {
1044             *wcp++ = (wchar_t) type;
1045             *wcp++ = expsign ? L_('-') : L_('+');
1046
1047             /* Find the magnitude of the exponent.      */
1048             expscale = 10;
1049             while (expscale <= exponent)
1050               expscale *= 10;
1051
1052             if (exponent < 10)
1053               /* Exponent always has at least two digits.  */
1054               *wcp++ = L_('0');
1055             else
1056               do
1057                 {
1058                   expscale /= 10;
1059                   *wcp++ = L_('0') + (exponent / expscale);
1060                   exponent %= expscale;
1061                 }
1062               while (expscale > 10);
1063             *wcp++ = L_('0') + exponent;
1064           }
1065       }
1066
1067     /* Compute number of characters which must be filled with the padding
1068        character.  */
1069     if (is_neg || info->showsign || info->space)
1070       --width;
1071     width -= wcp - wstartp;
1072
1073     if (!info->left && info->pad != '0' && width > 0)
1074       PADN (info->pad, width);
1075
1076     if (is_neg)
1077       outchar ('-');
1078     else if (info->showsign)
1079       outchar ('+');
1080     else if (info->space)
1081       outchar (' ');
1082
1083     if (!info->left && info->pad == '0' && width > 0)
1084       PADN ('0', width);
1085
1086     {
1087       char *buffer = NULL;
1088       char *buffer_end __attribute__((__unused__)) = NULL;
1089       char *cp = NULL;
1090       char *tmpptr;
1091
1092       if (! wide)
1093         {
1094           /* Create the single byte string.  */
1095           size_t decimal_len;
1096           size_t thousands_sep_len;
1097           wchar_t *copywc;
1098 #ifdef USE_LOCALE_SUPPORT
1099           size_t factor = ((info->i18n && USE_I18N_NUMBER_H)
1100                            ? nl_langinfo_wc (_NL_CTYPE_MB_CUR_MAX)
1101                            : 1);
1102 #else
1103           size_t factor = 1;
1104 #endif
1105
1106           decimal_len = strlen (decimal);
1107
1108           if (thousands_sep == NULL)
1109             thousands_sep_len = 0;
1110           else
1111             thousands_sep_len = strlen (thousands_sep);
1112
1113           size_t nbuffer = (2 + chars_needed * factor + decimal_len
1114                             + ngroups * thousands_sep_len);
1115           if (__builtin_expect (buffer_malloced, 0))
1116             {
1117               buffer = (char *) malloc (nbuffer);
1118               if (buffer == NULL)
1119                 {
1120                   /* Signal an error to the caller.  */
1121                   free (wbuffer);
1122                   return -1;
1123                 }
1124             }
1125           else
1126             buffer = (char *) alloca (nbuffer);
1127           buffer_end = buffer + nbuffer;
1128
1129           /* Now copy the wide character string.  Since the character
1130              (except for the decimal point and thousands separator) must
1131              be coming from the ASCII range we can esily convert the
1132              string without mapping tables.  */
1133           for (cp = buffer, copywc = wstartp; copywc < wcp; ++copywc)
1134             if (*copywc == decimalwc)
1135               memcpy (cp, decimal, decimal_len), cp += decimal_len;
1136             else if (*copywc == thousands_sepwc)
1137               mempcpy (cp, thousands_sep, thousands_sep_len), cp += thousands_sep_len;
1138             else
1139               *cp++ = (char) *copywc;
1140         }
1141
1142       tmpptr = buffer;
1143 #if USE_I18N_NUMBER_H
1144       if (__builtin_expect (info->i18n, 0))
1145         {
1146           tmpptr = _i18n_number_rewrite (tmpptr, cp, buffer_end);
1147           cp = buffer_end;
1148           assert ((uintptr_t) buffer <= (uintptr_t) tmpptr);
1149           assert ((uintptr_t) tmpptr < (uintptr_t) buffer_end);
1150         }
1151 #endif
1152
1153       PRINT (tmpptr, wstartp, wide ? wcp - wstartp : cp - tmpptr);
1154
1155       /* Free the memory if necessary.  */
1156       if (__builtin_expect (buffer_malloced, 0))
1157         {
1158           free (buffer);
1159           free (wbuffer);
1160         }
1161     }
1162
1163     if (info->left && width > 0)
1164       PADN (info->pad, width);
1165   }
1166   return done;
1167 }
1168 \f
1169 /* Return the number of extra grouping characters that will be inserted
1170    into a number with INTDIG_MAX integer digits.  */
1171
1172 static unsigned int
1173 guess_grouping (unsigned int intdig_max, const char *grouping)
1174 {
1175   unsigned int groups;
1176
1177   /* We treat all negative values like CHAR_MAX.  */
1178
1179   if (*grouping == CHAR_MAX || *grouping <= 0)
1180     /* No grouping should be done.  */
1181     return 0;
1182
1183   groups = 0;
1184   while (intdig_max > (unsigned int) *grouping)
1185     {
1186       ++groups;
1187       intdig_max -= *grouping++;
1188
1189       if (*grouping == 0)
1190         {
1191           /* Same grouping repeats.  */
1192           groups += (intdig_max - 1) / grouping[-1];
1193           break;
1194         }
1195       else if (*grouping == CHAR_MAX || *grouping <= 0)
1196         /* No more grouping should be done.  */
1197         break;
1198     }
1199
1200   return groups;
1201 }
1202
1203 /* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
1204    There is guaranteed enough space past BUFEND to extend it.
1205    Return the new end of buffer.  */
1206
1207 static wchar_t *
1208 group_number (wchar_t *buf, wchar_t *bufend, unsigned int intdig_no,
1209               const char *grouping, wchar_t thousands_sep, int ngroups)
1210 {
1211   wchar_t *p;
1212
1213   if (ngroups == 0)
1214     return bufend;
1215
1216   /* Move the fractional part down.  */
1217   memmove (buf + intdig_no + ngroups, buf + intdig_no,
1218            (bufend - (buf + intdig_no)) * sizeof (wchar_t));
1219
1220   p = buf + intdig_no + ngroups - 1;
1221   do
1222     {
1223       unsigned int len = *grouping++;
1224       do
1225         *p-- = buf[--intdig_no];
1226       while (--len > 0);
1227       *p-- = thousands_sep;
1228
1229       if (*grouping == 0)
1230         /* Same grouping repeats.  */
1231         --grouping;
1232       else if (*grouping == CHAR_MAX || *grouping <= 0)
1233         /* No more grouping should be done.  */
1234         break;
1235     } while (intdig_no > (unsigned int) *grouping);
1236
1237   /* Copy the remaining ungrouped digits.  */
1238   do
1239     *p-- = buf[--intdig_no];
1240   while (p > buf);
1241
1242   return bufend + ngroups;
1243 }