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