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");