1 /* Copyright (C
) 2007-2018 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");