1 /* Copyright (C
) 2007, 2008, 2009, 2010 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
, sign_bit
);
117 sign
= calculate_sign (dtp
, 0);
119 /* Handle special cases.
*/
123 /* For this one we choose to not output a decimal point.
125 if (w
== 1 && ft
== FMT_F
)
127 out
= write_block (dtp
, w
);
136 /* Normalize the fractional component.
*/
137 buffer
[2] = buffer
[1];
140 /* Figure out where to place the decimal point.
*/
144 if (d
== 0 && e
<= 0 && dtp
->u.p.scale_factor
== 0)
146 memmove (digits
+ 1, digits
, ndigits
- 1);
151 nbefore
= e
+ dtp
->u.p.scale_factor
;
171 i
= dtp
->u.p.scale_factor
;
172 if (d
<= 0 && i
== 0)
174 generate_error (&dtp
->common
, LIBERROR_FORMAT
, "Precision not "
175 "greater than zero in format specifier 'E' or 'D'");
178 if (i
<= -d || i
>= d
+ 2)
180 generate_error (&dtp
->common
, LIBERROR_FORMAT
, "Scale factor "
181 "out of range in format specifier 'E' or 'D'");
197 nafter
= (d
- i
) + 1;
213 /* The exponent must be a multiple of three
, with
1-3 digits before
214 the decimal point.
*/
223 nbefore
= 3 - nbefore
;
242 /* Should never happen.
*/
243 internal_error (&dtp
->common
, "Unexpected format token");
246 /* Round the value. The value being rounded is an unsigned magnitude.
247 The ROUND_COMPATIBLE is rounding away from zero when there is a tie.
*/
248 switch (dtp
->u.p.current_unit
->round_status
)
250 case ROUND_ZERO
: /* Do nothing and truncation occurs.
*/
263 /* Round compatible unless there is a tie. A tie is a
5 with
264 all trailing zero
's. */
265 i = nafter + nbefore;
266 if (digits[i] == '5')
268 for(i++ ; i < ndigits; i++)
270 if (digits[i] != '0')
273 /* It is a tie so round to even. */
274 switch (digits[nafter + nbefore - 1])
281 /* If odd, round away from zero to even. */
284 /* If even, skip rounding, truncate to even. */
289 case ROUND_PROCDEFINED:
290 case ROUND_UNSPECIFIED:
291 case ROUND_COMPATIBLE:
293 /* Just fall through and do the actual rounding. */
298 if (nbefore + nafter == 0)
301 if (nzero_real == d && digits[0] >= rchar)
303 /* We rounded to zero but shouldn't have
*/
310 else
if (nbefore
+ nafter
< ndigits
)
312 ndigits
= nbefore
+ nafter
;
314 if (digits
[i
] >= rchar
)
316 /* Propagate the carry.
*/
317 for (i
--; i
>= 0; i
--)
319 if (digits
[i
] != '9')
329 /* The carry overflowed. Fortunately we have some spare
330 space at the start of the buffer. We may discard some
331 digits
, but this is ok because we already know they are
345 else
if (ft
== FMT_EN
)
362 /* Calculate the format of the exponent field.
*/
366 for (i
= abs (e
); i
>= 10; i
/= 10)
371 /* Width not specified. Must be no more than
3 digits.
*/
372 if (e
> 999 || e
< -999)
377 if (e
> 99 || e
< -99)
383 /* Exponent width specified
, check it is wide enough.
*/
384 if (edigits
> f
->u.real.e
)
387 edigits
= f
->u.real.e
+ 2;
393 /* Zero values always output as positive
, even if the value was negative
395 for (i
= 0; i
< ndigits
; i
++)
397 if (digits
[i
] != '0')
402 /* The output is zero
, so set the sign according to the sign bit unless
403 -fno
-sign
-zero was specified.
*/
404 if (compile_options.sign_zero
== 1)
405 sign
= calculate_sign (dtp
, sign_bit
);
407 sign
= calculate_sign (dtp
, 0);
410 /* Pick a field size if none was specified.
*/
412 w
= nbefore
+ nzero
+ nafter
+ (sign
!= S_NONE ?
2 : 1);
414 /* Work out how much padding is needed.
*/
415 nblanks
= w
- (nbefore
+ nzero
+ nafter
+ edigits
+ 1);
419 if (dtp
->u.p.g0_no_blanks
)
425 /* Create the ouput buffer.
*/
426 out
= write_block (dtp
, w
);
430 /* Check the value fits in the specified field width.
*/
431 if (nblanks
< 0 || edigits
== -1)
437 /* See if we have space for a zero before the decimal point.
*/
438 if (nbefore
== 0 && nblanks
> 0)
446 /* Pad to full field width.
*/
448 if ( ( nblanks
> 0 ) && !dtp
->u.p.no_leading_blank
)
450 memset (out
, ' ', nblanks
);
454 /* Output the initial
sign (if any
).
*/
457 else
if (sign
== S_MINUS
)
460 /* Output an optional leading zero.
*/
464 /* Output the part before the decimal point
, padding with zeros.
*/
467 if (nbefore
> ndigits
)
470 memcpy (out
, digits
, i
);
478 memcpy (out
, digits
, i
);
486 /* Output the decimal point.
*/
487 *(out
++) = dtp
->u.p.current_unit
->decimal_status
== DECIMAL_POINT ?
'.' : ',';
489 /* Output leading zeros after the decimal point.
*/
492 for (i
= 0; i
< nzero
; i
++)
496 /* Output digits after the decimal point
, padding with zeros.
*/
499 if (nafter
> ndigits
)
504 memcpy (out
, digits
, i
);
513 /* Output the exponent.
*/
522 snprintf (buffer
, size
, "%+0*d", edigits
, e
);
524 sprintf (buffer
, "%+0*d", edigits
, e
);
526 memcpy (out
, buffer
, edigits
);
529 if (dtp
->u.p.no_leading_blank
)
532 memset( out
, ' ' , nblanks
);
533 dtp
->u.p.no_leading_blank
= 0;
538 #undef MIN_FIELD_WIDTH
542 /* Write
"Infinite" or
"Nan" as appropriate for the given format.
*/
545 write_infnan (st_parameter_dt
*dtp
, const fnode
*f
, int isnan_flag
, int sign_bit
)
550 if (f
->format
!= FMT_B
&& f
->format
!= FMT_O
&& f
->format
!= FMT_Z
)
554 /* If the field width is zero
, the processor must select a width
555 not zero.
4 is chosen to allow output of
'-Inf' or
'+Inf' */
558 p
= write_block (dtp
, nb
);
573 /* If the sign is negative and the width is
3, there is
574 insufficient room to output
'-Inf', so output asterisks
*/
582 /* The negative sign is mandatory
*/
588 /* The positive sign is optional
, but we output it for
594 /* We have room
, so output
'Infinity' */
595 memcpy(p
+ nb
- 8, "Infinity", 8);
598 /* For the case of width equals
8, there is not enough room
599 for the sign and
'Infinity' so we go with
'Inf' */
600 memcpy(p
+ nb
- 3, "Inf", 3);
602 if (nb
< 9 && nb
> 3)
603 p
[nb
- 4] = fin
; /* Put the sign in front of Inf
*/
605 p
[nb
- 9] = fin
; /* Put the sign in front of Infinity
*/
608 memcpy(p
+ nb
- 3, "NaN", 3);
614 /* Returns the value of
10**d.
*/
616 #define
CALCULATE_EXP(x
) \
617 inline static GFC_REAL_ ## x \
618 calculate_exp_ ##
x (int d
)\
621 GFC_REAL_ ## x r
= 1.0;\
622 for (i
= 0; i
< (d
>= 0 ? d
: -d
); i
++)\
624 r
= (d
>= 0) ? r
: 1.0 / r
;\
632 #ifdef HAVE_GFC_REAL_10
636 #ifdef HAVE_GFC_REAL_16
641 /* Generate corresponding I
/O format for FMT_G and output.
642 The rules to translate FMT_G to FMT_E or FMT_F from
DEC fortran
643 LRM (table
11-2, Chapter
11, "I/O Formatting", P11
-25) is
:
645 Data Magnitude Equivalent Conversion
646 0< m
< 0.1-0.5*10**(-d
-1) Ew.d
[Ee
]
647 m
= 0 F(w
-n
).
(d
-1), n
' '
648 0.1-0.5*10**(-d
-1)<= m
< 1-0.5*10**(-d
) F(w
-n
).d
, n
' '
649 1-0.5*10**(-d
)<= m
< 10-0.5*10**(-d
+1) F(w
-n
).
(d
-1), n
' '
650 10-0.5*10**(-d
+1)<= m
< 100-0.5*10**(-d
+2) F(w
-n
).
(d
-2), n
' '
651 ................ ..........
652 10**(d
-1)-0.5*10**(-1)<= m
<10**d
-0.5 F(w
-n
).0,n(' ')
653 m
>= 10**d
-0.5 Ew.d
[Ee
]
655 notes
: for Gw.d
, n
' ' means
4 blanks
656 for Gw.dEe
, n
' ' means e
+2 blanks
*/
658 #define
OUTPUT_FLOAT_FMT_G(x
) \
660 output_float_FMT_G_ ##
x (st_parameter_dt
*dtp
, const fnode
*f
, \
661 GFC_REAL_ ## x m
, char
*buffer
, size_t size
, \
662 int sign_bit
, bool zero_flag
, int ndigits
, int edigits
) \
664 int e
= f
->u.real.e
;\
665 int d
= f
->u.real.d
;\
666 int w
= f
->u.real.w
;\
668 GFC_REAL_ ## x rexp_d
;\
672 int save_scale_factor
, nb
= 0;\
674 save_scale_factor
= dtp
->u.p.scale_factor
;\
675 newf
= (fnode *)
get_mem (sizeof (fnode
));\
677 rexp_d
= calculate_exp_ ##
x (-d
);\
678 if ((m
> 0.0 && m
< 0.1 - 0.05 * rexp_d
) ||
(rexp_d
* (m
+ 0.5) >= 1.0) ||\
679 ((m
== 0.0) && !(compile_options.allow_std
& GFC_STD_F2003
)))\
681 newf
->format
= FMT_E
;\
697 GFC_REAL_ ## x temp
;\
698 mid
= (low
+ high
) / 2;\
700 temp
= (calculate_exp_ ##
x (mid
- 1) * (1 - 0.5 * rexp_d
));\
705 if (ubound
== lbound
+ 1)\
712 if (ubound
== lbound
+ 1)\
731 newf
->format
= FMT_F
;\
732 newf
->u.real.w
= f
->u.real.w
- nb
;\
735 newf
->u.real.d
= d
- 1;\
737 newf
->u.real.d
= - (mid
- d
- 1);\
739 dtp
->u.p.scale_factor
= 0;\
742 output_float (dtp
, newf
, buffer
, size
, sign_bit
, zero_flag
, ndigits
, \
744 dtp
->u.p.scale_factor
= save_scale_factor
;\
748 if (nb
> 0 && !dtp
->u.p.g0_no_blanks
)\
750 p
= write_block (dtp
, nb
);\
753 memset (p
, ' ', nb
);\
757 OUTPUT_FLOAT_FMT_G(4)
759 OUTPUT_FLOAT_FMT_G(8)
761 #ifdef HAVE_GFC_REAL_10
762 OUTPUT_FLOAT_FMT_G(10)
765 #ifdef HAVE_GFC_REAL_16
766 OUTPUT_FLOAT_FMT_G(16)
769 #undef OUTPUT_FLOAT_FMT_G
772 /* Define a macro to build code for write_float.
*/
774 /* Note
: Before output_float is called
, sprintf is used to print to buffer the
775 number in the format
+D.DDDDe
+ddd. For an N digit exponent
, this gives us
776 (MIN_FIELD_WIDTH
-5)-N digits after the decimal point
, plus another one
777 before the decimal point.
779 # The result will always contain a decimal point
, even if no
782 - The converted value is to be left adjusted on the field boundary
784 + A
sign (+ or
-) always be placed before a number
786 MIN_FIELD_WIDTH minimum field width
788 * (ndigits
-1) is used as the precision
790 e format
: [-]d.ddde±dd where there is one digit before the
791 decimal
-point character and the number of digits after it is
792 equal to the precision. The exponent always contains at least two
793 digits
; if the value is zero
, the exponent is
00.
*/
798 snprintf (buffer
, size
, "%+-#" STR(MIN_FIELD_WIDTH
) ".*" \
799 "e", ndigits
- 1, tmp
);
802 snprintf (buffer
, size
, "%+-#" STR(MIN_FIELD_WIDTH
) ".*" \
803 "Le", ndigits
- 1, tmp
);
808 sprintf (buffer
, "%+-#" STR(MIN_FIELD_WIDTH
) ".*" \
809 "e", ndigits
- 1, tmp
);
812 sprintf (buffer
, "%+-#" STR(MIN_FIELD_WIDTH
) ".*" \
813 "Le", ndigits
- 1, tmp
);
817 #define
WRITE_FLOAT(x
,y
)\
820 tmp
= * (GFC_REAL_ ## x *)source
;\
821 sign_bit
= __builtin_signbit (tmp
);\
822 if (!isfinite (tmp
))\
824 write_infnan (dtp
, f
, isnan (tmp
), sign_bit
);\
827 tmp
= sign_bit ?
-tmp
: tmp
;\
828 zero_flag
= (tmp
== 0.0);\
832 if (f
->format
!= FMT_G
)\
833 output_float (dtp
, f
, buffer
, size
, sign_bit
, zero_flag
, ndigits
, \
836 output_float_FMT_G_ ##
x (dtp
, f
, tmp
, buffer
, size
, sign_bit
, \
837 zero_flag
, ndigits
, edigits
);\
840 /* Output a real number according to its format.
*/
843 write_float (st_parameter_dt
*dtp
, const fnode
*f
, const char
*source
, int len
)
846 #if
defined(HAVE_GFC_REAL_16
) && __LDBL_DIG__
> 18
847 # define MIN_FIELD_WIDTH
46
849 # define MIN_FIELD_WIDTH
31
851 #define
STR(x
) STR1(x
)
854 /* This must be large enough to accurately hold any value.
*/
855 char buffer
[MIN_FIELD_WIDTH
+1];
856 int sign_bit
, ndigits
, edigits
;
860 size
= MIN_FIELD_WIDTH
+1;
862 /* printf pads blanks for us on the exponent so we just need it big enough
863 to handle the largest number of exponent digits expected.
*/
866 if (f
->format
== FMT_F || f
->format
== FMT_EN || f
->format
== FMT_G
867 ||
((f
->format
== FMT_D || f
->format
== FMT_E
)
868 && dtp
->u.p.scale_factor
!= 0))
870 /* Always convert at full precision to avoid double rounding.
*/
871 ndigits
= MIN_FIELD_WIDTH
- 4 - edigits
;
875 /* The number of digits is known
, so let printf do the rounding.
*/
876 if (f
->format
== FMT_ES
)
877 ndigits
= f
->u.real.d
+ 1;
879 ndigits
= f
->u.real.d
;
880 if (ndigits
> MIN_FIELD_WIDTH
- 4 - edigits
)
881 ndigits
= MIN_FIELD_WIDTH
- 4 - edigits
;
894 #ifdef HAVE_GFC_REAL_10
899 #ifdef HAVE_GFC_REAL_16
905 internal_error (NULL
, "bad real kind");