libgfortran/io/write_float.def

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