libgfortran/io/write_float.def

   1 /* Copyright (C) 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   2    Contributed by Andy Vaught
   3    Write float code factoring to this file by Jerry DeLisle
   4    F2003 I/O support contributed by Jerry DeLisle
   5
   6 This file is part of the GNU Fortran runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU General Public License as published by
  10 the Free Software Foundation; either version 3, or (at your option)
  11 any later version.
  12
  13 Libgfortran is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27 #include "config.h"
  28
  29 typedef enum
  30 { S_NONE, S_MINUS, S_PLUS }
  31 sign_t;
  32
  33 /* Given a flag that indicates if a value is negative or not, return a
  34    sign_t that gives the sign that we need to produce.  */
  35
  36 static sign_t
  37 calculate_sign (st_parameter_dt *dtp, int negative_flag)
  38 {
  39   sign_t s = S_NONE;
  40
  41   if (negative_flag)
  42     s = S_MINUS;
  43   else
  44     switch (dtp->u.p.sign_status)
  45       {
  46       case SIGN_SP:     /* Show sign. */
  47         s = S_PLUS;
  48         break;
  49       case SIGN_SS:     /* Suppress sign. */
  50         s = S_NONE;
  51         break;
  52       case SIGN_S:      /* Processor defined. */
  53       case SIGN_UNSPECIFIED:
  54         s = options.optional_plus ? S_PLUS : S_NONE;
  55         break;
  56       }
  57
  58   return s;
  59 }
  60
  61
  62 /* Output a real number according to its format which is FMT_G free.  */
  63
  64 static try
  65 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
  66               int sign_bit, bool zero_flag, int ndigits, int edigits)
  67 {
  68   char *out;
  69   char *digits;
  70   int e;
  71   char expchar, rchar;
  72   format_token ft;
  73   int w;
  74   int d;
  75   /* Number of digits before the decimal point.  */
  76   int nbefore;
  77   /* Number of zeros after the decimal point.  */
  78   int nzero;
  79   /* Number of digits after the decimal point.  */
  80   int nafter;
  81   /* Number of zeros after the decimal point, whatever the precision.  */
  82   int nzero_real;
  83   int leadzero;
  84   int nblanks;
  85   int i;
  86   sign_t sign;
  87
  88   ft = f->format;
  89   w = f->u.real.w;
  90   d = f->u.real.d;
  91
  92   rchar = '5';
  93   nzero_real = -1;
  94
  95   /* We should always know the field width and precision.  */
  96   if (d < 0)
  97     internal_error (&dtp->common, "Unspecified precision");
  98
  99   sign = calculate_sign (dtp, sign_bit);
 100
 101   /* The following code checks the given string has punctuation in the correct
 102      places.  Uncomment if needed for debugging.
 103      if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
 104                     || buffer[ndigits + 2] != 'e'))
 105        internal_error (&dtp->common, "printf is broken");  */
 106
 107   /* Read the exponent back in.  */
 108   e = atoi (&buffer[ndigits + 3]) + 1;
 109
 110   /* Make sure zero comes out as 0.0e0.   */
 111   if (zero_flag)
 112     {
 113       e = 0;
 114       if (compile_options.sign_zero != 1)
 115         sign = calculate_sign (dtp, 0);
 116
 117       /* Handle special cases.  */
 118       if (w == 0)
 119         w = d + 1;
 120
 121       /* For this one we choose to not output a decimal point.
 122          F95 10.5.1.2.1  */
 123       if (w == 1 && ft == FMT_F)
 124         {
 125           out = write_block (dtp, w);
 126           if (out == NULL)
 127             return FAILURE;
 128
 129           if (unlikely (is_char4_unit (dtp)))
 130             {
 131               gfc_char4_t *out4 = (gfc_char4_t *) out;
 132               *out4 = '0';
 133               return SUCCESS;
 134             }
 135
 136           *out = '0';
 137           return SUCCESS;
 138         }
 139     }
 140
 141   /* Normalize the fractional component.  */
 142   buffer[2] = buffer[1];
 143   digits = &buffer[2];
 144
 145   /* Figure out where to place the decimal point.  */
 146   switch (ft)
 147     {
 148     case FMT_F:
 149       if (d == 0 && e <= 0 && dtp->u.p.scale_factor == 0)
 150         {
 151           memmove (digits + 1, digits, ndigits - 1);
 152           digits[0] = '0';
 153           e++;
 154         }
 155
 156       nbefore = e + dtp->u.p.scale_factor;
 157       if (nbefore < 0)
 158         {
 159           nzero = -nbefore;
 160           nzero_real = nzero;
 161           if (nzero > d)
 162             nzero = d;
 163           nafter = d - nzero;
 164           nbefore = 0;
 165         }
 166       else
 167         {
 168           nzero = 0;
 169           nafter = d;
 170         }
 171       expchar = 0;
 172       break;
 173
 174     case FMT_E:
 175     case FMT_D:
 176       i = dtp->u.p.scale_factor;
 177       if (d <= 0 && i == 0)
 178         {
 179           generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
 180                           "greater than zero in format specifier 'E' or 'D'");
 181           return FAILURE;
 182         }
 183       if (i <= -d || i >= d + 2)
 184         {
 185           generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
 186                           "out of range in format specifier 'E' or 'D'");
 187           return FAILURE;
 188         }
 189
 190       if (!zero_flag)
 191         e -= i;
 192       if (i < 0)
 193         {
 194           nbefore = 0;
 195           nzero = -i;
 196           nafter = d + i;
 197         }
 198       else if (i > 0)
 199         {
 200           nbefore = i;
 201           nzero = 0;
 202           nafter = (d - i) + 1;
 203         }
 204       else /* i == 0 */
 205         {
 206           nbefore = 0;
 207           nzero = 0;
 208           nafter = d;
 209         }
 210
 211       if (ft == FMT_E)
 212         expchar = 'E';
 213       else
 214         expchar = 'D';
 215       break;
 216
 217     case FMT_EN:
 218       /* The exponent must be a multiple of three, with 1-3 digits before
 219          the decimal point.  */
 220       if (!zero_flag)
 221         e--;
 222       if (e >= 0)
 223         nbefore = e % 3;
 224       else
 225         {
 226           nbefore = (-e) % 3;
 227           if (nbefore != 0)
 228             nbefore = 3 - nbefore;
 229         }
 230       e -= nbefore;
 231       nbefore++;
 232       nzero = 0;
 233       nafter = d;
 234       expchar = 'E';
 235       break;
 236
 237     case FMT_ES:
 238       if (!zero_flag)
 239         e--;
 240       nbefore = 1;
 241       nzero = 0;
 242       nafter = d;
 243       expchar = 'E';
 244       break;
 245
 246     default:
 247       /* Should never happen.  */
 248       internal_error (&dtp->common, "Unexpected format token");
 249     }
 250
 251   /* Round the value.  The value being rounded is an unsigned magnitude.
 252      The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
 253   switch (dtp->u.p.current_unit->round_status)
 254     {
 255       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
 256         goto skip;
 257       case ROUND_UP:
 258         if (sign_bit)
 259           goto skip;
 260         rchar = '0';
 261         break;
 262       case ROUND_DOWN:
 263         if (!sign_bit)
 264           goto skip;
 265         rchar = '0';
 266         break;
 267       case ROUND_NEAREST:
 268         /* Round compatible unless there is a tie. A tie is a 5 with
 269            all trailing zero's.  */
 270         i = nafter + nbefore;
 271         if (digits[i] == '5')
 272           {
 273             for(i++ ; i < ndigits; i++)
 274               {
 275                 if (digits[i] != '0')
 276                   goto do_rnd;
 277               }
 278             /* It is a  tie so round to even.  */
 279             switch (digits[nafter + nbefore - 1])
 280               {
 281                 case '1':
 282                 case '3':
 283                 case '5':
 284                 case '7':
 285                 case '9':
 286                   /* If odd, round away from zero to even.  */
 287                   break;
 288                 default:
 289                   /* If even, skip rounding, truncate to even.  */
 290                   goto skip;
 291               }
 292           }
 293          /* Fall through.  */
 294       case ROUND_PROCDEFINED:
 295       case ROUND_UNSPECIFIED:
 296       case ROUND_COMPATIBLE:
 297         rchar = '5';
 298         /* Just fall through and do the actual rounding.  */
 299     }
 300
 301   do_rnd:
 302
 303   if (nbefore + nafter == 0)
 304     {
 305       ndigits = 0;
 306       if (nzero_real == d && digits[0] >= rchar)
 307         {
 308           /* We rounded to zero but shouldn't have */
 309           nzero--;
 310           nafter = 1;
 311           digits[0] = '1';
 312           ndigits = 1;
 313         }
 314     }
 315   else if (nbefore + nafter < ndigits)
 316     {
 317       ndigits = nbefore + nafter;
 318       i = ndigits;
 319       if (digits[i] >= rchar)
 320         {
 321           /* Propagate the carry.  */
 322           for (i--; i >= 0; i--)
 323             {
 324               if (digits[i] != '9')
 325                 {
 326                   digits[i]++;
 327                   break;
 328                 }
 329               digits[i] = '0';
 330             }
 331
 332           if (i < 0)
 333             {
 334               /* The carry overflowed.  Fortunately we have some spare
 335                  space at the start of the buffer.  We may discard some
 336                  digits, but this is ok because we already know they are
 337                  zero.  */
 338               digits--;
 339               digits[0] = '1';
 340               if (ft == FMT_F)
 341                 {
 342                   if (nzero > 0)
 343                     {
 344                       nzero--;
 345                       nafter++;
 346                     }
 347                   else
 348                     nbefore++;
 349                 }
 350               else if (ft == FMT_EN)
 351                 {
 352                   nbefore++;
 353                   if (nbefore == 4)
 354                     {
 355                       nbefore = 1;
 356                       e += 3;
 357                     }
 358                 }
 359               else
 360                 e++;
 361             }
 362         }
 363     }
 364
 365   skip:
 366
 367   /* Calculate the format of the exponent field.  */
 368   if (expchar)
 369     {
 370       edigits = 1;
 371       for (i = abs (e); i >= 10; i /= 10)
 372         edigits++;
 373
 374       if (f->u.real.e < 0)
 375         {
 376           /* Width not specified.  Must be no more than 3 digits.  */
 377           if (e > 999 || e < -999)
 378             edigits = -1;
 379           else
 380             {
 381               edigits = 4;
 382               if (e > 99 || e < -99)
 383                 expchar = ' ';
 384             }
 385         }
 386       else
 387         {
 388           /* Exponent width specified, check it is wide enough.  */
 389           if (edigits > f->u.real.e)
 390             edigits = -1;
 391           else
 392             edigits = f->u.real.e + 2;
 393         }
 394     }
 395   else
 396     edigits = 0;
 397
 398   /* Zero values always output as positive, even if the value was negative
 399      before rounding.  */
 400   for (i = 0; i < ndigits; i++)
 401     {
 402       if (digits[i] != '0')
 403         break;
 404     }
 405   if (i == ndigits)
 406     {
 407       /* The output is zero, so set the sign according to the sign bit unless
 408          -fno-sign-zero was specified.  */
 409       if (compile_options.sign_zero == 1)
 410         sign = calculate_sign (dtp, sign_bit);
 411       else
 412         sign = calculate_sign (dtp, 0);
 413     }
 414
 415   /* Pick a field size if none was specified.  */
 416   if (w <= 0)
 417     {
 418       w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
 419       w = w == 1 ? 2 : w;
 420     }
 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 || w == 1 || (w == 2 && sign != S_NONE))
 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   sign_t sign;
 664   int mark;
 665
 666   if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
 667     {
 668       sign = calculate_sign (dtp, sign_bit);
 669       mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
 670
 671       nb =  f->u.real.w;
 672
 673       /* If the field width is zero, the processor must select a width
 674          not zero.  4 is chosen to allow output of '-Inf' or '+Inf' */
 675
 676       if (nb == 0)
 677         {
 678           if (isnan_flag)
 679             nb = 3;
 680           else
 681             nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
 682         }
 683       p = write_block (dtp, nb);
 684       if (p == NULL)
 685         return;
 686       if (nb < 3)
 687         {
 688           if (unlikely (is_char4_unit (dtp)))
 689             {
 690               gfc_char4_t *p4 = (gfc_char4_t *) p;
 691               memset4 (p4, '*', nb);
 692             }
 693           else
 694             memset (p, '*', nb);
 695           return;
 696         }
 697
 698       if (unlikely (is_char4_unit (dtp)))
 699         {
 700           gfc_char4_t *p4 = (gfc_char4_t *) p;
 701           memset4 (p4, ' ', nb);
 702         }
 703       else
 704         memset(p, ' ', nb);
 705
 706       if (!isnan_flag)
 707         {
 708           if (sign_bit)
 709             {
 710               /* If the sign is negative and the width is 3, there is
 711                  insufficient room to output '-Inf', so output asterisks */
 712               if (nb == 3)
 713                 {
 714                   if (unlikely (is_char4_unit (dtp)))
 715                     {
 716                       gfc_char4_t *p4 = (gfc_char4_t *) p;
 717                       memset4 (p4, '*', nb);
 718                     }
 719                   else
 720                     memset (p, '*', nb);
 721                   return;
 722                 }
 723               /* The negative sign is mandatory */
 724               fin = '-';
 725             }
 726           else
 727             /* The positive sign is optional, but we output it for
 728                consistency */
 729             fin = '+';
 730
 731           if (unlikely (is_char4_unit (dtp)))
 732             {
 733               gfc_char4_t *p4 = (gfc_char4_t *) p;
 734
 735               if (nb > mark)
 736                 /* We have room, so output 'Infinity' */
 737                 memcpy4 (p4 + nb - 8, "Infinity", 8);
 738               else
 739                 /* For the case of width equals mark, there is not enough room
 740                    for the sign and 'Infinity' so we go with 'Inf' */
 741                 memcpy4 (p4 + nb - 3, "Inf", 3);
 742
 743               if (sign == S_PLUS || sign == S_MINUS)
 744                 {
 745                   if (nb < 9 && nb > 3)
 746                     /* Put the sign in front of Inf */
 747                     p4[nb - 4] = (gfc_char4_t) fin;
 748                   else if (nb > 8)
 749                     /* Put the sign in front of Infinity */
 750                     p4[nb - 9] = (gfc_char4_t) fin;
 751                 }
 752               return;
 753             }
 754
 755           if (nb > mark)
 756             /* We have room, so output 'Infinity' */
 757             memcpy(p + nb - 8, "Infinity", 8);
 758           else
 759             /* For the case of width equals 8, there is not enough room
 760                for the sign and 'Infinity' so we go with 'Inf' */
 761             memcpy(p + nb - 3, "Inf", 3);
 762
 763           if (sign == S_PLUS || sign == S_MINUS)
 764             {
 765               if (nb < 9 && nb > 3)
 766                 p[nb - 4] = fin;  /* Put the sign in front of Inf */
 767               else if (nb > 8)
 768                 p[nb - 9] = fin;  /* Put the sign in front of Infinity */
 769             }
 770         }
 771       else
 772         {
 773           if (unlikely (is_char4_unit (dtp)))
 774             {
 775               gfc_char4_t *p4 = (gfc_char4_t *) p;
 776               memcpy4 (p4 + nb - 3, "NaN", 3);
 777             }
 778           else
 779             memcpy(p + nb - 3, "NaN", 3);
 780         }
 781       return;
 782     }
 783 }
 784
 785
 786 /* Returns the value of 10**d.  */
 787
 788 #define CALCULATE_EXP(x) \
 789 inline static GFC_REAL_ ## x \
 790 calculate_exp_ ## x  (int d)\
 791 {\
 792   int i;\
 793   GFC_REAL_ ## x r = 1.0;\
 794   for (i = 0; i< (d >= 0 ? d : -d); i++)\
 795     r *= 10;\
 796   r = (d >= 0) ? r : 1.0 / r;\
 797   return r;\
 798 }
 799
 800 CALCULATE_EXP(4)
 801
 802 CALCULATE_EXP(8)
 803
 804 #ifdef HAVE_GFC_REAL_10
 805 CALCULATE_EXP(10)
 806 #endif
 807
 808 #ifdef HAVE_GFC_REAL_16
 809 CALCULATE_EXP(16)
 810 #endif
 811 #undef CALCULATE_EXP
 812
 813 /* Generate corresponding I/O format for FMT_G and output.
 814    The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
 815    LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
 816
 817    Data Magnitude                              Equivalent Conversion
 818    0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
 819    m = 0                                       F(w-n).(d-1), n' '
 820    0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
 821    1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
 822    10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
 823    ................                           ..........
 824    10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
 825    m >= 10**d-0.5                              Ew.d[Ee]
 826
 827    notes: for Gw.d ,  n' ' means 4 blanks
 828           for Gw.dEe, n' ' means e+2 blanks  */
 829
 830 #define OUTPUT_FLOAT_FMT_G(x) \
 831 static void \
 832 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
 833                       GFC_REAL_ ## x m, char *buffer, size_t size, \
 834                       int sign_bit, bool zero_flag, int ndigits, int edigits) \
 835 { \
 836   int e = f->u.real.e;\
 837   int d = f->u.real.d;\
 838   int w = f->u.real.w;\
 839   fnode *newf;\
 840   GFC_REAL_ ## x rexp_d;\
 841   int low, high, mid;\
 842   int ubound, lbound;\
 843   char *p, pad = ' ';\
 844   int save_scale_factor, nb = 0;\
 845   try result;\
 846 \
 847   save_scale_factor = dtp->u.p.scale_factor;\
 848   newf = (fnode *) get_mem (sizeof (fnode));\
 849 \
 850   rexp_d = calculate_exp_ ## x (-d);\
 851   if ((m > 0.0 && m < 0.1 - 0.05 * rexp_d) || (rexp_d * (m + 0.5) >= 1.0) ||\
 852       ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
 853     { \
 854       newf->format = FMT_E;\
 855       newf->u.real.w = w;\
 856       newf->u.real.d = d;\
 857       newf->u.real.e = e;\
 858       nb = 0;\
 859       goto finish;\
 860     }\
 861 \
 862   mid = 0;\
 863   low = 0;\
 864   high = d + 1;\
 865   lbound = 0;\
 866   ubound = d + 1;\
 867 \
 868   while (low <= high)\
 869     { \
 870       GFC_REAL_ ## x temp;\
 871       mid = (low + high) / 2;\
 872 \
 873       temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
 874 \
 875       if (m < temp)\
 876         { \
 877           ubound = mid;\
 878           if (ubound == lbound + 1)\
 879             break;\
 880           high = mid - 1;\
 881         }\
 882       else if (m > temp)\
 883         { \
 884           lbound = mid;\
 885           if (ubound == lbound + 1)\
 886             { \
 887               mid ++;\
 888               break;\
 889             }\
 890           low = mid + 1;\
 891         }\
 892       else\
 893         {\
 894           mid++;\
 895           break;\
 896         }\
 897     }\
 898 \
 899   if (e > 4)\
 900     e = 4;\
 901   if (e < 0)\
 902     nb = 4;\
 903   else\
 904     nb = e + 2;\
 905 \
 906   nb = nb >= w ? 0 : nb;\
 907   newf->format = FMT_F;\
 908   newf->u.real.w = f->u.real.w - nb;\
 909 \
 910   if (m == 0.0)\
 911     newf->u.real.d = d - 1;\
 912   else\
 913     newf->u.real.d = - (mid - d - 1);\
 914 \
 915   dtp->u.p.scale_factor = 0;\
 916 \
 917  finish:\
 918   result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
 919                          ndigits, edigits);\
 920   dtp->u.p.scale_factor = save_scale_factor;\
 921 \
 922   free (newf);\
 923 \
 924   if (nb > 0 && !dtp->u.p.g0_no_blanks)\
 925     {\
 926       p = write_block (dtp, nb);\
 927       if (p == NULL)\
 928         return;\
 929       if (result == FAILURE)\
 930         pad = '*';\
 931       if (unlikely (is_char4_unit (dtp)))\
 932         {\
 933           gfc_char4_t *p4 = (gfc_char4_t *) p;\
 934           memset4 (p4, pad, nb);\
 935         }\
 936       else\
 937         memset (p, pad, nb);\
 938     }\
 939 }\
 940
 941 OUTPUT_FLOAT_FMT_G(4)
 942
 943 OUTPUT_FLOAT_FMT_G(8)
 944
 945 #ifdef HAVE_GFC_REAL_10
 946 OUTPUT_FLOAT_FMT_G(10)
 947 #endif
 948
 949 #ifdef HAVE_GFC_REAL_16
 950 OUTPUT_FLOAT_FMT_G(16)
 951 #endif
 952
 953 #undef OUTPUT_FLOAT_FMT_G
 954
 955
 956 /* Define a macro to build code for write_float.  */
 957
 958   /* Note: Before output_float is called, sprintf is used to print to buffer the
 959      number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
 960      (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
 961      before the decimal point.
 962
 963      #   The result will always contain a decimal point, even if no
 964          digits follow it
 965
 966      -   The converted value is to be left adjusted on the field boundary
 967
 968      +   A sign (+ or -) always be placed before a number
 969
 970      MIN_FIELD_WIDTH  minimum field width
 971
 972      *   (ndigits-1) is used as the precision
 973
 974      e format: [-]d.dddeÂ±dd where there is one digit before the
 975        decimal-point character and the number of digits after it is
 976        equal to the precision. The exponent always contains at least two
 977        digits; if the value is zero, the exponent is 00.  */
 978
 979 #ifdef HAVE_SNPRINTF
 980
 981 #define DTOA \
 982 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 983           "e", ndigits - 1, tmp);
 984
 985 #define DTOAL \
 986 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 987           "Le", ndigits - 1, tmp);
 988
 989 #else
 990
 991 #define DTOA \
 992 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 993          "e", ndigits - 1, tmp);
 994
 995 #define DTOAL \
 996 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
 997          "Le", ndigits - 1, tmp);
 998
 999 #endif
