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