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
6 This file is part of the GNU Fortran runtime library (libgfortran).
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)
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.
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.
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/>. */
30 { S_NONE, S_MINUS, S_PLUS }
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. */
37 calculate_sign (st_parameter_dt *dtp, int negative_flag)
44 switch (dtp->u.p.sign_status)
46 case SIGN_SP: /* Show sign. */
49 case SIGN_SS: /* Suppress sign. */
52 case SIGN_S: /* Processor defined. */
53 case SIGN_UNSPECIFIED:
54 s = options.optional_plus ? S_PLUS : S_NONE;
62 /* Output a real number according to its format which is FMT_G free. */
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)
75 /* Number of digits before the decimal point. */
77 /* Number of zeros after the decimal point. */
79 /* Number of digits after the decimal point. */
81 /* Number of zeros after the decimal point, whatever the precision. */
95 /* We should always know the field width and precision. */
97 internal_error (&dtp->common, "Unspecified precision");
99 sign = calculate_sign (dtp, sign_bit);
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"); */
107 /* Read the exponent back in. */
108 e = atoi (&buffer[ndigits + 3]) + 1;
110 /* Make sure zero comes out as 0.0e0. */
114 if (compile_options.sign_zero != 1)
115 sign = calculate_sign (dtp, 0);
117 /* Handle special cases. */
121 /* For this one we choose to not output a decimal point.
123 if (w == 1 && ft == FMT_F)
125 out = write_block (dtp, w);
129 if (unlikely (is_char4_unit (dtp)))
131 gfc_char4_t *out4 = (gfc_char4_t *) out;
141 /* Normalize the fractional component. */
142 buffer[2] = buffer[1];
145 /* Figure out where to place the decimal point. */
149 if (d == 0 && e <= 0 && dtp->u.p.scale_factor == 0)
151 memmove (digits + 1, digits, ndigits - 1);
156 nbefore = e + dtp->u.p.scale_factor;
176 i = dtp->u.p.scale_factor;
177 if (d <= 0 && i == 0)
179 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
180 "greater than zero in format specifier 'E' or 'D'");
183 if (i <= -d || i >= d + 2)
185 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
186 "out of range in format specifier 'E' or 'D'");
202 nafter = (d - i) + 1;
218 /* The exponent must be a multiple of three, with 1-3 digits before
219 the decimal point. */
228 nbefore = 3 - nbefore;
247 /* Should never happen. */
248 internal_error (&dtp->common, "Unexpected format token");
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)
255 case ROUND_ZERO: /* Do nothing and truncation occurs. */
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')
273 for(i++ ; i < ndigits; i++)
275 if (digits[i] != '0')
278 /* It is a tie so round to even. */
279 switch (digits[nafter + nbefore - 1])
286 /* If odd, round away from zero to even. */
289 /* If even, skip rounding, truncate to even. */
294 case ROUND_PROCDEFINED:
295 case ROUND_UNSPECIFIED:
296 case ROUND_COMPATIBLE:
298 /* Just fall through and do the actual rounding. */
303 if (nbefore + nafter == 0)
306 if (nzero_real == d && digits[0] >= rchar)
308 /* We rounded to zero but shouldn't have */
315 else if (nbefore + nafter < ndigits)
317 ndigits = nbefore + nafter;
319 if (digits[i] >= rchar)
321 /* Propagate the carry. */
322 for (i--; i >= 0; i--)
324 if (digits[i] != '9')
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
350 else if (ft == FMT_EN)
367 /* Calculate the format of the exponent field. */
371 for (i = abs (e); i >= 10; i /= 10)
376 /* Width not specified. Must be no more than 3 digits. */
377 if (e > 999 || e < -999)
382 if (e > 99 || e < -99)
388 /* Exponent width specified, check it is wide enough. */
389 if (edigits > f->u.real.e)
392 edigits = f->u.real.e + 2;
398 /* Zero values always output as positive, even if the value was negative
400 for (i = 0; i < ndigits; i++)
402 if (digits[i] != '0')
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);
412 sign = calculate_sign (dtp, 0);
415 /* Pick a field size if none was specified. */
418 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
422 /* Work out how much padding is needed. */
423 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
427 if (dtp->u.p.g0_no_blanks)
433 /* Create the ouput buffer. */
434 out = write_block (dtp, w);
438 /* Check the value fits in the specified field width. */
439 if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
441 if (unlikely (is_char4_unit (dtp)))
443 gfc_char4_t *out4 = (gfc_char4_t *) out;
444 memset4 (out4, '*', w);
451 /* See if we have space for a zero before the decimal point. */
452 if (nbefore == 0 && nblanks > 0)
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)))
465 gfc_char4_t *out4 = (gfc_char4_t *) out;
466 /* Pad to full field width. */
468 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
470 memset4 (out4, ' ', nblanks);
474 /* Output the initial sign (if any). */
477 else if (sign == S_MINUS)
480 /* Output an optional leading zero. */
484 /* Output the part before the decimal point, padding with zeros. */
487 if (nbefore > ndigits)
490 memcpy4 (out4, digits, i);
498 memcpy4 (out4, digits, i);
506 /* Output the decimal point. */
507 *(out4++) = dtp->u.p.current_unit->decimal_status
508 == DECIMAL_POINT ? '.' : ',';
510 /* Output leading zeros after the decimal point. */
513 for (i = 0; i < nzero; i++)
517 /* Output digits after the decimal point, padding with zeros. */
520 if (nafter > ndigits)
525 memcpy4 (out4, digits, i);
534 /* Output the exponent. */
543 snprintf (buffer, size, "%+0*d", edigits, e);
545 sprintf (buffer, "%+0*d", edigits, e);
547 memcpy4 (out4, buffer, edigits);
550 if (dtp->u.p.no_leading_blank)
553 memset4 (out4, ' ' , nblanks);
554 dtp->u.p.no_leading_blank = 0;
557 } /* End of character(kind=4) internal unit code. */
559 /* Pad to full field width. */
561 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
563 memset (out, ' ', nblanks);
567 /* Output the initial sign (if any). */
570 else if (sign == S_MINUS)
573 /* Output an optional leading zero. */
577 /* Output the part before the decimal point, padding with zeros. */
580 if (nbefore > ndigits)
583 memcpy (out, digits, i);
591 memcpy (out, digits, i);
599 /* Output the decimal point. */
600 *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
602 /* Output leading zeros after the decimal point. */
605 for (i = 0; i < nzero; i++)
609 /* Output digits after the decimal point, padding with zeros. */
612 if (nafter > ndigits)
617 memcpy (out, digits, i);
626 /* Output the exponent. */
635 snprintf (buffer, size, "%+0*d", edigits, e);
637 sprintf (buffer, "%+0*d", edigits, e);
639 memcpy (out, buffer, edigits);
642 if (dtp->u.p.no_leading_blank)
645 memset( out , ' ' , nblanks );
646 dtp->u.p.no_leading_blank = 0;
651 #undef MIN_FIELD_WIDTH
656 /* Write "Infinite" or "Nan" as appropriate for the given format. */
659 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
666 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
668 sign = calculate_sign (dtp, sign_bit);
669 mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
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' */
681 nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
683 p = write_block (dtp, nb);
688 if (unlikely (is_char4_unit (dtp)))
690 gfc_char4_t *p4 = (gfc_char4_t *) p;
691 memset4 (p4, '*', nb);
698 if (unlikely (is_char4_unit (dtp)))
700 gfc_char4_t *p4 = (gfc_char4_t *) p;
701 memset4 (p4, ' ', nb);
710 /* If the sign is negative and the width is 3, there is
711 insufficient room to output '-Inf', so output asterisks */
714 if (unlikely (is_char4_unit (dtp)))
716 gfc_char4_t *p4 = (gfc_char4_t *) p;
717 memset4 (p4, '*', nb);
723 /* The negative sign is mandatory */
727 /* The positive sign is optional, but we output it for
731 if (unlikely (is_char4_unit (dtp)))
733 gfc_char4_t *p4 = (gfc_char4_t *) p;
736 /* We have room, so output 'Infinity' */
737 memcpy4 (p4 + nb - 8, "Infinity", 8);
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);
743 if (sign == S_PLUS || sign == S_MINUS)
745 if (nb < 9 && nb > 3)
746 /* Put the sign in front of Inf */
747 p4[nb - 4] = (gfc_char4_t) fin;
749 /* Put the sign in front of Infinity */
750 p4[nb - 9] = (gfc_char4_t) fin;
756 /* We have room, so output 'Infinity' */
757 memcpy(p + nb - 8, "Infinity", 8);
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);
763 if (sign == S_PLUS || sign == S_MINUS)
765 if (nb < 9 && nb > 3)
766 p[nb - 4] = fin; /* Put the sign in front of Inf */
768 p[nb - 9] = fin; /* Put the sign in front of Infinity */
773 if (unlikely (is_char4_unit (dtp)))
775 gfc_char4_t *p4 = (gfc_char4_t *) p;
776 memcpy4 (p4 + nb - 3, "NaN", 3);
779 memcpy(p + nb - 3, "NaN", 3);
786 /* Returns the value of 10**d. */
788 #define CALCULATE_EXP(x) \
789 inline static GFC_REAL_ ## x \
790 calculate_exp_ ## x (int d)\
793 GFC_REAL_ ## x r = 1.0;\
794 for (i = 0; i< (d >= 0 ? d : -d); i++)\
796 r = (d >= 0) ? r : 1.0 / r;\
804 #ifdef HAVE_GFC_REAL_10
808 #ifdef HAVE_GFC_REAL_16
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:
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]
827 notes: for Gw.d , n' ' means 4 blanks
828 for Gw.dEe, n' ' means e+2 blanks */
830 #define OUTPUT_FLOAT_FMT_G(x) \
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) \
836 int e = f->u.real.e;\
837 int d = f->u.real.d;\
838 int w = f->u.real.w;\
840 GFC_REAL_ ## x rexp_d;\
844 int save_scale_factor, nb = 0;\
847 save_scale_factor = dtp->u.p.scale_factor;\
848 newf = (fnode *) get_mem (sizeof (fnode));\
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)))\
854 newf->format = FMT_E;\
870 GFC_REAL_ ## x temp;\
871 mid = (low + high) / 2;\
873 temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
878 if (ubound == lbound + 1)\
885 if (ubound == lbound + 1)\
906 nb = nb >= w ? 0 : nb;\
907 newf->format = FMT_F;\
908 newf->u.real.w = f->u.real.w - nb;\
911 newf->u.real.d = d - 1;\
913 newf->u.real.d = - (mid - d - 1);\
915 dtp->u.p.scale_factor = 0;\
918 result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
920 dtp->u.p.scale_factor = save_scale_factor;\
924 if (nb > 0 && !dtp->u.p.g0_no_blanks)\
926 p = write_block (dtp, nb);\
929 if (result == FAILURE)\
931 if (unlikely (is_char4_unit (dtp)))\
933 gfc_char4_t *p4 = (gfc_char4_t *) p;\
934 memset4 (p4, pad, nb);\
937 memset (p, pad, nb);\
941 OUTPUT_FLOAT_FMT_G(4)
943 OUTPUT_FLOAT_FMT_G(8)
945 #ifdef HAVE_GFC_REAL_10
946 OUTPUT_FLOAT_FMT_G(10)
949 #ifdef HAVE_GFC_REAL_16
950 OUTPUT_FLOAT_FMT_G(16)
953 #undef OUTPUT_FLOAT_FMT_G
956 /* Define a macro to build code for write_float. */
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.
963 # The result will always contain a decimal point, even if no
966 - The converted value is to be left adjusted on the field boundary
968 + A sign (+ or -) always be placed before a number
970 MIN_FIELD_WIDTH minimum field width
972 * (ndigits-1) is used as the precision
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. */
982 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
983 "e", ndigits - 1, tmp);
986 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
987 "Le", ndigits - 1, tmp);
992 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
993 "e", ndigits - 1, tmp);
996 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
997 "Le", ndigits - 1, tmp);
1001 #if defined(GFC_REAL_16_IS_FLOAT128)
1003 __qmath_(quadmath_snprintf) (buffer, sizeof buffer, \
1004 "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
1005 "Qe", ndigits - 1, tmp);
1008 #define WRITE_FLOAT(x,y)\
1010 GFC_REAL_ ## x tmp;\
1011 tmp = * (GFC_REAL_ ## x *)source;\
1012 sign_bit = signbit (tmp);\
1013 if (!isfinite (tmp))\
1015 write_infnan (dtp, f, isnan (tmp), sign_bit);\
1018 tmp = sign_bit ? -tmp : tmp;\
1019 zero_flag = (tmp == 0.0);\
1023 if (f->format != FMT_G)\
1024 output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
1027 output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
1028 zero_flag, ndigits, edigits);\
1031 /* Output a real number according to its format. */
1034 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
1037 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18
1038 # define MIN_FIELD_WIDTH 46
1040 # define MIN_FIELD_WIDTH 31
1042 #define STR(x) STR1(x)
1045 /* This must be large enough to accurately hold any value. */
1046 char buffer[MIN_FIELD_WIDTH+1];
1047 int sign_bit, ndigits, edigits;
1051 size = MIN_FIELD_WIDTH+1;
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. */
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))
1061 /* Always convert at full precision to avoid double rounding. */
1062 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
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;
1070 ndigits = f->u.real.d;
1071 if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
1072 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1085 #ifdef HAVE_GFC_REAL_10
1090 #ifdef HAVE_GFC_REAL_16
1092 # ifdef GFC_REAL_16_IS_FLOAT128
1100 internal_error (NULL, "bad real kind");