libgfortran/io/write_float.def

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