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 void
  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;
 130
 131           if (unlikely (is_char4_unit (dtp)))
 132             {
 133               gfc_char4_t *out4 = (gfc_char4_t *) out;
 134               *out4 = '0';
 135               return;
 136             }
 137
 138           *out = '0';
 139           return;
 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;
 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;
 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;
 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;
 446         }
 447       star_fill (out, w);
 448       return;
 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;
 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 }
 653
 654
 655 /* Write "Infinite" or "Nan" as appropriate for the given format.  */
 656
 657 static void
 658 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
 659 {
 660   char * p, fin;
 661   int nb = 0;
 662
 663   if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
 664     {
 665       nb =  f->u.real.w;
 666
 667       /* If the field width is zero, the processor must select a width
 668          not zero.  4 is chosen to allow output of '-Inf' or '+Inf' */
 669
 670       if (nb == 0) nb = 4;
 671       p = write_block (dtp, nb);
 672       if (p == NULL)
 673         return;
 674       if (nb < 3)
 675         {
 676           if (unlikely (is_char4_unit (dtp)))
 677             {
 678               gfc_char4_t *p4 = (gfc_char4_t *) p;
 679               memset4 (p4, '*', nb);
 680             }
 681           else
 682             memset (p, '*', nb);
 683           return;
 684         }
 685
 686       if (unlikely (is_char4_unit (dtp)))
 687         {
 688           gfc_char4_t *p4 = (gfc_char4_t *) p;
 689           memset4 (p4, ' ', nb);
 690         }
 691       else
 692         memset(p, ' ', nb);
 693
 694       if (!isnan_flag)
 695         {
 696           if (sign_bit)
 697             {
 698               /* If the sign is negative and the width is 3, there is
 699                  insufficient room to output '-Inf', so output asterisks */
 700               if (nb == 3)
 701                 {
 702                   if (unlikely (is_char4_unit (dtp)))
 703                     {
 704                       gfc_char4_t *p4 = (gfc_char4_t *) p;
 705                       memset4 (p4, '*', nb);
 706                     }
 707                   else
 708                     memset (p, '*', nb);
 709                   return;
 710                 }
 711               /* The negative sign is mandatory */
 712               fin = '-';
 713             }
 714           else
 715             /* The positive sign is optional, but we output it for
 716                consistency */
 717             fin = '+';
 718
 719           if (unlikely (is_char4_unit (dtp)))
 720             {
 721               gfc_char4_t *p4 = (gfc_char4_t *) p;
 722               if (nb > 8)
 723                 /* We have room, so output 'Infinity' */
 724                 memcpy4 (p4 + nb - 8, "Infinity", 8);
 725               else
 726                 /* For the case of width equals 8, there is not enough room
 727                    for the sign and 'Infinity' so we go with 'Inf' */
 728                 memcpy4 (p4 + nb - 3, "Inf", 3);
 729
 730               if (nb < 9 && nb > 3)
 731                 /* Put the sign in front of Inf */
 732                 p4[nb - 4] = (gfc_char4_t) fin;
 733               else if (nb > 8)
 734                 /* Put the sign in front of Infinity */
 735                 p4[nb - 9] = (gfc_char4_t) fin;
 736               return;
 737             }
 738
 739           if (nb > 8)
 740             /* We have room, so output 'Infinity' */
 741             memcpy(p + nb - 8, "Infinity", 8);
 742           else
 743             /* For the case of width equals 8, there is not enough room
 744                for the sign and 'Infinity' so we go with 'Inf' */
 745             memcpy(p + nb - 3, "Inf", 3);
 746
 747           if (nb < 9 && nb > 3)
 748             p[nb - 4] = fin;  /* Put the sign in front of Inf */
 749           else if (nb > 8)
 750             p[nb - 9] = fin;  /* Put the sign in front of Infinity */
 751         }
 752       else
 753         {
 754           if (unlikely (is_char4_unit (dtp)))
 755             {
 756               gfc_char4_t *p4 = (gfc_char4_t *) p;
 757               memcpy4 (p4 + nb - 3, "NaN", 3);
 758             }
 759           else
 760             memcpy(p + nb - 3, "NaN", 3);
 761         }
 762       return;
 763     }
 764 }
 765
 766
 767 /* Returns the value of 10**d.  */
 768
 769 #define CALCULATE_EXP(x) \
 770 inline static GFC_REAL_ ## x \
 771 calculate_exp_ ## x  (int d)\
 772 {\
 773   int i;\
 774   GFC_REAL_ ## x r = 1.0;\
 775   for (i = 0; i< (d >= 0 ? d : -d); i++)\
 776     r *= 10;\
 777   r = (d >= 0) ? r : 1.0 / r;\
 778   return r;\
 779 }
 780
 781 CALCULATE_EXP(4)
 782
 783 CALCULATE_EXP(8)
 784
 785 #ifdef HAVE_GFC_REAL_10
 786 CALCULATE_EXP(10)
 787 #endif
 788
 789 #ifdef HAVE_GFC_REAL_16
 790 CALCULATE_EXP(16)
 791 #endif
 792 #undef CALCULATE_EXP
 793
 794 /* Generate corresponding I/O format for FMT_G and output.
 795    The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
 796    LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
 797
 798    Data Magnitude                              Equivalent Conversion
 799    0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
 800    m = 0                                       F(w-n).(d-1), n' '
 801    0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
 802    1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
 803    10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
 804    ................                           ..........
 805    10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
 806    m >= 10**d-0.5                              Ew.d[Ee]
 807
 808    notes: for Gw.d ,  n' ' means 4 blanks
 809           for Gw.dEe, n' ' means e+2 blanks  */
 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;\
 822   int low, high, mid;\
 823   int ubound, lbound;\
 824   char *p;\
 825   int save_scale_factor, nb = 0;\
 826 \
 827   save_scale_factor = dtp->u.p.scale_factor;\
 828   newf = (fnode *) get_mem (sizeof (fnode));\
 829 \
 830   rexp_d = calculate_exp_ ## x (-d);\
 831   if ((m > 0.0 && m < 0.1 - 0.05 * rexp_d) || (rexp_d * (m + 0.5) >= 1.0) ||\
 832       ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
 833     { \
 834       newf->format = FMT_E;\
 835       newf->u.real.w = w;\
 836       newf->u.real.d = d;\
 837       newf->u.real.e = e;\
 838       nb = 0;\
 839       goto finish;\
 840     }\
 841 \
 842   mid = 0;\
 843   low = 0;\
 844   high = d + 1;\
 845   lbound = 0;\
 846   ubound = d + 1;\
 847 \
 848   while (low <= high)\
 849     { \
 850       GFC_REAL_ ## x temp;\
 851       mid = (low + high) / 2;\
 852 \
 853       temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
 854 \
 855       if (m < temp)\
 856         { \
 857           ubound = mid;\
 858           if (ubound == lbound + 1)\
 859             break;\
 860           high = mid - 1;\
 861         }\
 862       else if (m > temp)\
 863         { \
 864           lbound = mid;\
 865           if (ubound == lbound + 1)\
 866             { \
 867               mid ++;\
 868               break;\
 869             }\
 870           low = mid + 1;\
 871         }\
 872       else\
 873         {\
 874           mid++;\
 875           break;\
 876         }\
 877     }\
 878 \
 879   if (e < 0)\
 880     nb = 4;\
 881   else\
 882     nb = e + 2;\
 883 \
 884   newf->format = FMT_F;\
 885   newf->u.real.w = f->u.real.w - nb;\
 886 \
 887   if (m == 0.0)\
 888     newf->u.real.d = d - 1;\
 889   else\
 890     newf->u.real.d = - (mid - d - 1);\
 891 \
 892   dtp->u.p.scale_factor = 0;\
 893 \
 894  finish:\
 895   output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
 896                 edigits);\
 897   dtp->u.p.scale_factor = save_scale_factor;\
 898 \
 899   free (newf);\
 900 \
 901   if (nb > 0 && !dtp->u.p.g0_no_blanks)\
 902     {\
 903       p = write_block (dtp, nb);\
 904       if (p == NULL)\
 905         return;\
 906       if (unlikely (is_char4_unit (dtp)))\
 907         {\
 908           gfc_char4_t *p4 = (gfc_char4_t *) p;\
 909           memset4 (p4, ' ', nb);\
 910         }\
 911       else\
 912         memset (p, ' ', nb);\
 913     }\
 914 }\
 915
 916 OUTPUT_FLOAT_FMT_G(4)
 917
 918 OUTPUT_FLOAT_FMT_G(8)
 919
 920 #ifdef HAVE_GFC_REAL_10
 921 OUTPUT_FLOAT_FMT_G(10)
 922 #endif
 923
 924 #ifdef HAVE_GFC_REAL_16
 925 OUTPUT_FLOAT_FMT_G(16)
 926 #endif
 927
 928 #undef OUTPUT_FLOAT_FMT_G
 929
 930
 931 /* Define a macro to build code for write_float.  */
 932
 933   /* Note: Before output_float is called, sprintf is used to print to buffer the
 934      number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
 935      (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
 936      before the decimal point.
 937
 938      #   The result will always contain a decimal point, even if no
 939          digits follow it
 940
 941      -   The converted value is to be left adjusted on the field boundary
 942
 943      +   A sign (+ or -) always be placed before a number
 944
 945      MIN_FIELD_WIDTH  minimum field width
 946
 947      *   (ndigits-1) is used as the precision
 948
 949      e format: [-]d.dddeÂ±dd where there is one digit before the
 950        decimal-point character and the number of digits after it is
 951        equal to the precision. The exponent always contains at least two
 952        digits; if the value is zero, the exponent is 00.  */
 953
 954 #ifdef HAVE_SNPRINTF
 955
 956 #define DTOA \
 957 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 958           "e", ndigits - 1, tmp);
 959
 960 #define DTOAL \
 961 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 962           "Le", ndigits - 1, tmp);
 963
 964 #else
 965
 966 #define DTOA \
 967 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 968          "e", ndigits - 1, tmp);
 969
 970 #define DTOAL \
 971 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 972          "Le", ndigits - 1, tmp);
 973
 974 #endif
 975
 976 #define WRITE_FLOAT(x,y)\
 977 {\
 978         GFC_REAL_ ## x tmp;\
 979         tmp = * (GFC_REAL_ ## x *)source;\
 980         sign_bit = __builtin_signbit (tmp);\
 981         if (!isfinite (tmp))\
 982           { \
 983             write_infnan (dtp, f, isnan (tmp), sign_bit);\
 984             return;\
 985           }\
 986         tmp = sign_bit ? -tmp : tmp;\
 987         zero_flag = (tmp == 0.0);\
 988 \
 989         DTOA ## y\
 990 \
 991         if (f->format != FMT_G)\
 992           output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
 993                         edigits);\
 994         else \
 995           output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
 996                                     zero_flag, ndigits, edigits);\
 997 }\
 998
 999 /* Output a real number according to its format.  */
