1 /* Copyright (C) 2007-2016 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 /* Determine the precision except for EN format. For G format,
63 determines an upper bound to be used for sizing the buffer. */
66 determine_precision (st_parameter_dt * dtp, const fnode * f, int len)
68 int precision = f->u.real.d;
74 precision += dtp->u.p.scale_factor;
77 /* Scale factor has no effect on output. */
81 /* See F2008 10.7.2.3.3.6 */
82 if (dtp->u.p.scale_factor <= 0)
83 precision += dtp->u.p.scale_factor - 1;
89 /* If the scale factor has a large negative value, we must do our
90 own rounding? Use ROUND='NEAREST', which should be what snprintf
93 (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
94 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
95 dtp->u.p.current_unit->round_status = ROUND_NEAREST;
97 /* Add extra guard digits up to at least full precision when we do
99 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
100 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
102 precision += 2 * len + 4;
111 /* Build a real number according to its format which is FMT_G free. */
114 build_float_string (st_parameter_dt *dtp, const fnode *f, char *buffer,
115 size_t size, int nprinted, int precision, int sign_bit,
116 bool zero_flag, int npad, char *result, size_t *len)
123 /* Number of digits before the decimal point. */
125 /* Number of zeros after the decimal point. */
127 /* Number of digits after the decimal point. */
131 int ndigits, edigits;
137 p = dtp->u.p.scale_factor;
141 /* We should always know the field width and precision. */
143 internal_error (&dtp->common, "Unspecified precision");
145 sign = calculate_sign (dtp, sign_bit);
147 /* Calculate total number of digits. */
149 ndigits = nprinted - 2;
151 ndigits = precision + 1;
153 /* Read the exponent back in. */
155 e = atoi (&buffer[ndigits + 3]) + 1;
159 /* Make sure zero comes out as 0.0e0. */
163 /* Normalize the fractional component. */
166 buffer[2] = buffer[1];
172 /* Figure out where to place the decimal point. */
176 nbefore = ndigits - precision;
177 /* Make sure the decimal point is a '.'; depending on the
178 locale, this might not be the case otherwise. */
179 digits[nbefore] = '.';
184 memmove (digits + nbefore, digits + nbefore + 1, p);
185 digits[nbefore + p] = '.';
192 if (nbefore + p >= 0)
195 memmove (digits + nbefore + p + 1, digits + nbefore + p, -p);
197 digits[nbefore] = '.';
202 nzero = -(nbefore + p);
203 memmove (digits + 1, digits, nbefore);
205 if (nafter == 0 && d > 0)
207 /* This is needed to get the correct rounding. */
208 memmove (digits + 1, digits, ndigits - 1);
215 /* Reset digits to 0 in order to get correct rounding
217 for (i = 0; i < ndigits; i++)
219 digits[ndigits - 1] = '1';
233 while (digits[0] == '0' && nbefore > 0)
241 /* If we need to do rounding ourselves, get rid of the dot by
242 moving the fractional part. */
243 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
244 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
245 memmove (digits + nbefore, digits + nbefore + 1, ndigits - nbefore);
250 i = dtp->u.p.scale_factor;
251 if (d <= 0 && p == 0)
253 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
254 "greater than zero in format specifier 'E' or 'D'");
257 if (p <= -d || p >= d + 2)
259 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
260 "out of range in format specifier 'E' or 'D'");
276 nafter = (d - p) + 1;
292 /* The exponent must be a multiple of three, with 1-3 digits before
293 the decimal point. */
302 nbefore = 3 - nbefore;
321 /* Should never happen. */
322 internal_error (&dtp->common, "Unexpected format token");
328 /* Round the value. The value being rounded is an unsigned magnitude. */
329 switch (dtp->u.p.current_unit->round_status)
331 /* For processor defined and unspecified rounding we use
332 snprintf to print the exact number of digits needed, and thus
333 let snprintf handle the rounding. On system claiming support
334 for IEEE 754, this ought to be round to nearest, ties to
335 even, corresponding to the Fortran ROUND='NEAREST'. */
336 case ROUND_PROCDEFINED:
337 case ROUND_UNSPECIFIED:
338 case ROUND_ZERO: /* Do nothing and truncation occurs. */
349 /* Round compatible unless there is a tie. A tie is a 5 with
350 all trailing zero's. */
351 i = nafter + nbefore;
352 if (digits[i] == '5')
354 for(i++ ; i < ndigits; i++)
356 if (digits[i] != '0')
359 /* It is a tie so round to even. */
360 switch (digits[nafter + nbefore - 1])
367 /* If odd, round away from zero to even. */
370 /* If even, skip rounding, truncate to even. */
375 /* The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */
376 case ROUND_COMPATIBLE:
384 if (ft != FMT_F && w > 0 && d == 0 && p == 0)
386 /* Scan for trailing zeros to see if we really need to round it. */
387 for(i = nbefore + nafter; i < ndigits; i++)
389 if (digits[i] != '0')
396 if (nbefore + nafter == 0)
397 /* Handle the case Fw.0 and value < 1.0 */
400 if (digits[0] >= rchar)
402 /* We rounded to zero but shouldn't have */
409 else if (nbefore + nafter < ndigits)
411 i = ndigits = nbefore + nafter;
412 if (digits[i] >= rchar)
414 /* Propagate the carry. */
415 for (i--; i >= 0; i--)
417 if (digits[i] != '9')
427 /* The carry overflowed. Fortunately we have some spare
428 space at the start of the buffer. We may discard some
429 digits, but this is ok because we already know they are
443 else if (ft == FMT_EN)
460 /* Calculate the format of the exponent field. */
461 if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
464 for (i = abs (e); i >= 10; i /= 10)
469 /* Width not specified. Must be no more than 3 digits. */
470 if (e > 999 || e < -999)
475 if (e > 99 || e < -99)
481 /* Exponent width specified, check it is wide enough. */
482 if (edigits > f->u.real.e)
485 edigits = f->u.real.e + 2;
491 /* Scan the digits string and count the number of zeros. If we make it
492 all the way through the loop, we know the value is zero after the
493 rounding completed above. */
495 for (i = 0; i < ndigits + hasdot; i++)
497 if (digits[i] == '.')
499 else if (digits[i] != '0')
503 /* To format properly, we need to know if the rounded result is zero and if
504 so, we set the zero_flag which may have been already set for
506 if (i == ndigits + hasdot)
509 /* The output is zero, so set the sign according to the sign bit unless
510 -fno-sign-zero was specified. */
511 if (compile_options.sign_zero == 1)
512 sign = calculate_sign (dtp, sign_bit);
514 sign = calculate_sign (dtp, 0);
517 /* Pick a field size if none was specified, taking into account small
518 values that may have been rounded to zero. */
522 w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
525 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
530 /* Work out how much padding is needed. */
531 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
535 /* See if we have space for a zero before the decimal point. */
536 if (nbefore == 0 && nblanks > 0)
544 if (dtp->u.p.g0_no_blanks)
550 /* Create the final float string. */
554 /* Check the value fits in the specified field width. */
555 if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
557 star_fill (put, *len);
561 /* Pad to full field width. */
562 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
564 memset (put, ' ', nblanks);
568 /* Set the initial sign (if any). */
571 else if (sign == S_MINUS)
574 /* Set an optional leading zero. */
578 /* Set the part before the decimal point, padding with zeros. */
581 if (nbefore > ndigits)
584 memcpy (put, digits, i);
592 memcpy (put, digits, i);
600 /* Set the decimal point. */
601 *(put++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
603 && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
604 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
607 /* Set leading zeros after the decimal point. */
610 for (i = 0; i < nzero; i++)
614 /* Set digits after the decimal point, padding with zeros. */
617 if (nafter > ndigits)
622 memcpy (put, digits, i);
631 /* Set the exponent. */
632 if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
639 snprintf (buffer, size, "%+0*d", edigits, e);
640 memcpy (put, buffer, edigits);
644 if (dtp->u.p.no_leading_blank)
646 memset (put , ' ' , nblanks);
647 dtp->u.p.no_leading_blank = 0;
651 if (npad > 0 && !dtp->u.p.g0_no_blanks)
653 memset (put , ' ' , npad);
657 /* NULL terminate the string. */
664 /* Write "Infinite" or "Nan" as appropriate for the given format. */
667 build_infnan_string (st_parameter_dt *dtp, const fnode *f, int isnan_flag,
668 int sign_bit, char *p, size_t *len)
675 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
677 sign = calculate_sign (dtp, sign_bit);
678 mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
683 /* If the field width is zero, the processor must select a width
684 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
686 if ((nb == 0) || dtp->u.p.g0_no_blanks)
691 nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
708 /* If the sign is negative and the width is 3, there is
709 insufficient room to output '-Inf', so output asterisks */
715 /* The negative sign is mandatory */
719 /* The positive sign is optional, but we output it for
724 /* We have room, so output 'Infinity' */
725 memcpy(p + nb - 8, "Infinity", 8);
727 /* For the case of width equals 8, there is not enough room
728 for the sign and 'Infinity' so we go with 'Inf' */
729 memcpy(p + nb - 3, "Inf", 3);
731 if (sign == S_PLUS || sign == S_MINUS)
733 if (nb < 9 && nb > 3)
734 p[nb - 4] = fin; /* Put the sign in front of Inf */
736 p[nb - 9] = fin; /* Put the sign in front of Infinity */
740 memcpy(p + nb - 3, "NaN", 3);
745 /* Returns the value of 10**d. */
747 #define CALCULATE_EXP(x) \
748 static GFC_REAL_ ## x \
749 calculate_exp_ ## x (int d)\
752 GFC_REAL_ ## x r = 1.0;\
753 for (i = 0; i< (d >= 0 ? d : -d); i++)\
755 r = (d >= 0) ? r : 1.0 / r;\
763 #ifdef HAVE_GFC_REAL_10
767 #ifdef HAVE_GFC_REAL_16
773 /* Define macros to build code for format_float. */
775 /* Note: Before output_float is called, snprintf is used to print to buffer the
776 number in the format +D.DDDDe+ddd.
778 # The result will always contain a decimal point, even if no
781 - The converted value is to be left adjusted on the field boundary
783 + A sign (+ or -) always be placed before a number
785 * prec is used as the precision
787 e format: [-]d.ddde±dd where there is one digit before the
788 decimal-point character and the number of digits after it is
789 equal to the precision. The exponent always contains at least two
790 digits; if the value is zero, the exponent is 00. */
793 #define TOKENPASTE(x, y) TOKENPASTE2(x, y)
794 #define TOKENPASTE2(x, y) x ## y
796 #define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val)
798 #define DTOA2(prec,val) \
799 snprintf (buffer, size, "%+-#.*e", (prec), (val))
801 #define DTOA2L(prec,val) \
802 snprintf (buffer, size, "%+-#.*Le", (prec), (val))
805 #if defined(GFC_REAL_16_IS_FLOAT128)
806 #define DTOA2Q(prec,val) \
807 quadmath_snprintf (buffer, size, "%+-#.*Qe", (prec), (val))
810 #define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val)
812 /* For F format, we print to the buffer with f format. */
813 #define FDTOA2(prec,val) \
814 snprintf (buffer, size, "%+-#.*f", (prec), (val))
816 #define FDTOA2L(prec,val) \
817 snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
820 #if defined(GFC_REAL_16_IS_FLOAT128)
821 #define FDTOA2Q(prec,val) \
822 quadmath_snprintf (buffer, size, "%+-#.*Qf", \
827 /* EN format is tricky since the number of significant digits depends
828 on the magnitude. Solve it by first printing a temporary value and
829 figure out the number of significant digits from the printed
830 exponent. Values y, 0.95*10.0**e <= y <10.0**e, are rounded to
831 10.0**e even when the final result will not be rounded to 10.0**e.
832 For these values the exponent returned by atoi has to be decremented
833 by one. The values y in the ranges
834 (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1))
835 (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2)
836 (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1)
837 are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)),
838 100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0
839 represents d zeroes, by the lines 279 to 297. */
840 #define EN_PREC(x,y)\
842 volatile GFC_REAL_ ## x tmp, one = 1.0;\
843 tmp = * (GFC_REAL_ ## x *)source;\
846 nprinted = DTOA(y,0,tmp);\
847 int e = atoi (&buffer[4]);\
848 if (buffer[1] == '1')\
850 tmp = (calculate_exp_ ## x (-e)) * tmp;\
851 tmp = one - (tmp < 0 ? -tmp : tmp);\
857 nbefore = 3 + nbefore;\
864 determine_en_precision (st_parameter_dt *dtp, const fnode *f,
865 const char *source, int len)
869 const size_t size = 10;
870 int nbefore; /* digits before decimal point - 1. */
882 #ifdef HAVE_GFC_REAL_10
887 #ifdef HAVE_GFC_REAL_16
889 # ifdef GFC_REAL_16_IS_FLOAT128
897 internal_error (NULL, "bad real kind");
903 int prec = f->u.real.d + nbefore;
904 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
905 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
911 /* Generate corresponding I/O format. and output.
912 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
913 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
915 Data Magnitude Equivalent Conversion
916 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
917 m = 0 F(w-n).(d-1), n' '
918 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
919 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
920 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
921 ................ ..........
922 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
923 m >= 10**d-0.5 Ew.d[Ee]
925 notes: for Gw.d , n' ' means 4 blanks
926 for Gw.dEe, n' ' means e+2 blanks
927 for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
928 the asm volatile is required for 32-bit x86 platforms. */
929 #define FORMAT_FLOAT(x,y)\
933 m = * (GFC_REAL_ ## x *)source;\
934 sign_bit = signbit (m);\
937 build_infnan_string (dtp, f, isnan (m), sign_bit, result, res_len);\
940 m = sign_bit ? -m : m;\
941 zero_flag = (m == 0.0);\
942 if (f->format == FMT_G)\
944 int e = f->u.real.e;\
945 int d = f->u.real.d;\
946 int w = f->u.real.w;\
948 GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
951 int save_scale_factor;\
952 volatile GFC_REAL_ ## x temp;\
953 save_scale_factor = dtp->u.p.scale_factor;\
954 switch (dtp->u.p.current_unit->round_status)\
957 r = sign_bit ? 1.0 : 0.0;\
968 exp_d = calculate_exp_ ## x (d);\
969 r_sc = (1 - r / exp_d);\
971 if ((m > 0.0 && ((m < temp) || (r >= (exp_d - m))))\
972 || ((m == 0.0) && !(compile_options.allow_std\
973 & (GFC_STD_F2003 | GFC_STD_F2008)))\
976 newf.format = FMT_E;\
978 newf.u.real.d = d - comp_d;\
981 precision = determine_precision (dtp, &newf, x);\
982 nprinted = DTOA(y,precision,m);\
993 mid = (low + high) / 2;\
994 temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
998 if (ubound == lbound + 1)\
1005 if (ubound == lbound + 1)\
1018 npad = e <= 0 ? 4 : e + 2;\
1019 npad = npad >= w ? w - 1 : npad;\
1020 npad = dtp->u.p.g0_no_blanks ? 0 : npad;\
1021 newf.format = FMT_F;\
1022 newf.u.real.w = w - npad;\
1023 newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
1024 dtp->u.p.scale_factor = 0;\
1025 precision = determine_precision (dtp, &newf, x);\
1026 nprinted = FDTOA(y,precision,m);\
1028 build_float_string (dtp, &newf, buffer, size, nprinted, precision,\
1029 sign_bit, zero_flag, npad, result, res_len);\
1030 dtp->u.p.scale_factor = save_scale_factor;\
1034 if (f->format == FMT_F)\
1035 nprinted = FDTOA(y,precision,m);\
1037 nprinted = DTOA(y,precision,m);\
1038 build_float_string (dtp, f, buffer, size, nprinted, precision,\
1039 sign_bit, zero_flag, npad, result, res_len);\
1043 /* Output a real number according to its format. */
1047 get_float_string (st_parameter_dt *dtp, const fnode *f, const char *source,
1048 int kind, int comp_d, char *buffer, int precision,
1049 size_t size, char *result, size_t *res_len)
1051 int sign_bit, nprinted;
1064 #ifdef HAVE_GFC_REAL_10
1069 #ifdef HAVE_GFC_REAL_16
1071 # ifdef GFC_REAL_16_IS_FLOAT128
1079 internal_error (NULL, "bad real kind");