stdio-common/printf_fp.c

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