1000
1001 static void
1002 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
1003 {
1004
1005 #if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
1006 # define MIN_FIELD_WIDTH 46
1007 #else
1008 # define MIN_FIELD_WIDTH 31
1009 #endif
1010 #define STR(x) STR1(x)
1011 #define STR1(x) #x
1012
1013   /* This must be large enough to accurately hold any value.  */
1014   char buffer[MIN_FIELD_WIDTH+1];
1015   int sign_bit, ndigits, edigits;
1016   bool zero_flag;
1017   size_t size;
1018
1019   size = MIN_FIELD_WIDTH+1;
1020
1021   /* printf pads blanks for us on the exponent so we just need it big enough
1022      to handle the largest number of exponent digits expected.  */
1023   edigits=4;
1024
1025   if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
1026       || ((f->format == FMT_D || f->format == FMT_E)
1027       && dtp->u.p.scale_factor != 0))
1028     {
1029       /* Always convert at full precision to avoid double rounding.  */
1030       ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1031     }
1032   else
1033     {
1034       /* The number of digits is known, so let printf do the rounding.  */
1035       if (f->format == FMT_ES)
1036         ndigits = f->u.real.d + 1;
1037       else
1038         ndigits = f->u.real.d;
1039       if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
1040         ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1041     }
1042
1043   switch (len)
1044     {
1045     case 4:
1046       WRITE_FLOAT(4,)
1047       break;
1048
1049     case 8:
1050       WRITE_FLOAT(8,)
1051       break;
1052
1053 #ifdef HAVE_GFC_REAL_10
1054     case 10:
1055       WRITE_FLOAT(10,L)
1056       break;
1057 #endif
1058 #ifdef HAVE_GFC_REAL_16
1059     case 16:
1060       WRITE_FLOAT(16,L)
1061       break;
1062 #endif
1063     default:
1064       internal_error (NULL, "bad real kind");
1065     }
1066 }