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       ndigits = nbefore + nafter;
 293       i = ndigits;
 294       if (digits[i] >= rchar)
 295         {
 296           /* Propagate the carry.  */
 297           for (i--; i >= 0; i--)
 298             {
 299               if (digits[i] != '9')
 300                 {
 301                   digits[i]++;
 302                   break;
 303                 }
 304               digits[i] = '0';
 305             }
 306
 307           if (i < 0)
 308             {
 309               /* The carry overflowed.  Fortunately we have some spare
 310                  space at the start of the buffer.  We may discard some
 311                  digits, but this is ok because we already know they are
 312                  zero.  */
 313               digits--;
 314               digits[0] = '1';
 315               if (ft == FMT_F)
 316                 {
 317                   if (nzero > 0)
 318                     {
 319                       nzero--;
 320                       nafter++;
 321                     }
 322                   else
 323                     nbefore++;
 324                 }
 325               else if (ft == FMT_EN)
 326                 {
 327                   nbefore++;
 328                   if (nbefore == 4)
 329                     {
 330                       nbefore = 1;
 331                       e += 3;
 332                     }
 333                 }
 334               else
 335                 e++;
 336             }
 337         }
 338     }
 339
 340   skip:
 341
 342   /* Calculate the format of the exponent field.  */
 343   if (expchar)
 344     {
 345       edigits = 1;
 346       for (i = abs (e); i >= 10; i /= 10)
 347         edigits++;
 348
 349       if (f->u.real.e < 0)
 350         {
 351           /* Width not specified.  Must be no more than 3 digits.  */
 352           if (e > 999 || e < -999)
 353             edigits = -1;
 354           else
 355             {
 356               edigits = 4;
 357               if (e > 99 || e < -99)
 358                 expchar = ' ';
 359             }
 360         }
 361       else
 362         {
 363           /* Exponent width specified, check it is wide enough.  */
 364           if (edigits > f->u.real.e)
 365             edigits = -1;
 366           else
 367             edigits = f->u.real.e + 2;
 368         }
 369     }
 370   else
 371     edigits = 0;
 372
 373   /* Scan the digits string and count the number of zeros.  If we make it
 374      all the way through the loop, we know the value is zero after the
 375      rounding completed above.  */
 376   for (i = 0; i < ndigits; i++)
 377     {
 378       if (digits[i] != '0')
 379         break;
 380     }
 381
 382   /* To format properly, we need to know if the rounded result is zero and if
 383      so, we set the zero_flag which may have been already set for
 384      actual zero.  */
 385   if (i == ndigits)
 386     {
 387       zero_flag = true;
 388       /* The output is zero, so set the sign according to the sign bit unless
 389          -fno-sign-zero was specified.  */
 390       if (compile_options.sign_zero == 1)
 391         sign = calculate_sign (dtp, sign_bit);
 392       else
 393         sign = calculate_sign (dtp, 0);
 394     }
 395
 396   /* Pick a field size if none was specified, taking into account small
 397      values that may have been rounded to zero.  */
 398   if (w <= 0)
 399     {
 400       if (zero_flag)
 401         w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
 402       else
 403         {
 404           w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
 405           w = w == 1 ? 2 : w;
 406         }
 407     }
 408
 409   /* Work out how much padding is needed.  */
 410   nblanks = w - (nbefore + nzero + nafter + edigits + 1);
 411   if (sign != S_NONE)
 412     nblanks--;
 413
 414   if (dtp->u.p.g0_no_blanks)
 415     {
 416       w -= nblanks;
 417       nblanks = 0;
 418     }
 419
 420   /* Create the ouput buffer.  */
 421   out = write_block (dtp, w);
 422   if (out == NULL)
 423     return FAILURE;
 424
 425   /* Check the value fits in the specified field width.  */
 426   if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
 427     {
 428       if (unlikely (is_char4_unit (dtp)))
 429         {
 430           gfc_char4_t *out4 = (gfc_char4_t *) out;
 431           memset4 (out4, '*', w);
 432           return FAILURE;
 433         }
 434       star_fill (out, w);
 435       return FAILURE;
 436     }
 437
 438   /* See if we have space for a zero before the decimal point.  */
 439   if (nbefore == 0 && nblanks > 0)
 440     {
 441       leadzero = 1;
 442       nblanks--;
 443     }
 444   else
 445     leadzero = 0;
 446
 447   /* For internal character(kind=4) units, we duplicate the code used for
 448      regular output slightly modified.  This needs to be maintained
 449      consistent with the regular code that follows this block.  */
 450   if (unlikely (is_char4_unit (dtp)))
 451     {
 452       gfc_char4_t *out4 = (gfc_char4_t *) out;
 453       /* Pad to full field width.  */
 454
 455       if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
 456         {
 457           memset4 (out4, ' ', nblanks);
 458           out4 += nblanks;
 459         }
 460
 461       /* Output the initial sign (if any).  */
 462       if (sign == S_PLUS)
 463         *(out4++) = '+';
 464       else if (sign == S_MINUS)
 465         *(out4++) = '-';
 466
 467       /* Output an optional leading zero.  */
 468       if (leadzero)
 469         *(out4++) = '0';
 470
 471       /* Output the part before the decimal point, padding with zeros.  */
 472       if (nbefore > 0)
 473         {
 474           if (nbefore > ndigits)
 475             {
 476               i = ndigits;
 477               memcpy4 (out4, digits, i);
 478               ndigits = 0;
 479               while (i < nbefore)
 480                 out4[i++] = '0';
 481             }
 482           else
 483             {
 484               i = nbefore;
 485               memcpy4 (out4, digits, i);
 486               ndigits -= i;
 487             }
 488
 489           digits += i;
 490           out4 += nbefore;
 491         }
 492
 493       /* Output the decimal point.  */
 494       *(out4++) = dtp->u.p.current_unit->decimal_status
 495                     == DECIMAL_POINT ? '.' : ',';
 496
 497       /* Output leading zeros after the decimal point.  */
 498       if (nzero > 0)
 499         {
 500           for (i = 0; i < nzero; i++)
 501             *(out4++) = '0';
 502         }
 503
 504       /* Output digits after the decimal point, padding with zeros.  */
 505       if (nafter > 0)
 506         {
 507           if (nafter > ndigits)
 508             i = ndigits;
 509           else
 510             i = nafter;
 511
 512           memcpy4 (out4, digits, i);
 513           while (i < nafter)
 514             out4[i++] = '0';
 515
 516           digits += i;
 517           ndigits -= i;
 518           out4 += nafter;
 519         }
 520
 521       /* Output the exponent.  */
 522       if (expchar)
 523         {
 524           if (expchar != ' ')
 525             {
 526               *(out4++) = expchar;
 527               edigits--;
 528             }
 529           snprintf (buffer, size, "%+0*d", edigits, e);
 530           memcpy4 (out4, buffer, edigits);
 531         }
 532
 533       if (dtp->u.p.no_leading_blank)
 534         {
 535           out4 += edigits;
 536           memset4 (out4, ' ' , nblanks);
 537           dtp->u.p.no_leading_blank = 0;
 538         }
 539       return SUCCESS;
 540     } /* End of character(kind=4) internal unit code.  */
 541
 542   /* Pad to full field width.  */
 543
 544   if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
 545     {
 546       memset (out, ' ', nblanks);
 547       out += nblanks;
 548     }
 549
 550   /* Output the initial sign (if any).  */
 551   if (sign == S_PLUS)
 552     *(out++) = '+';
 553   else if (sign == S_MINUS)
 554     *(out++) = '-';
 555
 556   /* Output an optional leading zero.  */
 557   if (leadzero)
 558     *(out++) = '0';
 559
 560   /* Output the part before the decimal point, padding with zeros.  */
 561   if (nbefore > 0)
 562     {
 563       if (nbefore > ndigits)
 564         {
 565           i = ndigits;
 566           memcpy (out, digits, i);
 567           ndigits = 0;
 568           while (i < nbefore)
 569             out[i++] = '0';
 570         }
 571       else
 572         {
 573           i = nbefore;
 574           memcpy (out, digits, i);
 575           ndigits -= i;
 576         }
 577
 578       digits += i;
 579       out += nbefore;
 580     }
 581
 582   /* Output the decimal point.  */
 583   *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
 584
 585   /* Output leading zeros after the decimal point.  */
 586   if (nzero > 0)
 587     {
 588       for (i = 0; i < nzero; i++)
 589         *(out++) = '0';
 590     }
 591
 592   /* Output digits after the decimal point, padding with zeros.  */
 593   if (nafter > 0)
 594     {
 595       if (nafter > ndigits)
 596         i = ndigits;
 597       else
 598         i = nafter;
 599
 600       memcpy (out, digits, i);
 601       while (i < nafter)
 602         out[i++] = '0';
 603
 604       digits += i;
 605       ndigits -= i;
 606       out += nafter;
 607     }
 608
 609   /* Output the exponent.  */
 610   if (expchar)
 611     {
 612       if (expchar != ' ')
 613         {
 614           *(out++) = expchar;
 615           edigits--;
 616         }
 617       snprintf (buffer, size, "%+0*d", edigits, e);
 618       memcpy (out, buffer, edigits);
 619     }
 620
 621   if (dtp->u.p.no_leading_blank)
 622     {
 623       out += edigits;
 624       memset( out , ' ' , nblanks );
 625       dtp->u.p.no_leading_blank = 0;
 626     }
 627
 628 #undef STR
 629 #undef STR1
 630 #undef MIN_FIELD_WIDTH
 631   return SUCCESS;
 632 }
 633
 634
 635 /* Write "Infinite" or "Nan" as appropriate for the given format.  */
 636
 637 static void
 638 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
 639 {
 640   char * p, fin;
 641   int nb = 0;
 642   sign_t sign;
 643   int mark;
 644
 645   if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
 646     {
 647       sign = calculate_sign (dtp, sign_bit);
 648       mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
 649
 650       nb =  f->u.real.w;
 651
 652       /* If the field width is zero, the processor must select a width
 653          not zero.  4 is chosen to allow output of '-Inf' or '+Inf' */
 654
 655       if ((nb == 0) || dtp->u.p.g0_no_blanks)
 656         {
 657           if (isnan_flag)
 658             nb = 3;
 659           else
 660             nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
 661         }
 662       p = write_block (dtp, nb);
 663       if (p == NULL)
 664         return;
 665       if (nb < 3)
 666         {
 667           if (unlikely (is_char4_unit (dtp)))
 668             {
 669               gfc_char4_t *p4 = (gfc_char4_t *) p;
 670               memset4 (p4, '*', nb);
 671             }
 672           else
 673             memset (p, '*', nb);
 674           return;
 675         }
 676
 677       if (unlikely (is_char4_unit (dtp)))
 678         {
 679           gfc_char4_t *p4 = (gfc_char4_t *) p;
 680           memset4 (p4, ' ', nb);
 681         }
 682       else
 683         memset(p, ' ', nb);
 684
 685       if (!isnan_flag)
 686         {
 687           if (sign_bit)
 688             {
 689               /* If the sign is negative and the width is 3, there is
 690                  insufficient room to output '-Inf', so output asterisks */
 691               if (nb == 3)
 692                 {
 693                   if (unlikely (is_char4_unit (dtp)))
 694                     {
 695                       gfc_char4_t *p4 = (gfc_char4_t *) p;
 696                       memset4 (p4, '*', nb);
 697                     }
 698                   else
 699                     memset (p, '*', nb);
 700                   return;
 701                 }
 702               /* The negative sign is mandatory */
 703               fin = '-';
 704             }
 705           else
 706             /* The positive sign is optional, but we output it for
 707                consistency */
 708             fin = '+';
 709
 710           if (unlikely (is_char4_unit (dtp)))
 711             {
 712               gfc_char4_t *p4 = (gfc_char4_t *) p;
 713
 714               if (nb > mark)
 715                 /* We have room, so output 'Infinity' */
 716                 memcpy4 (p4 + nb - 8, "Infinity", 8);
 717               else
 718                 /* For the case of width equals mark, there is not enough room
 719                    for the sign and 'Infinity' so we go with 'Inf' */
 720                 memcpy4 (p4 + nb - 3, "Inf", 3);
 721
 722               if (sign == S_PLUS || sign == S_MINUS)
 723                 {
 724                   if (nb < 9 && nb > 3)
 725                     /* Put the sign in front of Inf */
 726                     p4[nb - 4] = (gfc_char4_t) fin;
 727                   else if (nb > 8)
 728                     /* Put the sign in front of Infinity */
 729                     p4[nb - 9] = (gfc_char4_t) fin;
 730                 }
 731               return;
 732             }
 733
 734           if (nb > mark)
 735             /* We have room, so output 'Infinity' */
 736             memcpy(p + nb - 8, "Infinity", 8);
 737           else
 738             /* For the case of width equals 8, there is not enough room
 739                for the sign and 'Infinity' so we go with 'Inf' */
 740             memcpy(p + nb - 3, "Inf", 3);
 741
 742           if (sign == S_PLUS || sign == S_MINUS)
 743             {
 744               if (nb < 9 && nb > 3)
 745                 p[nb - 4] = fin;  /* Put the sign in front of Inf */
 746               else if (nb > 8)
 747                 p[nb - 9] = fin;  /* Put the sign in front of Infinity */
 748             }
 749         }
 750       else
 751         {
 752           if (unlikely (is_char4_unit (dtp)))
 753             {
 754               gfc_char4_t *p4 = (gfc_char4_t *) p;
 755               memcpy4 (p4 + nb - 3, "NaN", 3);
 756             }
 757           else
 758             memcpy(p + nb - 3, "NaN", 3);
 759         }
 760       return;
 761     }
 762 }
 763
 764
 765 /* Returns the value of 10**d.  */
 766
 767 #define CALCULATE_EXP(x) \
 768 inline static GFC_REAL_ ## x \
 769 calculate_exp_ ## x  (int d)\
 770 {\
 771   int i;\
 772   GFC_REAL_ ## x r = 1.0;\
 773   for (i = 0; i< (d >= 0 ? d : -d); i++)\
 774     r *= 10;\
 775   r = (d >= 0) ? r : 1.0 / r;\
 776   return r;\
 777 }
 778
 779 CALCULATE_EXP(4)
 780
 781 CALCULATE_EXP(8)
 782
 783 #ifdef HAVE_GFC_REAL_10
 784 CALCULATE_EXP(10)
 785 #endif
 786
 787 #ifdef HAVE_GFC_REAL_16
 788 CALCULATE_EXP(16)
 789 #endif
 790 #undef CALCULATE_EXP
 791
 792 /* Generate corresponding I/O format for FMT_G and output.
 793    The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
 794    LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
 795
 796    Data Magnitude                              Equivalent Conversion
 797    0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
 798    m = 0                                       F(w-n).(d-1), n' '
 799    0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
 800    1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
 801    10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
 802    ................                           ..........
 803    10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
 804    m >= 10**d-0.5                              Ew.d[Ee]
 805
 806    notes: for Gw.d ,  n' ' means 4 blanks
 807           for Gw.dEe, n' ' means e+2 blanks
 808           for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
 809           the asm volatile is required for 32-bit x86 platforms.  */
 810
 811 #define OUTPUT_FLOAT_FMT_G(x) \
 812 static void \
 813 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
 814                       GFC_REAL_ ## x m, char *buffer, size_t size, \
 815                       int sign_bit, bool zero_flag, int ndigits, int edigits) \
 816 { \
 817   int e = f->u.real.e;\
 818   int d = f->u.real.d;\
 819   int w = f->u.real.w;\
 820   fnode *newf;\
 821   GFC_REAL_ ## x rexp_d, r = 0.5;\
 822   int low, high, mid;\
 823   int ubound, lbound;\
 824   char *p, pad = ' ';\
 825   int save_scale_factor, nb = 0;\
 826   try result;\
 827 \
 828   save_scale_factor = dtp->u.p.scale_factor;\
 829   newf = (fnode *) get_mem (sizeof (fnode));\
 830 \
 831   switch (dtp->u.p.current_unit->round_status)\
 832     {\
 833       case ROUND_ZERO:\
 834         r = sign_bit ? 1.0 : 0.0;\
 835         break;\
 836       case ROUND_UP:\
 837         r = 1.0;\
 838         break;\
 839       case ROUND_DOWN:\
 840         r = 0.0;\
 841         break;\
 842       default:\
 843         break;\
 844     }\
 845 \
 846   rexp_d = calculate_exp_ ## x (-d);\
 847   if ((m > 0.0 && ((m < 0.1 - 0.1 * r * rexp_d) || (rexp_d * (m + r) >= 1.0)))\
 848       || ((m == 0.0) && !(compile_options.allow_std\
 849                           & (GFC_STD_F2003 | GFC_STD_F2008))))\
 850     { \
 851       newf->format = FMT_E;\
 852       newf->u.real.w = w;\
 853       newf->u.real.d = d;\
 854       newf->u.real.e = e;\
 855       nb = 0;\
 856       goto finish;\
 857     }\
 858 \
 859   mid = 0;\
 860   low = 0;\
 861   high = d + 1;\
 862   lbound = 0;\
 863   ubound = d + 1;\
 864 \
 865   while (low <= high)\
 866     { \
 867       GFC_REAL_ ## x temp;\
 868       mid = (low + high) / 2;\
 869 \
 870       temp = (calculate_exp_ ## x (mid - 1) * (1 - r * rexp_d));\
 871       asm volatile ("" : "+m" (temp));\
 872 \
 873       if (m < temp)\
 874         { \
 875           ubound = mid;\
 876           if (ubound == lbound + 1)\
 877             break;\
 878           high = mid - 1;\
 879         }\
 880       else if (m > temp)\
 881         { \
 882           lbound = mid;\
 883           if (ubound == lbound + 1)\
 884             { \
 885               mid ++;\
 886               break;\
 887             }\
 888           low = mid + 1;\
 889         }\
 890       else\
 891         {\
 892           mid++;\
 893           break;\
 894         }\
 895     }\
 896 \
 897   if (e > 4)\
 898     e = 4;\
 899   if (e < 0)\
 900     nb = 4;\
 901   else\
 902     nb = e + 2;\
 903 \
 904   nb = nb >= w ? 0 : nb;\
 905   newf->format = FMT_F;\
 906   newf->u.real.w = f->u.real.w - nb;\
 907 \
 908   if (m == 0.0)\
 909     newf->u.real.d = d - 1;\
 910   else\
 911     newf->u.real.d = - (mid - d - 1);\
 912 \
 913   dtp->u.p.scale_factor = 0;\
 914 \
 915  finish:\
 916   result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
 917                          ndigits, edigits);\
 918   dtp->u.p.scale_factor = save_scale_factor;\
 919 \
 920   free (newf);\
 921 \
 922   if (nb > 0 && !dtp->u.p.g0_no_blanks)\
 923     {\
 924       p = write_block (dtp, nb);\
 925       if (p == NULL)\
 926         return;\
 927       if (result == FAILURE)\
 928         pad = '*';\
 929       if (unlikely (is_char4_unit (dtp)))\
 930         {\
 931           gfc_char4_t *p4 = (gfc_char4_t *) p;\
 932           memset4 (p4, pad, nb);\
 933         }\
 934       else\
 935         memset (p, pad, nb);\
 936     }\
 937 }\
 938
 939 OUTPUT_FLOAT_FMT_G(4)
 940
 941 OUTPUT_FLOAT_FMT_G(8)
 942
 943 #ifdef HAVE_GFC_REAL_10
 944 OUTPUT_FLOAT_FMT_G(10)
 945 #endif
 946
 947 #ifdef HAVE_GFC_REAL_16
 948 OUTPUT_FLOAT_FMT_G(16)
 949 #endif
 950
 951 #undef OUTPUT_FLOAT_FMT_G
 952
 953
 954 /* Define a macro to build code for write_float.  */
 955
 956   /* Note: Before output_float is called, snprintf is used to print to buffer the
 957      number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
 958      (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
 959      before the decimal point.
 960
 961      #   The result will always contain a decimal point, even if no
 962          digits follow it
 963
 964      -   The converted value is to be left adjusted on the field boundary
 965
 966      +   A sign (+ or -) always be placed before a number
 967
 968      MIN_FIELD_WIDTH  minimum field width
 969
 970      *   (ndigits-1) is used as the precision
 971
 972      e format: [-]d.dddeÂ±dd where there is one digit before the
 973        decimal-point character and the number of digits after it is
 974        equal to the precision. The exponent always contains at least two
 975        digits; if the value is zero, the exponent is 00.  */
 976
 977 #define DTOA \
 978 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 979           "e", ndigits - 1, tmp);
 980
 981 #define DTOAL \
 982 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 983           "Le", ndigits - 1, tmp);
 984
 985
 986 #if defined(GFC_REAL_16_IS_FLOAT128)
 987 #define DTOAQ \
 988 __qmath_(quadmath_snprintf) (buffer, sizeof buffer, \
 989                              "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 990                              "Qe", ndigits - 1, tmp);
 991 #endif
 992
 993 #define WRITE_FLOAT(x,y)\
 994 {\
 995         GFC_REAL_ ## x tmp;\
 996         tmp = * (GFC_REAL_ ## x *)source;\
 997         sign_bit = signbit (tmp);\
 998         if (!isfinite (tmp))\
 999           { \
1000             write_infnan (dtp, f, isnan (tmp), sign_bit);\
1001             return;\
1002           }\
1003         tmp = sign_bit ? -tmp : tmp;\
1004         zero_flag = (tmp == 0.0);\
1005 \
1006         DTOA ## y\
1007 \
1008         if (f->format != FMT_G)\
1009           output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
1010                         edigits);\
1011         else \
1012           output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
1013                                     zero_flag, ndigits, edigits);\
1014 }\
1015
1016 /* Output a real number according to its format.  */
1017
1018 static void
1019 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
1020 {
1021
1022 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18
1023 # define MIN_FIELD_WIDTH 48
1024 #else
1025 # define MIN_FIELD_WIDTH 31
1026 #endif
1027 #define STR(x) STR1(x)
1028 #define STR1(x) #x
1029
1030   /* This must be large enough to accurately hold any value.  */
1031   char buffer[MIN_FIELD_WIDTH+1];
1032   int sign_bit, ndigits, edigits;
1033   bool zero_flag;
1034   size_t size;
1035
1036   size = MIN_FIELD_WIDTH+1;
1037
1038   /* printf pads blanks for us on the exponent so we just need it big enough
1039      to handle the largest number of exponent digits expected.  */
1040   edigits=4;
1041
1042   if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
1043       || ((f->format == FMT_D || f->format == FMT_E)
1044       && dtp->u.p.scale_factor != 0))
1045     {
1046       /* Always convert at full precision to avoid double rounding.  */
1047       ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1048     }
1049   else
1050     {
1051       /* The number of digits is known, so let printf do the rounding.  */
1052       if (f->format == FMT_ES)
1053         ndigits = f->u.real.d + 1;
1054       else
1055         ndigits = f->u.real.d;
1056       if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
1057         ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1058     }
1059
1060   switch (len)
1061     {
1062     case 4:
1063       WRITE_FLOAT(4,)
1064       break;
1065
1066     case 8:
1067       WRITE_FLOAT(8,)
1068       break;
1069
1070 #ifdef HAVE_GFC_REAL_10
1071     case 10:
1072       WRITE_FLOAT(10,L)
1073       break;
1074 #endif
1075 #ifdef HAVE_GFC_REAL_16
1076     case 16:
1077 # ifdef GFC_REAL_16_IS_FLOAT128
1078       WRITE_FLOAT(16,Q)
1079 # else
1080       WRITE_FLOAT(16,L)
1081 # endif
1082       break;
1083 #endif
1084     default:
1085       internal_error (NULL, "bad real kind");
1086     }
1087 }