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 if ((w > 0) && (nbefore > (int) size))
180 star_fill (result, w);
184 /* Make sure the decimal point is a '.'; depending on the
185 locale, this might not be the case otherwise. */
186 digits[nbefore] = '.';
191 memmove (digits + nbefore, digits + nbefore + 1, p);
192 digits[nbefore + p] = '.';
199 if (nbefore + p >= 0)
202 memmove (digits + nbefore + p + 1, digits + nbefore + p, -p);
204 digits[nbefore] = '.';
209 nzero = -(nbefore + p);
210 memmove (digits + 1, digits, nbefore);
212 if (nafter == 0 && d > 0)
214 /* This is needed to get the correct rounding. */
215 memmove (digits + 1, digits, ndigits - 1);
222 /* Reset digits to 0 in order to get correct rounding
224 for (i = 0; i < ndigits; i++)
226 digits[ndigits - 1] = '1';
240 while (digits[0] == '0' && nbefore > 0)
248 /* If we need to do rounding ourselves, get rid of the dot by
249 moving the fractional part. */
250 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
251 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
252 memmove (digits + nbefore, digits + nbefore + 1, ndigits - nbefore);
257 i = dtp->u.p.scale_factor;
258 if (d <= 0 && p == 0)
260 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
261 "greater than zero in format specifier 'E' or 'D'");
264 if (p <= -d || p >= d + 2)
266 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
267 "out of range in format specifier 'E' or 'D'");
283 nafter = (d - p) + 1;
299 /* The exponent must be a multiple of three, with 1-3 digits before
300 the decimal point. */
309 nbefore = 3 - nbefore;
328 /* Should never happen. */
329 internal_error (&dtp->common, "Unexpected format token");
335 /* Round the value. The value being rounded is an unsigned magnitude. */
336 switch (dtp->u.p.current_unit->round_status)
338 /* For processor defined and unspecified rounding we use
339 snprintf to print the exact number of digits needed, and thus
340 let snprintf handle the rounding. On system claiming support
341 for IEEE 754, this ought to be round to nearest, ties to
342 even, corresponding to the Fortran ROUND='NEAREST'. */
343 case ROUND_PROCDEFINED:
344 case ROUND_UNSPECIFIED:
345 case ROUND_ZERO: /* Do nothing and truncation occurs. */
356 /* Round compatible unless there is a tie. A tie is a 5 with
357 all trailing zero's. */
358 i = nafter + nbefore;
359 if (digits[i] == '5')
361 for(i++ ; i < ndigits; i++)
363 if (digits[i] != '0')
366 /* It is a tie so round to even. */
367 switch (digits[nafter + nbefore - 1])
374 /* If odd, round away from zero to even. */
377 /* If even, skip rounding, truncate to even. */
382 /* The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */
383 case ROUND_COMPATIBLE:
391 if (ft != FMT_F && w > 0 && d == 0 && p == 0)
393 /* Scan for trailing zeros to see if we really need to round it. */
394 for(i = nbefore + nafter; i < ndigits; i++)
396 if (digits[i] != '0')
403 if (nbefore + nafter == 0)
404 /* Handle the case Fw.0 and value < 1.0 */
407 if (digits[0] >= rchar)
409 /* We rounded to zero but shouldn't have */
416 else if (nbefore + nafter < ndigits)
418 i = ndigits = nbefore + nafter;
419 if (digits[i] >= rchar)
421 /* Propagate the carry. */
422 for (i--; i >= 0; i--)
424 if (digits[i] != '9')
434 /* The carry overflowed. Fortunately we have some spare
435 space at the start of the buffer. We may discard some
436 digits, but this is ok because we already know they are
450 else if (ft == FMT_EN)
467 /* Calculate the format of the exponent field. */
468 if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
471 for (i = abs (e); i >= 10; i /= 10)
476 /* Width not specified. Must be no more than 3 digits. */
477 if (e > 999 || e < -999)
482 if (e > 99 || e < -99)
488 /* Exponent width specified, check it is wide enough. */
489 if (edigits > f->u.real.e)
492 edigits = f->u.real.e + 2;
498 /* Scan the digits string and count the number of zeros. If we make it
499 all the way through the loop, we know the value is zero after the
500 rounding completed above. */
502 for (i = 0; i < ndigits + hasdot; i++)
504 if (digits[i] == '.')
506 else if (digits[i] != '0')
510 /* To format properly, we need to know if the rounded result is zero and if
511 so, we set the zero_flag which may have been already set for
513 if (i == ndigits + hasdot)
516 /* The output is zero, so set the sign according to the sign bit unless
517 -fno-sign-zero was specified. */
518 if (compile_options.sign_zero == 1)
519 sign = calculate_sign (dtp, sign_bit);
521 sign = calculate_sign (dtp, 0);
524 /* Pick a field size if none was specified, taking into account small
525 values that may have been rounded to zero. */
529 w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
532 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
537 /* Work out how much padding is needed. */
538 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
542 /* See if we have space for a zero before the decimal point. */
543 if (nbefore == 0 && nblanks > 0)
551 if (dtp->u.p.g0_no_blanks)
557 /* Create the final float string. */
561 /* Check the value fits in the specified field width. */
562 if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
564 star_fill (put, *len);
568 /* Pad to full field width. */
569 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
571 memset (put, ' ', nblanks);
575 /* Set the initial sign (if any). */
578 else if (sign == S_MINUS)
581 /* Set an optional leading zero. */
585 /* Set the part before the decimal point, padding with zeros. */
588 if (nbefore > ndigits)
591 memcpy (put, digits, i);
599 memcpy (put, digits, i);
607 /* Set the decimal point. */
608 *(put++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
610 && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
611 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
614 /* Set leading zeros after the decimal point. */
617 for (i = 0; i < nzero; i++)
621 /* Set digits after the decimal point, padding with zeros. */
624 if (nafter > ndigits)
629 memcpy (put, digits, i);
638 /* Set the exponent. */
639 if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
646 snprintf (buffer, size, "%+0*d", edigits, e);
647 memcpy (put, buffer, edigits);
651 if (dtp->u.p.no_leading_blank)
653 memset (put , ' ' , nblanks);
654 dtp->u.p.no_leading_blank = 0;
658 if (npad > 0 && !dtp->u.p.g0_no_blanks)
660 memset (put , ' ' , npad);
664 /* NULL terminate the string. */
671 /* Write "Infinite" or "Nan" as appropriate for the given format. */
674 build_infnan_string (st_parameter_dt *dtp, const fnode *f, int isnan_flag,
675 int sign_bit, char *p, size_t *len)
682 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
684 sign = calculate_sign (dtp, sign_bit);
685 mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
690 /* If the field width is zero, the processor must select a width
691 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
693 if ((nb == 0) || dtp->u.p.g0_no_blanks)
698 nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
715 /* If the sign is negative and the width is 3, there is
716 insufficient room to output '-Inf', so output asterisks */
722 /* The negative sign is mandatory */
726 /* The positive sign is optional, but we output it for
731 /* We have room, so output 'Infinity' */
732 memcpy(p + nb - 8, "Infinity", 8);
734 /* For the case of width equals 8, there is not enough room
735 for the sign and 'Infinity' so we go with 'Inf' */
736 memcpy(p + nb - 3, "Inf", 3);
738 if (sign == S_PLUS || sign == S_MINUS)
740 if (nb < 9 && nb > 3)
741 p[nb - 4] = fin; /* Put the sign in front of Inf */
743 p[nb - 9] = fin; /* Put the sign in front of Infinity */
747 memcpy(p + nb - 3, "NaN", 3);
752 /* Returns the value of 10**d. */
754 #define CALCULATE_EXP(x) \
755 static GFC_REAL_ ## x \
756 calculate_exp_ ## x (int d)\
759 GFC_REAL_ ## x r = 1.0;\
760 for (i = 0; i< (d >= 0 ? d : -d); i++)\
762 r = (d >= 0) ? r : 1.0 / r;\
770 #ifdef HAVE_GFC_REAL_10
774 #ifdef HAVE_GFC_REAL_16
780 /* Define macros to build code for format_float. */
782 /* Note: Before output_float is called, snprintf is used to print to buffer the
783 number in the format +D.DDDDe+ddd.
785 # The result will always contain a decimal point, even if no
788 - The converted value is to be left adjusted on the field boundary
790 + A sign (+ or -) always be placed before a number
792 * prec is used as the precision
794 e format: [-]d.ddde±dd where there is one digit before the
795 decimal-point character and the number of digits after it is
796 equal to the precision. The exponent always contains at least two
797 digits; if the value is zero, the exponent is 00. */
800 #define TOKENPASTE(x, y) TOKENPASTE2(x, y)
801 #define TOKENPASTE2(x, y) x ## y
803 #define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val)
805 #define DTOA2(prec,val) \
806 snprintf (buffer, size, "%+-#.*e", (prec), (val))
808 #define DTOA2L(prec,val) \
809 snprintf (buffer, size, "%+-#.*Le", (prec), (val))
812 #if defined(GFC_REAL_16_IS_FLOAT128)
813 #define DTOA2Q(prec,val) \
814 quadmath_snprintf (buffer, size, "%+-#.*Qe", (prec), (val))
817 #define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val)
819 /* For F format, we print to the buffer with f format. */
820 #define FDTOA2(prec,val) \
821 snprintf (buffer, size, "%+-#.*f", (prec), (val))
823 #define FDTOA2L(prec,val) \
824 snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
827 #if defined(GFC_REAL_16_IS_FLOAT128)
828 #define FDTOA2Q(prec,val) \
829 quadmath_snprintf (buffer, size, "%+-#.*Qf", \
834 /* EN format is tricky since the number of significant digits depends
835 on the magnitude. Solve it by first printing a temporary value and
836 figure out the number of significant digits from the printed
837 exponent. Values y, 0.95*10.0**e <= y <10.0**e, are rounded to
838 10.0**e even when the final result will not be rounded to 10.0**e.
839 For these values the exponent returned by atoi has to be decremented
840 by one. The values y in the ranges
841 (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1))
842 (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2)
843 (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1)
844 are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)),
845 100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0
846 represents d zeroes, by the lines 279 to 297. */
847 #define EN_PREC(x,y)\
849 volatile GFC_REAL_ ## x tmp, one = 1.0;\
850 tmp = * (GFC_REAL_ ## x *)source;\
853 nprinted = DTOA(y,0,tmp);\
854 int e = atoi (&buffer[4]);\
855 if (buffer[1] == '1')\
857 tmp = (calculate_exp_ ## x (-e)) * tmp;\
858 tmp = one - (tmp < 0 ? -tmp : tmp);\
864 nbefore = 3 + nbefore;\
871 determine_en_precision (st_parameter_dt *dtp, const fnode *f,
872 const char *source, int len)
876 const size_t size = 10;
877 int nbefore; /* digits before decimal point - 1. */
889 #ifdef HAVE_GFC_REAL_10
894 #ifdef HAVE_GFC_REAL_16
896 # ifdef GFC_REAL_16_IS_FLOAT128
904 internal_error (NULL, "bad real kind");
910 int prec = f->u.real.d + nbefore;
911 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
912 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
918 /* Generate corresponding I/O format. and output.
919 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
920 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
922 Data Magnitude Equivalent Conversion
923 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
924 m = 0 F(w-n).(d-1), n' '
925 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
926 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
927 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
928 ................ ..........
929 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
930 m >= 10**d-0.5 Ew.d[Ee]
932 notes: for Gw.d , n' ' means 4 blanks
933 for Gw.dEe, n' ' means e+2 blanks
934 for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
935 the asm volatile is required for 32-bit x86 platforms. */
936 #define FORMAT_FLOAT(x,y)\
940 m = * (GFC_REAL_ ## x *)source;\
941 sign_bit = signbit (m);\
944 build_infnan_string (dtp, f, isnan (m), sign_bit, result, res_len);\
947 m = sign_bit ? -m : m;\
948 zero_flag = (m == 0.0);\
949 if (f->format == FMT_G)\
951 int e = f->u.real.e;\
952 int d = f->u.real.d;\
953 int w = f->u.real.w;\
955 GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
958 int save_scale_factor;\
959 volatile GFC_REAL_ ## x temp;\
960 save_scale_factor = dtp->u.p.scale_factor;\
961 switch (dtp->u.p.current_unit->round_status)\
964 r = sign_bit ? 1.0 : 0.0;\
975 exp_d = calculate_exp_ ## x (d);\
976 r_sc = (1 - r / exp_d);\
978 if ((m > 0.0 && ((m < temp) || (r >= (exp_d - m))))\
979 || ((m == 0.0) && !(compile_options.allow_std\
980 & (GFC_STD_F2003 | GFC_STD_F2008)))\
983 newf.format = FMT_E;\
985 newf.u.real.d = d - comp_d;\
988 precision = determine_precision (dtp, &newf, x);\
989 nprinted = DTOA(y,precision,m);\
1000 mid = (low + high) / 2;\
1001 temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
1005 if (ubound == lbound + 1)\
1012 if (ubound == lbound + 1)\
1025 npad = e <= 0 ? 4 : e + 2;\
1026 npad = npad >= w ? w - 1 : npad;\
1027 npad = dtp->u.p.g0_no_blanks ? 0 : npad;\
1028 newf.format = FMT_F;\
1029 newf.u.real.w = w - npad;\
1030 newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
1031 dtp->u.p.scale_factor = 0;\
1032 precision = determine_precision (dtp, &newf, x);\
1033 nprinted = FDTOA(y,precision,m);\
1035 build_float_string (dtp, &newf, buffer, size, nprinted, precision,\
1036 sign_bit, zero_flag, npad, result, res_len);\
1037 dtp->u.p.scale_factor = save_scale_factor;\
1041 if (f->format == FMT_F)\
1042 nprinted = FDTOA(y,precision,m);\
1044 nprinted = DTOA(y,precision,m);\
1045 build_float_string (dtp, f, buffer, size, nprinted, precision,\
1046 sign_bit, zero_flag, npad, result, res_len);\
1050 /* Output a real number according to its format. */
1054 get_float_string (st_parameter_dt *dtp, const fnode *f, const char *source,
1055 int kind, int comp_d, char *buffer, int precision,
1056 size_t size, char *result, size_t *res_len)
1058 int sign_bit, nprinted;
1071 #ifdef HAVE_GFC_REAL_10
1076 #ifdef HAVE_GFC_REAL_16
1078 # ifdef GFC_REAL_16_IS_FLOAT128
1086 internal_error (NULL, "bad real kind");