libgfortran/io/write_float.def

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