libgfortran/io/write_float.def

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