1/* Copyright (C) 2007-2021 Free Software Foundation, Inc. 2 Contributed by Andy Vaught 3 Write float code factoring to this file by Jerry DeLisle 4 F2003 I/O support contributed by Jerry DeLisle 5 6This file is part of the GNU Fortran runtime library (libgfortran). 7 8Libgfortran is free software; you can redistribute it and/or modify 9it under the terms of the GNU General Public License as published by 10the Free Software Foundation; either version 3, or (at your option) 11any later version. 12 13Libgfortran is distributed in the hope that it will be useful, 14but WITHOUT ANY WARRANTY; without even the implied warranty of 15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16GNU General Public License for more details. 17 18Under Section 7 of GPL version 3, you are granted additional 19permissions described in the GCC Runtime Library Exception, version 203.1, as published by the Free Software Foundation. 21 22You should have received a copy of the GNU General Public License and 23a copy of the GCC Runtime Library Exception along with this program; 24see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 25<http://www.gnu.org/licenses/>. */ 26 27#include "config.h" 28 29typedef enum 30{ S_NONE, S_MINUS, S_PLUS } 31sign_t; 32 33/* Given a flag that indicates if a value is negative or not, return a 34 sign_t that gives the sign that we need to produce. */ 35 36static sign_t 37calculate_sign (st_parameter_dt *dtp, int negative_flag) 38{ 39 sign_t s = S_NONE; 40 41 if (negative_flag) 42 s = S_MINUS; 43 else 44 switch (dtp->u.p.sign_status) 45 { 46 case SIGN_SP: /* Show sign. */ 47 s = S_PLUS; 48 break; 49 case SIGN_SS: /* Suppress sign. */ 50 s = S_NONE; 51 break; 52 case SIGN_S: /* Processor defined. */ 53 case SIGN_UNSPECIFIED: 54 s = options.optional_plus ? S_PLUS : S_NONE; 55 break; 56 } 57 58 return s; 59} 60 61 62/* Determine the precision except for EN format. For G format, 63 determines an upper bound to be used for sizing the buffer. */ 64 65static int 66determine_precision (st_parameter_dt * dtp, const fnode * f, int len) 67{ 68 int precision = f->u.real.d; 69 70 switch (f->format) 71 { 72 case FMT_F: 73 case FMT_G: 74 precision += dtp->u.p.scale_factor; 75 break; 76 case FMT_ES: 77 /* Scale factor has no effect on output. */ 78 break; 79 case FMT_E: 80 case FMT_D: 81 /* See F2008 10.7.2.3.3.6 */ 82 if (dtp->u.p.scale_factor <= 0) 83 precision += dtp->u.p.scale_factor - 1; 84 break; 85 default: 86 return -1; 87 } 88 89 /* If the scale factor has a large negative value, we must do our 90 own rounding? Use ROUND='NEAREST', which should be what snprintf 91 is using as well. */ 92 if (precision < 0 && 93 (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED 94 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED)) 95 dtp->u.p.current_unit->round_status = ROUND_NEAREST; 96 97 /* Add extra guard digits up to at least full precision when we do 98 our own rounding. */ 99 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED 100 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED) 101 { 102 precision += 2 * len + 4; 103 if (precision < 0) 104 precision = 0; 105 } 106 107 return precision; 108} 109 110 111/* Build a real number according to its format which is FMT_G free. */ 112 113static void 114build_float_string (st_parameter_dt *dtp, const fnode *f, char *buffer, 115 size_t size, int nprinted, int precision, int sign_bit, 116 bool zero_flag, int npad, int default_width, char *result, 117 size_t *len) 118{ 119 char *put; 120 char *digits; 121 int e, w, d, p, i; 122 char expchar, rchar; 123 format_token ft; 124 /* Number of digits before the decimal point. */ 125 int nbefore; 126 /* Number of zeros after the decimal point. */ 127 int nzero; 128 /* Number of digits after the decimal point. */ 129 int nafter; 130 int leadzero; 131 int nblanks; 132 int ndigits, edigits; 133 sign_t sign; 134 135 ft = f->format; 136 if (f->u.real.w == DEFAULT_WIDTH) 137 /* This codepath can only be reached with -fdec-format-defaults. */ 138 { 139 w = default_width; 140 d = precision; 141 } 142 else 143 { 144 w = f->u.real.w; 145 d = f->u.real.d; 146 } 147 p = dtp->u.p.scale_factor; 148 *len = 0; 149 150 rchar = '5'; 151 152 /* We should always know the field width and precision. */ 153 if (d < 0) 154 internal_error (&dtp->common, "Unspecified precision"); 155 156 sign = calculate_sign (dtp, sign_bit); 157 158 /* Calculate total number of digits. */ 159 if (ft == FMT_F) 160 ndigits = nprinted - 2; 161 else 162 ndigits = precision + 1; 163 164 /* Read the exponent back in. */ 165 if (ft != FMT_F) 166 e = atoi (&buffer[ndigits + 3]) + 1; 167 else 168 e = 0; 169 170 /* Make sure zero comes out as 0.0e0. */ 171 if (zero_flag) 172 e = 0; 173 174 /* Normalize the fractional component. */ 175 if (ft != FMT_F) 176 { 177 buffer[2] = buffer[1]; 178 digits = &buffer[2]; 179 } 180 else 181 digits = &buffer[1]; 182 183 /* Figure out where to place the decimal point. */ 184 switch (ft) 185 { 186 case FMT_F: 187 nbefore = ndigits - precision; 188 if ((w > 0) && (nbefore > (int) size)) 189 { 190 *len = w; 191 star_fill (result, w); 192 result[w] = '\0'; 193 return; 194 } 195 /* Make sure the decimal point is a '.'; depending on the 196 locale, this might not be the case otherwise. */ 197 digits[nbefore] = '.'; 198 if (p != 0) 199 { 200 if (p > 0) 201 { 202 memmove (digits + nbefore, digits + nbefore + 1, p); 203 digits[nbefore + p] = '.'; 204 nbefore += p; 205 nafter = d; 206 nzero = 0; 207 } 208 else /* p < 0 */ 209 { 210 if (nbefore + p >= 0) 211 { 212 nzero = 0; 213 memmove (digits + nbefore + p + 1, digits + nbefore + p, -p); 214 nbefore += p; 215 digits[nbefore] = '.'; 216 nafter = d; 217 } 218 else 219 { 220 nzero = -(nbefore + p); 221 memmove (digits + 1, digits, nbefore); 222 nafter = d - nzero; 223 if (nafter == 0 && d > 0) 224 { 225 /* This is needed to get the correct rounding. */ 226 memmove (digits + 1, digits, ndigits - 1); 227 digits[1] = '0'; 228 nafter = 1; 229 nzero = d - 1; 230 } 231 else if (nafter < 0) 232 { 233 /* Reset digits to 0 in order to get correct rounding 234 towards infinity. */ 235 for (i = 0; i < ndigits; i++) 236 digits[i] = '0'; 237 digits[ndigits - 1] = '1'; 238 nafter = d; 239 nzero = 0; 240 } 241 nbefore = 0; 242 } 243 } 244 } 245 else 246 { 247 nzero = 0; 248 nafter = d; 249 } 250 251 while (digits[0] == '0' && nbefore > 0) 252 { 253 digits++; 254 nbefore--; 255 ndigits--; 256 } 257 258 expchar = 0; 259 /* If we need to do rounding ourselves, get rid of the dot by 260 moving the fractional part. */ 261 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED 262 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED) 263 memmove (digits + nbefore, digits + nbefore + 1, ndigits - nbefore); 264 break; 265 266 case FMT_E: 267 case FMT_D: 268 i = dtp->u.p.scale_factor; 269 if (d < 0 && p == 0) 270 { 271 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not " 272 "greater than zero in format specifier 'E' or 'D'"); 273 return; 274 } 275 if (p <= -d || p >= d + 2) 276 { 277 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor " 278 "out of range in format specifier 'E' or 'D'"); 279 return; 280 } 281 282 if (!zero_flag) 283 e -= p; 284 if (p < 0) 285 { 286 nbefore = 0; 287 nzero = -p; 288 nafter = d + p; 289 } 290 else if (p > 0) 291 { 292 nbefore = p; 293 nzero = 0; 294 nafter = (d - p) + 1; 295 } 296 else /* p == 0 */ 297 { 298 nbefore = 0; 299 nzero = 0; 300 nafter = d; 301 } 302 303 if (ft == FMT_E) 304 expchar = 'E'; 305 else 306 expchar = 'D'; 307 break; 308 309 case FMT_EN: 310 /* The exponent must be a multiple of three, with 1-3 digits before 311 the decimal point. */ 312 if (!zero_flag) 313 e--; 314 if (e >= 0) 315 nbefore = e % 3; 316 else 317 { 318 nbefore = (-e) % 3; 319 if (nbefore != 0) 320 nbefore = 3 - nbefore; 321 } 322 e -= nbefore; 323 nbefore++; 324 nzero = 0; 325 nafter = d; 326 expchar = 'E'; 327 break; 328 329 case FMT_ES: 330 if (!zero_flag) 331 e--; 332 nbefore = 1; 333 nzero = 0; 334 nafter = d; 335 expchar = 'E'; 336 break; 337 338 default: 339 /* Should never happen. */ 340 internal_error (&dtp->common, "Unexpected format token"); 341 } 342 343 if (zero_flag) 344 goto skip; 345 346 /* Round the value. The value being rounded is an unsigned magnitude. */ 347 switch (dtp->u.p.current_unit->round_status) 348 { 349 /* For processor defined and unspecified rounding we use 350 snprintf to print the exact number of digits needed, and thus 351 let snprintf handle the rounding. On system claiming support 352 for IEEE 754, this ought to be round to nearest, ties to 353 even, corresponding to the Fortran ROUND='NEAREST'. */ 354 case ROUND_PROCDEFINED: 355 case ROUND_UNSPECIFIED: 356 case ROUND_ZERO: /* Do nothing and truncation occurs. */ 357 goto skip; 358 case ROUND_UP: 359 if (sign_bit) 360 goto skip; 361 goto updown; 362 case ROUND_DOWN: 363 if (!sign_bit) 364 goto skip; 365 goto updown; 366 case ROUND_NEAREST: 367 /* Round compatible unless there is a tie. A tie is a 5 with 368 all trailing zero's. */ 369 i = nafter + nbefore; 370 if (digits[i] == '5') 371 { 372 for(i++ ; i < ndigits; i++) 373 { 374 if (digits[i] != '0') 375 goto do_rnd; 376 } 377 /* It is a tie so round to even. */ 378 switch (digits[nafter + nbefore - 1]) 379 { 380 case '1': 381 case '3': 382 case '5': 383 case '7': 384 case '9': 385 /* If odd, round away from zero to even. */ 386 break; 387 default: 388 /* If even, skip rounding, truncate to even. */ 389 goto skip; 390 } 391 } 392 /* Fall through. */ 393 /* The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */ 394 case ROUND_COMPATIBLE: 395 rchar = '5'; 396 goto do_rnd; 397 } 398 399 updown: 400 401 rchar = '0'; 402 /* Do not reset nbefore for FMT_F and FMT_EN. */ 403 if (ft != FMT_F && ft !=FMT_EN && w > 0 && d == 0 && p == 0) 404 nbefore = 1; 405 /* Scan for trailing zeros to see if we really need to round it. */ 406 for(i = nbefore + nafter; i < ndigits; i++) 407 { 408 if (digits[i] != '0') 409 goto do_rnd; 410 } 411 goto skip; 412 413 do_rnd: 414 415 if (nbefore + nafter == 0) 416 /* Handle the case Fw.0 and value < 1.0 */ 417 { 418 ndigits = 0; 419 if (digits[0] >= rchar) 420 { 421 /* We rounded to zero but shouldn't have */ 422 nbefore = 1; 423 digits--; 424 digits[0] = '1'; 425 ndigits = 1; 426 } 427 } 428 else if (nbefore + nafter < ndigits) 429 { 430 i = ndigits = nbefore + nafter; 431 if (digits[i] >= rchar) 432 { 433 /* Propagate the carry. */ 434 for (i--; i >= 0; i--) 435 { 436 if (digits[i] != '9') 437 { 438 digits[i]++; 439 break; 440 } 441 digits[i] = '0'; 442 } 443 444 if (i < 0) 445 { 446 /* The carry overflowed. Fortunately we have some spare 447 space at the start of the buffer. We may discard some 448 digits, but this is ok because we already know they are 449 zero. */ 450 digits--; 451 digits[0] = '1'; 452 if (ft == FMT_F) 453 { 454 if (nzero > 0) 455 { 456 nzero--; 457 nafter++; 458 } 459 else 460 nbefore++; 461 } 462 else if (ft == FMT_EN) 463 { 464 nbefore++; 465 if (nbefore == 4) 466 { 467 nbefore = 1; 468 e += 3; 469 } 470 } 471 else 472 e++; 473 } 474 } 475 } 476 477 skip: 478 479 /* Calculate the format of the exponent field. */ 480 if (expchar && !(dtp->u.p.g0_no_blanks && e == 0)) 481 { 482 edigits = 1; 483 for (i = abs (e); i >= 10; i /= 10) 484 edigits++; 485 486 if (f->u.real.e < 0) 487 { 488 /* Width not specified. Must be no more than 3 digits. */ 489 if (e > 999 || e < -999) 490 edigits = -1; 491 else 492 { 493 edigits = 4; 494 if (e > 99 || e < -99) 495 expchar = ' '; 496 } 497 } 498 else if (f->u.real.e == 0) 499 { 500 /* Zero width specified, no leading zeros in exponent */ 501 if (e > 999 || e < -999) 502 edigits = 6; 503 else if (e > 99 || e < -99) 504 edigits = 5; 505 else if (e > 9 || e < -9) 506 edigits = 4; 507 else 508 edigits = 3; 509 } 510 else 511 { 512 /* Exponent width specified, check it is wide enough. */ 513 if (edigits > f->u.real.e) 514 edigits = -1; 515 else 516 edigits = f->u.real.e + 2; 517 } 518 } 519 else 520 edigits = 0; 521 522 /* Scan the digits string and count the number of zeros. If we make it 523 all the way through the loop, we know the value is zero after the 524 rounding completed above. */ 525 int hasdot = 0; 526 for (i = 0; i < ndigits + hasdot; i++) 527 { 528 if (digits[i] == '.') 529 hasdot = 1; 530 else if (digits[i] != '0') 531 break; 532 } 533 534 /* To format properly, we need to know if the rounded result is zero and if 535 so, we set the zero_flag which may have been already set for 536 actual zero. */ 537 if (i == ndigits + hasdot) 538 { 539 zero_flag = true; 540 /* The output is zero, so set the sign according to the sign bit unless 541 -fno-sign-zero was specified. */ 542 if (compile_options.sign_zero == 1) 543 sign = calculate_sign (dtp, sign_bit); 544 else 545 sign = calculate_sign (dtp, 0); 546 } 547 548 /* Pick a field size if none was specified, taking into account small 549 values that may have been rounded to zero. */ 550 if (w <= 0) 551 { 552 if (zero_flag) 553 w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0); 554 else 555 { 556 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1); 557 w = w == 1 ? 2 : w; 558 } 559 } 560 561 /* Work out how much padding is needed. */ 562 nblanks = w - (nbefore + nzero + nafter + edigits + 1); 563 if (sign != S_NONE) 564 nblanks--; 565 566 /* See if we have space for a zero before the decimal point. */ 567 if (nbefore == 0 && nblanks > 0) 568 { 569 leadzero = 1; 570 nblanks--; 571 } 572 else 573 leadzero = 0; 574 575 if (dtp->u.p.g0_no_blanks) 576 { 577 w -= nblanks; 578 nblanks = 0; 579 } 580 581 /* Create the final float string. */ 582 *len = w + npad; 583 put = result; 584 585 /* Check the value fits in the specified field width. */ 586 if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE)) 587 { 588 star_fill (put, *len); 589 return; 590 } 591 592 /* Pad to full field width. */ 593 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank) 594 { 595 memset (put, ' ', nblanks); 596 put += nblanks; 597 } 598 599 /* Set the initial sign (if any). */ 600 if (sign == S_PLUS) 601 *(put++) = '+'; 602 else if (sign == S_MINUS) 603 *(put++) = '-'; 604 605 /* Set an optional leading zero. */ 606 if (leadzero) 607 *(put++) = '0'; 608 609 /* Set the part before the decimal point, padding with zeros. */ 610 if (nbefore > 0) 611 { 612 if (nbefore > ndigits) 613 { 614 i = ndigits; 615 memcpy (put, digits, i); 616 ndigits = 0; 617 while (i < nbefore) 618 put[i++] = '0'; 619 } 620 else 621 { 622 i = nbefore; 623 memcpy (put, digits, i); 624 ndigits -= i; 625 } 626 627 digits += i; 628 put += nbefore; 629 } 630 631 /* Set the decimal point. */ 632 *(put++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ','; 633 if (ft == FMT_F 634 && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED 635 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED)) 636 digits++; 637 638 /* Set leading zeros after the decimal point. */ 639 if (nzero > 0) 640 { 641 for (i = 0; i < nzero; i++) 642 *(put++) = '0'; 643 } 644 645 /* Set digits after the decimal point, padding with zeros. */ 646 if (ndigits >= 0 && nafter > 0) 647 { 648 if (nafter > ndigits) 649 i = ndigits; 650 else 651 i = nafter; 652 653 if (i > 0) 654 memcpy (put, digits, i); 655 while (i < nafter) 656 put[i++] = '0'; 657 658 digits += i; 659 ndigits -= i; 660 put += nafter; 661 } 662 663 /* Set the exponent. */ 664 if (expchar && !(dtp->u.p.g0_no_blanks && e == 0)) 665 { 666 if (expchar != ' ') 667 { 668 *(put++) = expchar; 669 edigits--; 670 } 671 snprintf (buffer, size, "%+0*d", edigits, e); 672 memcpy (put, buffer, edigits); 673 put += edigits; 674 } 675 676 if (dtp->u.p.no_leading_blank) 677 { 678 memset (put , ' ' , nblanks); 679 dtp->u.p.no_leading_blank = 0; 680 put += nblanks; 681 } 682 683 if (npad > 0 && !dtp->u.p.g0_no_blanks) 684 { 685 memset (put , ' ' , npad); 686 put += npad; 687 } 688 689 /* NULL terminate the string. */ 690 *put = '\0'; 691 692 return; 693} 694 695 696/* Write "Infinite" or "Nan" as appropriate for the given format. */ 697 698static void 699build_infnan_string (st_parameter_dt *dtp, const fnode *f, int isnan_flag, 700 int sign_bit, char *p, size_t *len) 701{ 702 char fin; 703 int nb = 0; 704 sign_t sign; 705 int mark; 706 707 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z) 708 { 709 sign = calculate_sign (dtp, sign_bit); 710 mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7; 711 712 nb = f->u.real.w; 713 *len = nb; 714 715 /* If the field width is zero, the processor must select a width 716 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */ 717 718 if ((nb == 0) || dtp->u.p.g0_no_blanks) 719 { 720 if (isnan_flag) 721 nb = 3; 722 else 723 nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3; 724 *len = nb; 725 } 726 727 p[*len] = '\0'; 728 if (nb < 3) 729 { 730 memset (p, '*', nb); 731 return; 732 } 733 734 memset(p, ' ', nb); 735 736 if (!isnan_flag) 737 { 738 if (sign_bit) 739 { 740 /* If the sign is negative and the width is 3, there is 741 insufficient room to output '-Inf', so output asterisks */ 742 if (nb == 3) 743 { 744 memset (p, '*', nb); 745 return; 746 } 747 /* The negative sign is mandatory */ 748 fin = '-'; 749 } 750 else 751 /* The positive sign is optional, but we output it for 752 consistency */ 753 fin = '+'; 754 755 if (nb > mark) 756 /* We have room, so output 'Infinity' */ 757 memcpy(p + nb - 8, "Infinity", 8); 758 else 759 /* For the case of width equals 8, there is not enough room 760 for the sign and 'Infinity' so we go with 'Inf' */ 761 memcpy(p + nb - 3, "Inf", 3); 762 763 if (sign == S_PLUS || sign == S_MINUS) 764 { 765 if (nb < 9 && nb > 3) 766 p[nb - 4] = fin; /* Put the sign in front of Inf */ 767 else if (nb > 8) 768 p[nb - 9] = fin; /* Put the sign in front of Infinity */ 769 } 770 } 771 else 772 memcpy(p + nb - 3, "NaN", 3); 773 } 774} 775 776 777/* Returns the value of 10**d. */ 778 779#define CALCULATE_EXP(x) \ 780static GFC_REAL_ ## x \ 781calculate_exp_ ## x (int d)\ 782{\ 783 int i;\ 784 GFC_REAL_ ## x r = 1.0;\ 785 for (i = 0; i< (d >= 0 ? d : -d); i++)\ 786 r *= 10;\ 787 r = (d >= 0) ? r : 1.0 / r;\ 788 return r;\ 789} 790 791CALCULATE_EXP(4) 792 793CALCULATE_EXP(8) 794 795#ifdef HAVE_GFC_REAL_10 796CALCULATE_EXP(10) 797#endif 798 799#ifdef HAVE_GFC_REAL_16 800CALCULATE_EXP(16) 801#endif 802#undef CALCULATE_EXP 803 804 805/* Define macros to build code for format_float. */ 806 807 /* Note: Before output_float is called, snprintf is used to print to buffer the 808 number in the format +D.DDDDe+ddd. 809 810 # The result will always contain a decimal point, even if no 811 digits follow it 812 813 - The converted value is to be left adjusted on the field boundary 814 815 + A sign (+ or -) always be placed before a number 816 817 * prec is used as the precision 818 819 e format: [-]d.ddde±dd where there is one digit before the 820 decimal-point character and the number of digits after it is 821 equal to the precision. The exponent always contains at least two 822 digits; if the value is zero, the exponent is 00. */ 823 824 825#define TOKENPASTE(x, y) TOKENPASTE2(x, y) 826#define TOKENPASTE2(x, y) x ## y 827 828#define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val) 829 830#define DTOA2(prec,val) \ 831snprintf (buffer, size, "%+-#.*e", (prec), (val)) 832 833#define DTOA2L(prec,val) \ 834snprintf (buffer, size, "%+-#.*Le", (prec), (val)) 835 836 837#if defined(GFC_REAL_16_IS_FLOAT128) 838#define DTOA2Q(prec,val) \ 839quadmath_snprintf (buffer, size, "%+-#.*Qe", (prec), (val)) 840#endif 841 842#define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val) 843 844/* For F format, we print to the buffer with f format. */ 845#define FDTOA2(prec,val) \ 846snprintf (buffer, size, "%+-#.*f", (prec), (val)) 847 848#define FDTOA2L(prec,val) \ 849snprintf (buffer, size, "%+-#.*Lf", (prec), (val)) 850 851 852#if defined(GFC_REAL_16_IS_FLOAT128) 853#define FDTOA2Q(prec,val) \ 854quadmath_snprintf (buffer, size, "%+-#.*Qf", \ 855 (prec), (val)) 856#endif 857 858 859/* EN format is tricky since the number of significant digits depends 860 on the magnitude. Solve it by first printing a temporary value and 861 figure out the number of significant digits from the printed 862 exponent. Values y, 0.95*10.0**e <= y <10.0**e, are rounded to 863 10.0**e even when the final result will not be rounded to 10.0**e. 864 For these values the exponent returned by atoi has to be decremented 865 by one. The values y in the ranges 866 (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1)) 867 (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2) 868 (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1) 869 are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)), 870 100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0 871 represents d zeroes, by the lines 279 to 297. */ 872#define EN_PREC(x,y)\ 873{\ 874 volatile GFC_REAL_ ## x tmp, one = 1.0;\ 875 tmp = * (GFC_REAL_ ## x *)source;\ 876 if (isfinite (tmp))\ 877 {\ 878 nprinted = DTOA(y,0,tmp);\ 879 int e = atoi (&buffer[4]);\ 880 if (buffer[1] == '1')\ 881 {\ 882 tmp = (calculate_exp_ ## x (-e)) * tmp;\ 883 tmp = one - (tmp < 0 ? -tmp : tmp);\ 884 if (tmp > 0)\ 885 e = e - 1;\ 886 }\ 887 nbefore = e%3;\ 888 if (nbefore < 0)\ 889 nbefore = 3 + nbefore;\ 890 }\ 891 else\ 892 nprinted = -1;\ 893}\ 894 895static int 896determine_en_precision (st_parameter_dt *dtp, const fnode *f, 897 const char *source, int len) 898{ 899 int nprinted; 900 char buffer[10]; 901 const size_t size = 10; 902 int nbefore; /* digits before decimal point - 1. */ 903 904 switch (len) 905 { 906 case 4: 907 EN_PREC(4,) 908 break; 909 910 case 8: 911 EN_PREC(8,) 912 break; 913 914#ifdef HAVE_GFC_REAL_10 915 case 10: 916 EN_PREC(10,L) 917 break; 918#endif 919#ifdef HAVE_GFC_REAL_16 920 case 16: 921# ifdef GFC_REAL_16_IS_FLOAT128 922 EN_PREC(16,Q) 923# else 924 EN_PREC(16,L) 925# endif 926 break; 927#endif 928 default: 929 internal_error (NULL, "bad real kind"); 930 } 931 932 if (nprinted == -1) 933 return -1; 934 935 int prec = f->u.real.d + nbefore; 936 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED 937 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED) 938 prec += 2 * len + 4; 939 return prec; 940} 941 942 943/* Generate corresponding I/O format. and output. 944 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran 945 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is: 946 947 Data Magnitude Equivalent Conversion 948 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee] 949 m = 0 F(w-n).(d-1), n' ' 950 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' ' 951 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' ' 952 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' ' 953 ................ .......... 954 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ') 955 m >= 10**d-0.5 Ew.d[Ee] 956 957 notes: for Gw.d , n' ' means 4 blanks 958 for Gw.dEe, n' ' means e+2 blanks 959 for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2 960 the asm volatile is required for 32-bit x86 platforms. */ 961#define FORMAT_FLOAT(x,y)\ 962{\ 963 int npad = 0;\ 964 GFC_REAL_ ## x m;\ 965 m = * (GFC_REAL_ ## x *)source;\ 966 sign_bit = signbit (m);\ 967 if (!isfinite (m))\ 968 { \ 969 build_infnan_string (dtp, f, isnan (m), sign_bit, result, res_len);\ 970 return;\ 971 }\ 972 m = sign_bit ? -m : m;\ 973 zero_flag = (m == 0.0);\ 974 if (f->format == FMT_G)\ 975 {\ 976 int e = f->u.real.e;\ 977 int d = f->u.real.d;\ 978 int w = f->u.real.w;\ 979 fnode newf;\ 980 GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\ 981 int low, high, mid;\ 982 int ubound, lbound;\ 983 int save_scale_factor;\ 984 volatile GFC_REAL_ ## x temp;\ 985 save_scale_factor = dtp->u.p.scale_factor;\ 986 if (w == DEFAULT_WIDTH)\ 987 {\ 988 w = default_width;\ 989 d = precision;\ 990 }\ 991 /* The switch between FMT_E and FMT_F is based on the absolute value. \ 992 Set r=0 for rounding toward zero and r = 1 otherwise. \ 993 If (exp_d - m) == 1 there is no rounding needed. */\ 994 switch (dtp->u.p.current_unit->round_status)\ 995 {\ 996 case ROUND_ZERO:\ 997 r = 0.0;\ 998 break;\ 999 case ROUND_UP:\ 1000 r = sign_bit ? 0.0 : 1.0;\ 1001 break;\ 1002 case ROUND_DOWN:\ 1003 r = sign_bit ? 1.0 : 0.0;\ 1004 break;\ 1005 default:\ 1006 break;\ 1007 }\ 1008 exp_d = calculate_exp_ ## x (d);\ 1009 r_sc = (1 - r / exp_d);\ 1010 temp = 0.1 * r_sc;\ 1011 if ((m > 0.0 && ((m < temp) || (r < 1 && r >= (exp_d - m))\ 1012 || (r == 1 && 1 > (exp_d - m))))\ 1013 || ((m == 0.0) && !(compile_options.allow_std\ 1014 & (GFC_STD_F2003 | GFC_STD_F2008)))\ 1015 || d == 0)\ 1016 { \ 1017 newf.format = FMT_E;\ 1018 newf.u.real.w = w;\ 1019 newf.u.real.d = d - comp_d;\ 1020 newf.u.real.e = e;\ 1021 npad = 0;\ 1022 precision = determine_precision (dtp, &newf, x);\ 1023 nprinted = DTOA(y,precision,m);\ 1024 }\ 1025 else \ 1026 {\ 1027 mid = 0;\ 1028 low = 0;\ 1029 high = d + 1;\ 1030 lbound = 0;\ 1031 ubound = d + 1;\ 1032 while (low <= high)\ 1033 {\ 1034 mid = (low + high) / 2;\ 1035 temp = (calculate_exp_ ## x (mid - 1) * r_sc);\ 1036 if (m < temp)\ 1037 { \ 1038 ubound = mid;\ 1039 if (ubound == lbound + 1)\ 1040 break;\ 1041 high = mid - 1;\ 1042 }\ 1043 else if (m > temp)\ 1044 { \ 1045 lbound = mid;\ 1046 if (ubound == lbound + 1)\ 1047 { \ 1048 mid ++;\ 1049 break;\ 1050 }\ 1051 low = mid + 1;\ 1052 }\ 1053 else\ 1054 {\ 1055 mid++;\ 1056 break;\ 1057 }\ 1058 }\ 1059 npad = e <= 0 ? 4 : e + 2;\ 1060 npad = npad >= w ? w - 1 : npad;\ 1061 npad = dtp->u.p.g0_no_blanks ? 0 : npad;\ 1062 newf.format = FMT_F;\ 1063 newf.u.real.w = w - npad;\ 1064 newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\ 1065 dtp->u.p.scale_factor = 0;\ 1066 precision = determine_precision (dtp, &newf, x);\ 1067 nprinted = FDTOA(y,precision,m);\ 1068 }\ 1069 build_float_string (dtp, &newf, buffer, size, nprinted, precision,\ 1070 sign_bit, zero_flag, npad, default_width,\ 1071 result, res_len);\ 1072 dtp->u.p.scale_factor = save_scale_factor;\ 1073 }\ 1074 else\ 1075 {\ 1076 if (f->format == FMT_F)\ 1077 nprinted = FDTOA(y,precision,m);\ 1078 else\ 1079 nprinted = DTOA(y,precision,m);\ 1080 build_float_string (dtp, f, buffer, size, nprinted, precision,\ 1081 sign_bit, zero_flag, npad, default_width,\ 1082 result, res_len);\ 1083 }\ 1084}\ 1085 1086/* Output a real number according to its format. */ 1087 1088 1089static void 1090get_float_string (st_parameter_dt *dtp, const fnode *f, const char *source, 1091 int kind, int comp_d, char *buffer, int precision, 1092 size_t size, char *result, size_t *res_len) 1093{ 1094 int sign_bit, nprinted; 1095 bool zero_flag; 1096 int default_width = 0; 1097 1098 if (f->u.real.w == DEFAULT_WIDTH) 1099 /* This codepath can only be reached with -fdec-format-defaults. The default 1100 * values are based on those used in the Oracle Fortran compiler. 1101 */ 1102 { 1103 default_width = default_width_for_float (kind); 1104 precision = default_precision_for_float (kind); 1105 } 1106 1107 switch (kind) 1108 { 1109 case 4: 1110 FORMAT_FLOAT(4,) 1111 break; 1112 1113 case 8: 1114 FORMAT_FLOAT(8,) 1115 break; 1116 1117#ifdef HAVE_GFC_REAL_10 1118 case 10: 1119 FORMAT_FLOAT(10,L) 1120 break; 1121#endif 1122#ifdef HAVE_GFC_REAL_16 1123 case 16: 1124# ifdef GFC_REAL_16_IS_FLOAT128 1125 FORMAT_FLOAT(16,Q) 1126# else 1127 FORMAT_FLOAT(16,L) 1128# endif 1129 break; 1130#endif 1131 default: 1132 internal_error (NULL, "bad real kind"); 1133 } 1134 return; 1135} 1136