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