1 /* Copyright (C) 2007, 2008, 2009 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 95 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 = 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");