libgfortran/io/write_float.def

   1 /* Copyright (C) 2007-2014 Free Software Foundation, Inc.
   2    Contributed by Andy Vaught
   3    Write float code factoring to this file by Jerry DeLisle
   4    F2003 I/O support contributed by Jerry DeLisle
   5
   6 This file is part of the GNU Fortran runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU General Public License as published by
  10 the Free Software Foundation; either version 3, or (at your option)
  11 any later version.
  12
  13 Libgfortran is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27 #include "config.h"
  28
  29 typedef enum
  30 { S_NONE, S_MINUS, S_PLUS }
  31 sign_t;
  32
  33 /* Given a flag that indicates if a value is negative or not, return a
  34    sign_t that gives the sign that we need to produce.  */
  35
  36 static sign_t
  37 calculate_sign (st_parameter_dt *dtp, int negative_flag)
  38 {
  39   sign_t s = S_NONE;
  40
  41   if (negative_flag)
  42     s = S_MINUS;
  43   else
  44     switch (dtp->u.p.sign_status)
  45       {
  46       case SIGN_SP:     /* Show sign. */
  47         s = S_PLUS;
  48         break;
  49       case SIGN_SS:     /* Suppress sign. */
  50         s = S_NONE;
  51         break;
  52       case SIGN_S:      /* Processor defined. */
  53       case SIGN_UNSPECIFIED:
  54         s = options.optional_plus ? S_PLUS : S_NONE;
  55         break;
  56       }
  57
  58   return s;
  59 }
  60
  61
  62 /* Determine the precision except for EN format. For G format,
  63    determines an upper bound to be used for sizing the buffer. */
  64
  65 static int
  66 determine_precision (st_parameter_dt * dtp, const fnode * f, int len)
  67 {
  68   int precision = f->u.real.d;
  69
  70   switch (f->format)
  71     {
  72     case FMT_F:
  73     case FMT_G:
  74       precision += dtp->u.p.scale_factor;
  75       break;
  76     case FMT_ES:
  77       /* Scale factor has no effect on output.  */
  78       break;
  79     case FMT_E:
  80     case FMT_D:
  81       /* See F2008 10.7.2.3.3.6 */
  82       if (dtp->u.p.scale_factor <= 0)
  83         precision += dtp->u.p.scale_factor - 1;
  84       break;
  85     default:
  86       return -1;
  87     }
  88
  89   /* If the scale factor has a large negative value, we must do our
  90      own rounding? Use ROUND='NEAREST', which should be what snprintf
  91      is using as well.  */
  92   if (precision < 0 &&
  93       (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
  94        || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
  95     dtp->u.p.current_unit->round_status = ROUND_NEAREST;
  96
  97   /* Add extra guard digits up to at least full precision when we do
  98      our own rounding.  */
  99   if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
 100       && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
 101     {
 102       precision += 2 * len + 4;
 103       if (precision < 0)
 104         precision = 0;
 105     }
 106
 107   return precision;
 108 }
 109
 110
 111 /* Output a real number according to its format which is FMT_G free.  */
 112
 113 static bool
 114 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
 115               int nprinted, int precision, int sign_bit, bool zero_flag)
 116 {
 117   char *out;
 118   char *digits;
 119   int e, w, d, p, i;
 120   char expchar, rchar;
 121   format_token ft;
 122   /* Number of digits before the decimal point.  */
 123   int nbefore;
 124   /* Number of zeros after the decimal point.  */
 125   int nzero;
 126   /* Number of digits after the decimal point.  */
 127   int nafter;
 128   /* 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 false;
 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 false;
 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 (ft != FMT_F && nbefore == 0 && 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     /* Handle the case Fw.0 and value < 1.0 */
 390     {
 391       ndigits = 0;
 392       if (nzero_real == d && digits[0] >= rchar)
 393         {
 394           /* We rounded to zero but shouldn't have */
 395           nbefore = 1;
 396           digits--;
 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 false;
 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 false;
 546         }
 547       star_fill (out, w);
 548       return false;
 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 true;
 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 true;
 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 #if defined(GFC_REAL_16_IS_FLOAT128)
 966 #define ISFINITE2Q(val) finiteq(val)
 967 #endif
 968 #define ISFINITE2(val) isfinite(val)
 969 #define ISFINITE2L(val) isfinite(val)
 970
 971 #define ISFINITE(suff,val) TOKENPASTE(ISFINITE2,suff)(val)
 972
 973
 974 #if defined(GFC_REAL_16_IS_FLOAT128)
 975 #define SIGNBIT2Q(val) signbitq(val)
 976 #endif
 977 #define SIGNBIT2(val) signbit(val)
 978 #define SIGNBIT2L(val) signbit(val)
 979
 980 #define SIGNBIT(suff,val) TOKENPASTE(SIGNBIT2,suff)(val)
 981
 982
 983 #if defined(GFC_REAL_16_IS_FLOAT128)
 984 #define ISNAN2Q(val) isnanq(val)
 985 #endif
 986 #define ISNAN2(val) isnan(val)
 987 #define ISNAN2L(val) isnan(val)
 988
 989 #define ISNAN(suff,val) TOKENPASTE(ISNAN2,suff)(val)
 990
 991
 992
 993 /* Generate corresponding I/O format for FMT_G and output.
 994    The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
 995    LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
 996
 997    Data Magnitude                              Equivalent Conversion
 998    0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
 999    m = 0                                       F(w-n).(d-1), n' '
1000    0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
1001    1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
1002    10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
1003    ................                           ..........
1004    10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
1005    m >= 10**d-0.5                              Ew.d[Ee]
1006
1007    notes: for Gw.d ,  n' ' means 4 blanks
1008           for Gw.dEe, n' ' means e+2 blanks
1009           for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
1010           the asm volatile is required for 32-bit x86 platforms.  */
1011
1012 #define OUTPUT_FLOAT_FMT_G(x,y)                 \
1013 static void \
1014 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
1015                       GFC_REAL_ ## x m, char *buffer, size_t size, \
1016                           int sign_bit, bool zero_flag, int comp_d) \
1017 { \
1018   int e = f->u.real.e;\
1019   int d = f->u.real.d;\
1020   int w = f->u.real.w;\
1021   fnode newf;\
1022   GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
1023   int low, high, mid;\
1024   int ubound, lbound;\
1025   char *p, pad = ' ';\
1026   int save_scale_factor, nb = 0;\
1027   bool result;\
1028   int nprinted, precision;\
1029   volatile GFC_REAL_ ## x temp;\
1030 \
1031   save_scale_factor = dtp->u.p.scale_factor;\
1032 \
1033   switch (dtp->u.p.current_unit->round_status)\
1034     {\
1035       case ROUND_ZERO:\
1036         r = sign_bit ? 1.0 : 0.0;\
1037         break;\
1038       case ROUND_UP:\
1039         r = 1.0;\
1040         break;\
1041       case ROUND_DOWN:\
1042         r = 0.0;\
1043         break;\
1044       default:\
1045         break;\
1046     }\
1047 \
1048   exp_d = calculate_exp_ ## x (d);\
1049   r_sc = (1 - r / exp_d);\
1050   temp = 0.1 * r_sc;\
1051   if ((m > 0.0 && ((m < temp) || (r >= (exp_d - m))))\
1052       || ((m == 0.0) && !(compile_options.allow_std\
1053                           & (GFC_STD_F2003 | GFC_STD_F2008)))\
1054       ||  d == 0)\
1055     { \
1056       newf.format = FMT_E;\
1057       newf.u.real.w = w;\
1058       newf.u.real.d = d - comp_d;\
1059       newf.u.real.e = e;\
1060       nb = 0;\
1061       precision = determine_precision (dtp, &newf, x);\
1062       nprinted = DTOA(y,precision,m);                        \
1063       goto finish;\
1064     }\
1065 \
1066   mid = 0;\
1067   low = 0;\
1068   high = d + 1;\
1069   lbound = 0;\
1070   ubound = d + 1;\
1071 \
1072   while (low <= high)\
1073     { \
1074       mid = (low + high) / 2;\
1075 \
1076       temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
1077 \
1078       if (m < temp)\
1079         { \
1080           ubound = mid;\
1081           if (ubound == lbound + 1)\
1082             break;\
1083           high = mid - 1;\
1084         }\
1085       else if (m > temp)\
1086         { \
1087           lbound = mid;\
1088           if (ubound == lbound + 1)\
1089             { \
1090               mid ++;\
1091               break;\
1092             }\
1093           low = mid + 1;\
1094         }\
1095       else\
1096         {\
1097           mid++;\
1098           break;\
1099         }\
1100     }\
1101 \
1102   nb = e <= 0 ? 4 : e + 2;\
1103   nb = nb >= w ? w - 1 : nb;\
1104   newf.format = FMT_F;\
1105   newf.u.real.w = w - nb;\
1106   newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
1107   dtp->u.p.scale_factor = 0;\
1108   precision = determine_precision (dtp, &newf, x);      \
1109   nprinted = FDTOA(y,precision,m);                                      \
1110 \
1111  finish:\
1112     result = output_float (dtp, &newf, buffer, size, nprinted, precision,\
1113                            sign_bit, zero_flag);\
1114   dtp->u.p.scale_factor = save_scale_factor;\
1115 \
1116 \
1117   if (nb > 0 && !dtp->u.p.g0_no_blanks)\
1118     {\
1119       p = write_block (dtp, nb);\
1120       if (p == NULL)\
1121         return;\
1122       if (!result)\
1123         pad = '*';\
1124       if (unlikely (is_char4_unit (dtp)))\
1125         {\
1126           gfc_char4_t *p4 = (gfc_char4_t *) p;\
1127           memset4 (p4, pad, nb);\
1128         }\
1129       else \
1130         memset (p, pad, nb);\
1131     }\
1132 }\
1133
1134 OUTPUT_FLOAT_FMT_G(4,)
1135
1136 OUTPUT_FLOAT_FMT_G(8,)
1137
1138 #ifdef HAVE_GFC_REAL_10
1139 OUTPUT_FLOAT_FMT_G(10,L)
1140 #endif
1141
1142 #ifdef HAVE_GFC_REAL_16
1143 # ifdef GFC_REAL_16_IS_FLOAT128
1144 OUTPUT_FLOAT_FMT_G(16,Q)
1145 #else
1146 OUTPUT_FLOAT_FMT_G(16,L)
1147 #endif
1148 #endif
1149
1150 #undef OUTPUT_FLOAT_FMT_G
1151
1152
1153 /* EN format is tricky since the number of significant digits depends
1154    on the magnitude.  Solve it by first printing a temporary value and
1155    figure out the number of significant digits from the printed
1156    exponent.  */
1157
1158 #define EN_PREC(x,y)\
1159 {\
1160     GFC_REAL_ ## x tmp;                         \
1161     tmp = * (GFC_REAL_ ## x *)source;                           \
1162     if (ISFINITE (y,tmp))                                       \
1163       nprinted = DTOA(y,0,tmp);                                 \
1164     else\
1165       nprinted = -1;\
1166 }\
1167
1168 static int
1169 determine_en_precision (st_parameter_dt *dtp, const fnode *f,
1170                         const char *source, int len)
1171 {
1172   int nprinted;
1173   char buffer[10];
1174   const size_t size = 10;
1175
1176   switch (len)
1177     {
1178     case 4:
1179       EN_PREC(4,)
1180       break;
1181
1182     case 8:
1183       EN_PREC(8,)
1184       break;
1185
1186 #ifdef HAVE_GFC_REAL_10
1187     case 10:
1188       EN_PREC(10,L)
1189       break;
1190 #endif
1191 #ifdef HAVE_GFC_REAL_16
1192     case 16:
1193 # ifdef GFC_REAL_16_IS_FLOAT128
1194       EN_PREC(16,Q)
1195 # else
1196       EN_PREC(16,L)
1197 # endif
1198       break;
1199 #endif
1200     default:
1201       internal_error (NULL, "bad real kind");
1202     }
1203
1204   if (nprinted == -1)
1205     return -1;
1206
1207   int e = atoi (&buffer[5]);
1208   int nbefore; /* digits before decimal point - 1.  */
1209   if (e >= 0)
1210     nbefore = e % 3;
1211   else
1212     {
1213       nbefore = (-e) % 3;
1214       if (nbefore != 0)
1215         nbefore = 3 - nbefore;
1216     }
1217   int prec = f->u.real.d + nbefore;
1218   if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
1219       && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
1220     prec += 2 * len + 4;
1221   return prec;
1222 }
1223
1224
1225 #define WRITE_FLOAT(x,y)\
1226 {\
1227         GFC_REAL_ ## x tmp;\
1228         tmp = * (GFC_REAL_ ## x *)source;\
1229         sign_bit = SIGNBIT (y,tmp);\
1230         if (!ISFINITE (y,tmp))\
1231           { \
1232             write_infnan (dtp, f, ISNAN (y,tmp), sign_bit);\
1233             return;\
1234           }\
1235         tmp = sign_bit ? -tmp : tmp;\
1236         zero_flag = (tmp == 0.0);\
1237         if (f->format == FMT_G)\
1238           output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
1239                                     zero_flag, comp_d);\
1240         else\
1241           {\
1242             if (f->format == FMT_F)\
1243               nprinted = FDTOA(y,precision,tmp);                \
1244             else\
1245               nprinted = DTOA(y,precision,tmp);                                 \
1246             output_float (dtp, f, buffer, size, nprinted, precision,\
1247                           sign_bit, zero_flag);\
1248           }\
1249 }\
1250
1251 /* Output a real number according to its format.  */
1252
1253 static void
1254 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
1255             int len, int comp_d)
1256 {
1257   int sign_bit, nprinted;
1258   int precision;  /* Precision for snprintf call.  */
1259   bool zero_flag;
1260
1261   if (f->format != FMT_EN)
1262     precision = determine_precision (dtp, f, len);
1263   else
1264     precision = determine_en_precision (dtp, f, source, len);
1265
1266   /* 4932 is the maximum exponent of long double and quad precision, 3
1267      extra characters for the sign, the decimal point, and the
1268      trailing null, and finally some extra digits depending on the
1269      requested precision.  */
1270   const size_t size = 4932 + 3 + precision;
1271   char buffer[size];
1272
1273   switch (len)
1274     {
1275     case 4:
1276       WRITE_FLOAT(4,)
1277       break;
1278
1279     case 8:
1280       WRITE_FLOAT(8,)
1281       break;
1282
1283 #ifdef HAVE_GFC_REAL_10
1284     case 10:
1285       WRITE_FLOAT(10,L)
1286       break;
1287 #endif
1288 #ifdef HAVE_GFC_REAL_16
1289     case 16:
1290 # ifdef GFC_REAL_16_IS_FLOAT128
1291       WRITE_FLOAT(16,Q)
1292 # else
1293       WRITE_FLOAT(16,L)
1294 # endif
1295       break;
1296 #endif
1297     default:
1298       internal_error (NULL, "bad real kind");
1299     }
1300 }