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