libgfortran/io/write_float.def

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