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