1000
1001 #if defined(GFC_REAL_16_IS_FLOAT128)
1002 #define DTOAQ \
1003 __qmath_(quadmath_snprintf) (buffer, sizeof buffer, \
1004                              "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
1005                              "Qe", ndigits - 1, tmp);
1006 #endif
1007
1008 #define WRITE_FLOAT(x,y)\
1009 {\
1010         GFC_REAL_ ## x tmp;\
1011         tmp = * (GFC_REAL_ ## x *)source;\
1012         sign_bit = signbit (tmp);\
1013         if (!isfinite (tmp))\
1014           { \
1015             write_infnan (dtp, f, isnan (tmp), sign_bit);\
1016             return;\
1017           }\
1018         tmp = sign_bit ? -tmp : tmp;\
1019         zero_flag = (tmp == 0.0);\
1020 \
1021         DTOA ## y\
1022 \
1023         if (f->format != FMT_G)\
1024           output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
1025                         edigits);\
1026         else \
1027           output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
1028                                     zero_flag, ndigits, edigits);\
1029 }\
1030
1031 /* Output a real number according to its format.  */
1032
1033 static void
1034 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
1035 {
1036
1037 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18
1038 # define MIN_FIELD_WIDTH 46
1039 #else
1040 # define MIN_FIELD_WIDTH 31
1041 #endif
1042 #define STR(x) STR1(x)
1043 #define STR1(x) #x
1044
1045   /* This must be large enough to accurately hold any value.  */
1046   char buffer[MIN_FIELD_WIDTH+1];
1047   int sign_bit, ndigits, edigits;
1048   bool zero_flag;
1049   size_t size;
1050
1051   size = MIN_FIELD_WIDTH+1;
1052
1053   /* printf pads blanks for us on the exponent so we just need it big enough
1054      to handle the largest number of exponent digits expected.  */
1055   edigits=4;
1056
1057   if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
1058       || ((f->format == FMT_D || f->format == FMT_E)
1059       && dtp->u.p.scale_factor != 0))
1060     {
1061       /* Always convert at full precision to avoid double rounding.  */
1062       ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1063     }
1064   else
1065     {
1066       /* The number of digits is known, so let printf do the rounding.  */
1067       if (f->format == FMT_ES)
1068         ndigits = f->u.real.d + 1;
1069       else
1070         ndigits = f->u.real.d;
1071       if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
1072         ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1073     }
1074
1075   switch (len)
1076     {
1077     case 4:
1078       WRITE_FLOAT(4,)
1079       break;
1080
1081     case 8:
1082       WRITE_FLOAT(8,)
1083       break;
1084
1085 #ifdef HAVE_GFC_REAL_10
1086     case 10:
1087       WRITE_FLOAT(10,L)
1088       break;
1089 #endif
1090 #ifdef HAVE_GFC_REAL_16
1091     case 16:
1092 # ifdef GFC_REAL_16_IS_FLOAT128
1093       WRITE_FLOAT(16,Q)
1094 # else
1095       WRITE_FLOAT(16,L)
1096 # endif
1097       break;
1098 #endif
1099     default:
1100       internal_error (NULL, "bad real kind");
1101     }
1102 }