This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[patch, libfortran] Implement KIND specific floating point formatted output
- From: Jerry DeLisle <jvdelisle at verizon dot net>
- To: Fortran List <fortran at gcc dot gnu dot org>
- Cc: gcc-patches <gcc-patches at gcc dot gnu dot org>
- Date: Tue, 21 Aug 2007 21:26:02 -0700
- Subject: [patch, libfortran] Implement KIND specific floating point formatted output
Hi all,
I was hesitating to submit this because I wanted to make sure there is at least
some benefit.
The attached patch re-factors the formatted write_float and associated functions
to not use GFC_REAL_LARGEST and to use KIND specific functions for each of the
four KINDs supported, 4,8,10, and 16.
I folded most of this into a new file 'write_float.def' which is included by
write.c. I factored as much code as I could that was not KIND dependent into
separate functions and then used macros to build KIND specific functions where
needed.
Performance for KIND=4 is improved some as would be expected. Larger KIND
performance remains about the same. I am concerned that there are some cases
where formatted output is a tad slower. I was using gfortran 4.2 unpatched for
comparison so that may be a difference between 4.2 and 4.3. (hmm is that a
performance regression?).
I have no way to test this on any system that supports KIND=16, though it should
be OK. This patch could be useful for resolving PR23138. Also the factoring
(re-organizing the code) may be useful toward later performance improvements.
Regression tested on x86-64-Gnu-Linux. No new test cases.
OK for trunk.
Regards,
Jerry
2007-08-21 Jerry DeLisle <jvdelisle@gcc.gnu.org>
* io/write.c (stdbool.h): Add include. (sign_t): Move typedef to
new file write_float.def. Include write_float.def.
(extract_real): Delete. (calculate_sign): Delete.
(calculate_exp): Delete. (calculate_G_format): Delete.
(output_float): Delete. (write_float): Delete.
* io/write_float.def (calculate_sign): Added.
(output_float): Refactored to be independent of kind and added to this
file for inclusion. (write_infnan): New function to write "Infinite" or
"NaN" depending on flags passed, independent of kind.
(CALCULATE_EXP): New macro to build kind specific functions. Use it.
(OUTPUT_FLOAT_FMT_G): New macro, likewise. Use it.
(DTOA, DTOAL): Macros to implement "decimal to ascii".
(WRITE_FLOAT): New macro for kind specific write_float functions.
(write_float): Revised function to determine kind and use WRITE_FLOAT
to implement kind specific output.
Index: write.c
===================================================================
--- write.c (revision 127265)
+++ write.c (working copy)
@@ -34,16 +34,13 @@ Boston, MA 02110-1301, USA. */
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
+#include <stdbool.h>
#include "libgfortran.h"
#include "io.h"
#define star_fill(p, n) memset(p, '*', n)
-
-typedef enum
-{ SIGN_NONE, SIGN_MINUS, SIGN_PLUS }
-sign_t;
-
+#include "write_float.def"
void
write_a (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
@@ -235,653 +232,6 @@ extract_uint (const void *p, int len)
return i;
}
-static GFC_REAL_LARGEST
-extract_real (const void *p, int len)
-{
- GFC_REAL_LARGEST i = 0;
- switch (len)
- {
- case 4:
- {
- GFC_REAL_4 tmp;
- memcpy ((void *) &tmp, p, len);
- i = tmp;
- }
- break;
- case 8:
- {
- GFC_REAL_8 tmp;
- memcpy ((void *) &tmp, p, len);
- i = tmp;
- }
- break;
-#ifdef HAVE_GFC_REAL_10
- case 10:
- {
- GFC_REAL_10 tmp;
- memcpy ((void *) &tmp, p, len);
- i = tmp;
- }
- break;
-#endif
-#ifdef HAVE_GFC_REAL_16
- case 16:
- {
- GFC_REAL_16 tmp;
- memcpy ((void *) &tmp, p, len);
- i = tmp;
- }
- break;
-#endif
- default:
- internal_error (NULL, "bad real kind");
- }
- return i;
-}
-
-
-/* Given a flag that indicate if a value is negative or not, return a
- sign_t that gives the sign that we need to produce. */
-
-static sign_t
-calculate_sign (st_parameter_dt *dtp, int negative_flag)
-{
- sign_t s = SIGN_NONE;
-
- if (negative_flag)
- s = SIGN_MINUS;
- else
- switch (dtp->u.p.sign_status)
- {
- case SIGN_SP:
- s = SIGN_PLUS;
- break;
- case SIGN_SS:
- s = SIGN_NONE;
- break;
- case SIGN_S:
- s = options.optional_plus ? SIGN_PLUS : SIGN_NONE;
- break;
- }
-
- return s;
-}
-
-
-/* Returns the value of 10**d. */
-
-static GFC_REAL_LARGEST
-calculate_exp (int d)
-{
- int i;
- GFC_REAL_LARGEST r = 1.0;
-
- for (i = 0; i< (d >= 0 ? d : -d); i++)
- r *= 10;
-
- r = (d >= 0) ? r : 1.0 / r;
-
- return r;
-}
-
-
-/* Generate corresponding I/O format for FMT_G output.
- The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
- LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
-
- Data Magnitude Equivalent Conversion
- 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
- m = 0 F(w-n).(d-1), n' '
- 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
- 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
- 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
- ................ ..........
- 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
- m >= 10**d-0.5 Ew.d[Ee]
-
- notes: for Gw.d , n' ' means 4 blanks
- for Gw.dEe, n' ' means e+2 blanks */
-
-static fnode *
-calculate_G_format (st_parameter_dt *dtp, const fnode *f,
- GFC_REAL_LARGEST value, int *num_blank)
-{
- int e = f->u.real.e;
- int d = f->u.real.d;
- int w = f->u.real.w;
- fnode *newf;
- GFC_REAL_LARGEST m, exp_d;
- int low, high, mid;
- int ubound, lbound;
-
- newf = get_mem (sizeof (fnode));
-
- /* Absolute value. */
- m = (value > 0.0) ? value : -value;
-
- /* In case of the two data magnitude ranges,
- generate E editing, Ew.d[Ee]. */
- exp_d = calculate_exp (d);
- if ((m > 0.0 && m < 0.1 - 0.05 / exp_d) || (m >= exp_d - 0.5 ) ||
- ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))
- {
- newf->format = FMT_E;
- newf->u.real.w = w;
- newf->u.real.d = d;
- newf->u.real.e = e;
- *num_blank = 0;
- return newf;
- }
-
- /* Use binary search to find the data magnitude range. */
- mid = 0;
- low = 0;
- high = d + 1;
- lbound = 0;
- ubound = d + 1;
-
- while (low <= high)
- {
- GFC_REAL_LARGEST temp;
- mid = (low + high) / 2;
-
- /* 0.1 * 10**mid - 0.5 * 10**(mid-d-1) */
- temp = 0.1 * calculate_exp (mid) - 0.5 * calculate_exp (mid - d - 1);
-
- if (m < temp)
- {
- ubound = mid;
- if (ubound == lbound + 1)
- break;
- high = mid - 1;
- }
- else if (m > temp)
- {
- lbound = mid;
- if (ubound == lbound + 1)
- {
- mid ++;
- break;
- }
- low = mid + 1;
- }
- else
- break;
- }
-
- /* Pad with blanks where the exponent would be. */
- if (e < 0)
- *num_blank = 4;
- else
- *num_blank = e + 2;
-
- /* Generate the F editing. F(w-n).(-(mid-d-1)), n' '. */
- newf->format = FMT_F;
- newf->u.real.w = f->u.real.w - *num_blank;
-
- /* Special case. */
- if (m == 0.0)
- newf->u.real.d = d - 1;
- else
- newf->u.real.d = - (mid - d - 1);
-
- /* For F editing, the scale factor is ignored. */
- dtp->u.p.scale_factor = 0;
- return newf;
-}
-
-
-/* Output a real number according to its format which is FMT_G free. */
-
-static void
-output_float (st_parameter_dt *dtp, const fnode *f, GFC_REAL_LARGEST value)
-{
-#if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
-# define MIN_FIELD_WIDTH 46
-#else
-# define MIN_FIELD_WIDTH 31
-#endif
-#define STR(x) STR1(x)
-#define STR1(x) #x
- /* This must be large enough to accurately hold any value. */
- char buffer[MIN_FIELD_WIDTH+1];
- char *out;
- char *digits;
- int e;
- char expchar;
- format_token ft;
- int w;
- int d;
- int edigits;
- int ndigits;
- /* Number of digits before the decimal point. */
- int nbefore;
- /* Number of zeros after the decimal point. */
- int nzero;
- /* Number of digits after the decimal point. */
- int nafter;
- /* Number of zeros after the decimal point, whatever the precision. */
- int nzero_real;
- int leadzero;
- int nblanks;
- int i;
- int sign_bit;
- sign_t sign;
-
- ft = f->format;
- w = f->u.real.w;
- d = f->u.real.d;
-
- nzero_real = -1;
-
-
- /* We should always know the field width and precision. */
- if (d < 0)
- internal_error (&dtp->common, "Unspecified precision");
-
- /* Use sprintf to print the number in the format +D.DDDDe+ddd
- For an N digit exponent, this gives us (MIN_FIELD_WIDTH-5)-N digits
- after the decimal point, plus another one before the decimal point. */
- sign = calculate_sign (dtp, value < 0.0);
- sign_bit = signbit (value);
- if (value < 0)
- value = -value;
-
- /* Special case when format specifies no digits after the decimal point. */
- if (d == 0 && ft == FMT_F)
- {
- if (value < 0.5)
- value = 0.0;
- else if (value < 1.0)
- value = value + 0.5;
- }
-
- /* printf pads blanks for us on the exponent so we just need it big enough
- to handle the largest number of exponent digits expected. */
- edigits=4;
-
- if (ft == FMT_F || ft == FMT_EN
- || ((ft == FMT_D || ft == FMT_E) && dtp->u.p.scale_factor != 0))
- {
- /* Always convert at full precision to avoid double rounding. */
- ndigits = MIN_FIELD_WIDTH - 4 - edigits;
- }
- else
- {
- /* We know the number of digits, so can let printf do the rounding
- for us. */
- if (ft == FMT_ES)
- ndigits = d + 1;
- else
- ndigits = d;
- if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
- ndigits = MIN_FIELD_WIDTH - 4 - edigits;
- }
-
- /* # The result will always contain a decimal point, even if no
- * digits follow it
- *
- * - The converted value is to be left adjusted on the field boundary
- *
- * + A sign (+ or -) always be placed before a number
- *
- * MIN_FIELD_WIDTH minimum field width
- *
- * * (ndigits-1) is used as the precision
- *
- * e format: [-]d.ddde±dd where there is one digit before the
- * decimal-point character and the number of digits after it is
- * equal to the precision. The exponent always contains at least two
- * digits; if the value is zero, the exponent is 00.
- */
-#ifdef HAVE_SNPRINTF
- snprintf (buffer, sizeof (buffer), "%+-#" STR(MIN_FIELD_WIDTH) ".*"
- GFC_REAL_LARGEST_FORMAT "e", ndigits - 1, value);
-#else
- sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*"
- GFC_REAL_LARGEST_FORMAT "e", ndigits - 1, value);
-#endif
-
- /* Check the resulting string has punctuation in the correct places. */
- if (d != 0 && (buffer[2] != '.' || buffer[ndigits + 2] != 'e'))
- internal_error (&dtp->common, "printf is broken");
-
- /* Read the exponent back in. */
- e = atoi (&buffer[ndigits + 3]) + 1;
-
- /* Make sure zero comes out as 0.0e0. */
- if (value == 0.0)
- {
- e = 0;
- if (compile_options.sign_zero == 1)
- sign = calculate_sign (dtp, sign_bit);
- else
- sign = calculate_sign (dtp, 0);
- }
-
- /* Normalize the fractional component. */
- buffer[2] = buffer[1];
- digits = &buffer[2];
-
- /* Figure out where to place the decimal point. */
- switch (ft)
- {
- case FMT_F:
- nbefore = e + dtp->u.p.scale_factor;
- if (nbefore < 0)
- {
- nzero = -nbefore;
- nzero_real = nzero;
- if (nzero > d)
- nzero = d;
- nafter = d - nzero;
- nbefore = 0;
- }
- else
- {
- nzero = 0;
- nafter = d;
- }
- expchar = 0;
- break;
-
- case FMT_E:
- case FMT_D:
- i = dtp->u.p.scale_factor;
- if (value != 0.0)
- e -= i;
- if (i < 0)
- {
- nbefore = 0;
- nzero = -i;
- nafter = d + i;
- }
- else if (i > 0)
- {
- nbefore = i;
- nzero = 0;
- nafter = (d - i) + 1;
- }
- else /* i == 0 */
- {
- nbefore = 0;
- nzero = 0;
- nafter = d;
- }
-
- if (ft == FMT_E)
- expchar = 'E';
- else
- expchar = 'D';
- break;
-
- case FMT_EN:
- /* The exponent must be a multiple of three, with 1-3 digits before
- the decimal point. */
- if (value != 0.0)
- e--;
- if (e >= 0)
- nbefore = e % 3;
- else
- {
- nbefore = (-e) % 3;
- if (nbefore != 0)
- nbefore = 3 - nbefore;
- }
- e -= nbefore;
- nbefore++;
- nzero = 0;
- nafter = d;
- expchar = 'E';
- break;
-
- case FMT_ES:
- if (value != 0.0)
- e--;
- nbefore = 1;
- nzero = 0;
- nafter = d;
- expchar = 'E';
- break;
-
- default:
- /* Should never happen. */
- internal_error (&dtp->common, "Unexpected format token");
- }
-
- /* Round the value. */
- if (nbefore + nafter == 0)
- {
- ndigits = 0;
- if (nzero_real == d && digits[0] >= '5')
- {
- /* We rounded to zero but shouldn't have */
- nzero--;
- nafter = 1;
- digits[0] = '1';
- ndigits = 1;
- }
- }
- else if (nbefore + nafter < ndigits)
- {
- ndigits = nbefore + nafter;
- i = ndigits;
- if (digits[i] >= '5')
- {
- /* Propagate the carry. */
- for (i--; i >= 0; i--)
- {
- if (digits[i] != '9')
- {
- digits[i]++;
- break;
- }
- digits[i] = '0';
- }
-
- if (i < 0)
- {
- /* The carry overflowed. Fortunately we have some spare space
- at the start of the buffer. We may discard some digits, but
- this is ok because we already know they are zero. */
- digits--;
- digits[0] = '1';
- if (ft == FMT_F)
- {
- if (nzero > 0)
- {
- nzero--;
- nafter++;
- }
- else
- nbefore++;
- }
- else if (ft == FMT_EN)
- {
- nbefore++;
- if (nbefore == 4)
- {
- nbefore = 1;
- e += 3;
- }
- }
- else
- e++;
- }
- }
- }
-
- /* Calculate the format of the exponent field. */
- if (expchar)
- {
- edigits = 1;
- for (i = abs (e); i >= 10; i /= 10)
- edigits++;
-
- if (f->u.real.e < 0)
- {
- /* Width not specified. Must be no more than 3 digits. */
- if (e > 999 || e < -999)
- edigits = -1;
- else
- {
- edigits = 4;
- if (e > 99 || e < -99)
- expchar = ' ';
- }
- }
- else
- {
- /* Exponent width specified, check it is wide enough. */
- if (edigits > f->u.real.e)
- edigits = -1;
- else
- edigits = f->u.real.e + 2;
- }
- }
- else
- edigits = 0;
-
- /* Pick a field size if none was specified. */
- if (w <= 0)
- w = nbefore + nzero + nafter + (sign != SIGN_NONE ? 2 : 1);
-
- /* Create the ouput buffer. */
- out = write_block (dtp, w);
- if (out == NULL)
- return;
-
- /* Zero values always output as positive, even if the value was negative
- before rounding. */
- for (i = 0; i < ndigits; i++)
- {
- if (digits[i] != '0')
- break;
- }
- if (i == ndigits)
- {
- /* The output is zero, so set the sign according to the sign bit unless
- -fno-sign-zero was specified. */
- if (compile_options.sign_zero == 1)
- sign = calculate_sign (dtp, sign_bit);
- else
- sign = calculate_sign (dtp, 0);
- }
-
- /* Work out how much padding is needed. */
- nblanks = w - (nbefore + nzero + nafter + edigits + 1);
- if (sign != SIGN_NONE)
- nblanks--;
-
- /* Check the value fits in the specified field width. */
- if (nblanks < 0 || edigits == -1)
- {
- star_fill (out, w);
- return;
- }
-
- /* See if we have space for a zero before the decimal point. */
- if (nbefore == 0 && nblanks > 0)
- {
- leadzero = 1;
- nblanks--;
- }
- else
- leadzero = 0;
-
- /* Pad to full field width. */
-
- if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
- {
- memset (out, ' ', nblanks);
- out += nblanks;
- }
-
- /* Output the initial sign (if any). */
- if (sign == SIGN_PLUS)
- *(out++) = '+';
- else if (sign == SIGN_MINUS)
- *(out++) = '-';
-
- /* Output an optional leading zero. */
- if (leadzero)
- *(out++) = '0';
-
- /* Output the part before the decimal point, padding with zeros. */
- if (nbefore > 0)
- {
- if (nbefore > ndigits)
- {
- i = ndigits;
- memcpy (out, digits, i);
- ndigits = 0;
- while (i < nbefore)
- out[i++] = '0';
- }
- else
- {
- i = nbefore;
- memcpy (out, digits, i);
- ndigits -= i;
- }
-
- digits += i;
- out += nbefore;
- }
- /* Output the decimal point. */
- *(out++) = '.';
-
- /* Output leading zeros after the decimal point. */
- if (nzero > 0)
- {
- for (i = 0; i < nzero; i++)
- *(out++) = '0';
- }
-
- /* Output digits after the decimal point, padding with zeros. */
- if (nafter > 0)
- {
- if (nafter > ndigits)
- i = ndigits;
- else
- i = nafter;
-
- memcpy (out, digits, i);
- while (i < nafter)
- out[i++] = '0';
-
- digits += i;
- ndigits -= i;
- out += nafter;
- }
-
- /* Output the exponent. */
- if (expchar)
- {
- if (expchar != ' ')
- {
- *(out++) = expchar;
- edigits--;
- }
-#if HAVE_SNPRINTF
- snprintf (buffer, sizeof (buffer), "%+0*d", edigits, e);
-#else
- sprintf (buffer, "%+0*d", edigits, e);
-#endif
- memcpy (out, buffer, edigits);
- }
-
- if (dtp->u.p.no_leading_blank)
- {
- out += edigits;
- memset( out , ' ' , nblanks );
- dtp->u.p.no_leading_blank = 0;
- }
-#undef STR
-#undef STR1
-#undef MIN_FIELD_WIDTH
-}
-
void
write_l (st_parameter_dt *dtp, const fnode *f, char *source, int len)
@@ -898,108 +248,6 @@ write_l (st_parameter_dt *dtp, const fno
p[f->u.w - 1] = (n) ? 'T' : 'F';
}
-/* Output a real number according to its format. */
-
-static void
-write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
-{
- GFC_REAL_LARGEST n;
- int nb =0, res, save_scale_factor;
- char * p, fin;
- fnode *f2 = NULL;
-
- n = extract_real (source, len);
-
- if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
- {
- res = isfinite (n);
- if (res == 0)
- {
- nb = f->u.real.w;
-
- /* If the field width is zero, the processor must select a width
- not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
-
- if (nb == 0) nb = 4;
- p = write_block (dtp, nb);
- if (p == NULL)
- return;
- if (nb < 3)
- {
- memset (p, '*',nb);
- return;
- }
-
- memset(p, ' ', nb);
- res = !isnan (n);
- if (res != 0)
- {
- if (signbit(n))
- {
-
- /* If the sign is negative and the width is 3, there is
- insufficient room to output '-Inf', so output asterisks */
-
- if (nb == 3)
- {
- memset (p, '*',nb);
- return;
- }
-
- /* The negative sign is mandatory */
-
- fin = '-';
- }
- else
-
- /* The positive sign is optional, but we output it for
- consistency */
-
- fin = '+';
-
- if (nb > 8)
-
- /* We have room, so output 'Infinity' */
-
- memcpy(p + nb - 8, "Infinity", 8);
- else
-
- /* For the case of width equals 8, there is not enough room
- for the sign and 'Infinity' so we go with 'Inf' */
-
- memcpy(p + nb - 3, "Inf", 3);
- if (nb < 9 && nb > 3)
- p[nb - 4] = fin; /* Put the sign in front of Inf */
- else if (nb > 8)
- p[nb - 9] = fin; /* Put the sign in front of Infinity */
- }
- else
- memcpy(p + nb - 3, "NaN", 3);
- return;
- }
- }
-
- if (f->format != FMT_G)
- output_float (dtp, f, n);
- else
- {
- save_scale_factor = dtp->u.p.scale_factor;
- f2 = calculate_G_format (dtp, f, n, &nb);
- output_float (dtp, f2, n);
- dtp->u.p.scale_factor = save_scale_factor;
- if (f2 != NULL)
- free_mem(f2);
-
- if (nb > 0)
- {
- p = write_block (dtp, nb);
- if (p == NULL)
- return;
- memset (p, ' ', nb);
- }
- }
-}
-
static void
write_int (st_parameter_dt *dtp, const fnode *f, const char *source, int len,
Index: write_float.def
===================================================================
--- write_float.def (revision 0)
+++ write_float.def (revision 0)
@@ -0,0 +1,808 @@
+/* Copyright (C) 2007 Free Software Foundation, Inc.
+ Contributed by Andy Vaught
+ Write float code factoring to this file by Jerry DeLisle
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file. (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Libgfortran; see the file COPYING. If not, write to
+the Free Software Foundation, 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA. */
+
+#include "config.h"
+
+typedef enum
+{ SIGN_NONE, SIGN_MINUS, SIGN_PLUS }
+sign_t;
+
+/* Given a flag that indicates if a value is negative or not, return a
+ sign_t that gives the sign that we need to produce. */
+
+static sign_t
+calculate_sign (st_parameter_dt *dtp, int negative_flag)
+{
+ sign_t s = SIGN_NONE;
+
+ if (negative_flag)
+ s = SIGN_MINUS;
+ else
+ switch (dtp->u.p.sign_status)
+ {
+ case SIGN_SP:
+ s = SIGN_PLUS;
+ break;
+ case SIGN_SS:
+ s = SIGN_NONE;
+ break;
+ case SIGN_S:
+ s = options.optional_plus ? SIGN_PLUS : SIGN_NONE;
+ break;
+ }
+
+ return s;
+}
+
+
+/* Output a real number according to its format which is FMT_G free. */
+
+static void
+output_float (st_parameter_dt *dtp, const fnode *f, char *buffer,
+ int sign_bit, bool zero_flag, int ndigits, int edigits)
+{
+ char *out;
+ char *digits;
+ int e;
+ char expchar;
+ format_token ft;
+ int w;
+ int d;
+ /* Number of digits before the decimal point. */
+ int nbefore;
+ /* Number of zeros after the decimal point. */
+ int nzero;
+ /* Number of digits after the decimal point. */
+ int nafter;
+ /* Number of zeros after the decimal point, whatever the precision. */
+ int nzero_real;
+ int leadzero;
+ int nblanks;
+ int i;
+ sign_t sign;
+
+ ft = f->format;
+ w = f->u.real.w;
+ d = f->u.real.d;
+
+ nzero_real = -1;
+
+ /* We should always know the field width and precision. */
+ if (d < 0)
+ internal_error (&dtp->common, "Unspecified precision");
+
+ /* Use sprintf to print the number in the format +D.DDDDe+ddd
+ For an N digit exponent, this gives us (MIN_FIELD_WIDTH-5)-N digits
+ after the decimal point, plus another one before the decimal point. */
+
+ sign = calculate_sign (dtp, sign_bit);
+
+ /* # The result will always contain a decimal point, even if no
+ * digits follow it
+ *
+ * - The converted value is to be left adjusted on the field boundary
+ *
+ * + A sign (+ or -) always be placed before a number
+ *
+ * MIN_FIELD_WIDTH minimum field width
+ *
+ * * (ndigits-1) is used as the precision
+ *
+ * e format: [-]d.ddde±dd where there is one digit before the
+ * decimal-point character and the number of digits after it is
+ * equal to the precision. The exponent always contains at least two
+ * digits; if the value is zero, the exponent is 00.
+ */
+
+ /* Check the given string has punctuation in the correct places. */
+ if (d != 0 && (buffer[2] != '.' || buffer[ndigits + 2] != 'e'))
+ internal_error (&dtp->common, "printf is broken");
+
+ /* Read the exponent back in. */
+ e = atoi (&buffer[ndigits + 3]) + 1;
+
+ /* Make sure zero comes out as 0.0e0. */
+ if (zero_flag)
+ {
+ e = 0;
+ if (compile_options.sign_zero == 1)
+ sign = calculate_sign (dtp, sign_bit);
+ else
+ sign = calculate_sign (dtp, 0);
+ }
+
+ /* Normalize the fractional component. */
+ buffer[2] = buffer[1];
+ digits = &buffer[2];
+
+ /* Figure out where to place the decimal point. */
+ switch (ft)
+ {
+ case FMT_F:
+ nbefore = e + dtp->u.p.scale_factor;
+ if (nbefore <= 0)
+ {
+ nzero = -nbefore;
+ nzero_real = nzero;
+ if (nzero > d)
+ nzero = d;
+ nafter = d - nzero;
+ nbefore = 0;
+ }
+ else
+ {
+ nzero = 0;
+ nafter = d;
+ }
+ expchar = 0;
+ break;
+
+ case FMT_E:
+ case FMT_D:
+ i = dtp->u.p.scale_factor;
+ if (!zero_flag)
+ e -= i;
+ if (i < 0)
+ {
+ nbefore = 0;
+ nzero = -i;
+ nafter = d + i;
+ }
+ else if (i > 0)
+ {
+ nbefore = i;
+ nzero = 0;
+ nafter = (d - i) + 1;
+ }
+ else /* i == 0 */
+ {
+ nbefore = 0;
+ nzero = 0;
+ nafter = d;
+ }
+
+ if (ft == FMT_E)
+ expchar = 'E';
+ else
+ expchar = 'D';
+ break;
+
+ case FMT_EN:
+ /* The exponent must be a multiple of three, with 1-3 digits before
+ the decimal point. */
+ if (!zero_flag)
+ e--;
+ if (e >= 0)
+ nbefore = e % 3;
+ else
+ {
+ nbefore = (-e) % 3;
+ if (nbefore != 0)
+ nbefore = 3 - nbefore;
+ }
+ e -= nbefore;
+ nbefore++;
+ nzero = 0;
+ nafter = d;
+ expchar = 'E';
+ break;
+
+ case FMT_ES:
+ if (!zero_flag)
+ e--;
+ nbefore = 1;
+ nzero = 0;
+ nafter = d;
+ expchar = 'E';
+ break;
+
+ default:
+ /* Should never happen. */
+ internal_error (&dtp->common, "Unexpected format token");
+ }
+
+ /* Round the value. */
+ if (nbefore + nafter == 0)
+ {
+ ndigits = 0;
+ if (nzero_real == d && digits[0] >= '5')
+ {
+ /* We rounded to zero but shouldn't have */
+ nzero--;
+ nafter = 1;
+ digits[0] = '1';
+ ndigits = 1;
+ }
+ }
+ else if (nbefore + nafter < ndigits)
+ {
+ ndigits = nbefore + nafter;
+ i = ndigits;
+ if (digits[i] >= '5')
+ {
+ /* Propagate the carry. */
+ for (i--; i >= 0; i--)
+ {
+ if (digits[i] != '9')
+ {
+ digits[i]++;
+ break;
+ }
+ digits[i] = '0';
+ }
+
+ if (i < 0)
+ {
+ /* The carry overflowed. Fortunately we have some spare space
+ at the start of the buffer. We may discard some digits, but
+ this is ok because we already know they are zero. */
+ digits--;
+ digits[0] = '1';
+ if (ft == FMT_F)
+ {
+ if (nzero > 0)
+ {
+ nzero--;
+ nafter++;
+ }
+ else
+ nbefore++;
+ }
+ else if (ft == FMT_EN)
+ {
+ nbefore++;
+ if (nbefore == 4)
+ {
+ nbefore = 1;
+ e += 3;
+ }
+ }
+ else
+ e++;
+ }
+ }
+ }
+
+ /* Calculate the format of the exponent field. */
+ if (expchar)
+ {
+ edigits = 1;
+ for (i = abs (e); i >= 10; i /= 10)
+ edigits++;
+
+ if (f->u.real.e < 0)
+ {
+ /* Width not specified. Must be no more than 3 digits. */
+ if (e > 999 || e < -999)
+ edigits = -1;
+ else
+ {
+ edigits = 4;
+ if (e > 99 || e < -99)
+ expchar = ' ';
+ }
+ }
+ else
+ {
+ /* Exponent width specified, check it is wide enough. */
+ if (edigits > f->u.real.e)
+ edigits = -1;
+ else
+ edigits = f->u.real.e + 2;
+ }
+ }
+ else
+ edigits = 0;
+
+ /* Pick a field size if none was specified. */
+ if (w <= 0)
+ w = nbefore + nzero + nafter + (sign != SIGN_NONE ? 2 : 1);
+
+ /* Create the ouput buffer. */
+ out = write_block (dtp, w);
+ if (out == NULL)
+ return;
+
+ /* Zero values always output as positive, even if the value was negative
+ before rounding. */
+ for (i = 0; i < ndigits; i++)
+ {
+ if (digits[i] != '0')
+ break;
+ }
+ if (i == ndigits)
+ {
+ /* The output is zero, so set the sign according to the sign bit unless
+ -fno-sign-zero was specified. */
+ if (compile_options.sign_zero == 1)
+ sign = calculate_sign (dtp, sign_bit);
+ else
+ sign = calculate_sign (dtp, 0);
+ }
+
+ /* Work out how much padding is needed. */
+ nblanks = w - (nbefore + nzero + nafter + edigits + 1);
+ if (sign != SIGN_NONE)
+ nblanks--;
+
+ /* Check the value fits in the specified field width. */
+ if (nblanks < 0 || edigits == -1)
+ {
+ star_fill (out, w);
+ return;
+ }
+
+ /* See if we have space for a zero before the decimal point. */
+ if (nbefore == 0 && nblanks > 0)
+ {
+ leadzero = 1;
+ nblanks--;
+ }
+ else
+ leadzero = 0;
+
+ /* Pad to full field width. */
+
+ if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+ {
+ memset (out, ' ', nblanks);
+ out += nblanks;
+ }
+
+ /* Output the initial sign (if any). */
+ if (sign == SIGN_PLUS)
+ *(out++) = '+';
+ else if (sign == SIGN_MINUS)
+ *(out++) = '-';
+
+ /* Output an optional leading zero. */
+ if (leadzero)
+ *(out++) = '0';
+
+ /* Output the part before the decimal point, padding with zeros. */
+ if (nbefore > 0)
+ {
+ if (nbefore > ndigits)
+ {
+ i = ndigits;
+ memcpy (out, digits, i);
+ ndigits = 0;
+ while (i < nbefore)
+ out[i++] = '0';
+ }
+ else
+ {
+ i = nbefore;
+ memcpy (out, digits, i);
+ ndigits -= i;
+ }
+
+ digits += i;
+ out += nbefore;
+ }
+ /* Output the decimal point. */
+ *(out++) = '.';
+
+ /* Output leading zeros after the decimal point. */
+ if (nzero > 0)
+ {
+ for (i = 0; i < nzero; i++)
+ *(out++) = '0';
+ }
+
+ /* Output digits after the decimal point, padding with zeros. */
+ if (nafter > 0)
+ {
+ if (nafter > ndigits)
+ i = ndigits;
+ else
+ i = nafter;
+
+ memcpy (out, digits, i);
+ while (i < nafter)
+ out[i++] = '0';
+
+ digits += i;
+ ndigits -= i;
+ out += nafter;
+ }
+
+ /* Output the exponent. */
+ if (expchar)
+ {
+ if (expchar != ' ')
+ {
+ *(out++) = expchar;
+ edigits--;
+ }
+#if HAVE_SNPRINTF
+ snprintf (buffer, sizeof (buffer), "%+0*d", edigits, e);
+#else
+ sprintf (buffer, "%+0*d", edigits, e);
+#endif
+ memcpy (out, buffer, edigits);
+ }
+
+ if (dtp->u.p.no_leading_blank)
+ {
+ out += edigits;
+ memset( out , ' ' , nblanks );
+ dtp->u.p.no_leading_blank = 0;
+ }
+#undef STR
+#undef STR1
+#undef MIN_FIELD_WIDTH
+}
+
+
+/* Write "Infinite" or "Nan" as appropriate for the given format. */
+
+static void
+write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
+{
+ char * p, fin;
+ int nb = 0;
+
+ if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
+ {
+ nb = f->u.real.w;
+
+ /* If the field width is zero, the processor must select a width
+ not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
+
+ if (nb == 0) nb = 4;
+ p = write_block (dtp, nb);
+ if (p == NULL)
+ return;
+ if (nb < 3)
+ {
+ memset (p, '*',nb);
+ return;
+ }
+
+ memset(p, ' ', nb);
+ if (!isnan_flag)
+ {
+ if (sign_bit)
+ {
+
+ /* If the sign is negative and the width is 3, there is
+ insufficient room to output '-Inf', so output asterisks */
+
+ if (nb == 3)
+ {
+ memset (p, '*',nb);
+ return;
+ }
+
+ /* The negative sign is mandatory */
+
+ fin = '-';
+ }
+ else
+
+ /* The positive sign is optional, but we output it for
+ consistency */
+ fin = '+';
+
+ if (nb > 8)
+
+ /* We have room, so output 'Infinity' */
+ memcpy(p + nb - 8, "Infinity", 8);
+ else
+
+ /* For the case of width equals 8, there is not enough room
+ for the sign and 'Infinity' so we go with 'Inf' */
+ memcpy(p + nb - 3, "Inf", 3);
+
+ if (nb < 9 && nb > 3)
+ p[nb - 4] = fin; /* Put the sign in front of Inf */
+ else if (nb > 8)
+ p[nb - 9] = fin; /* Put the sign in front of Infinity */
+ }
+ else
+ memcpy(p + nb - 3, "NaN", 3);
+ return;
+ }
+ }
+
+
+/* Returns the value of 10**d. */
+
+#define CALCULATE_EXP(x) \
+inline static GFC_REAL_ ## x \
+calculate_exp_ ## x (int d)\
+{\
+ int i;\
+ GFC_REAL_ ## x r = 1.0;\
+ for (i = 0; i< (d >= 0 ? d : -d); i++)\
+ r *= 10;\
+ r = (d >= 0) ? r : 1.0 / r;\
+ return r;\
+}
+
+CALCULATE_EXP(4)
+
+CALCULATE_EXP(8)
+
+#ifdef HAVE_GFC_REAL_10
+CALCULATE_EXP(10)
+#endif
+
+#ifdef HAVE_GFC_REAL_16
+CALCULATE_EXP(16)
+#endif
+#undef CALCULATE_EXP
+
+/* Generate corresponding I/O format for FMT_G and output.
+ The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
+ LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
+
+ Data Magnitude Equivalent Conversion
+ 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
+ m = 0 F(w-n).(d-1), n' '
+ 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
+ 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
+ 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
+ ................ ..........
+ 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
+ m >= 10**d-0.5 Ew.d[Ee]
+
+ notes: for Gw.d , n' ' means 4 blanks
+ for Gw.dEe, n' ' means e+2 blanks */
+
+#define OUTPUT_FLOAT_FMT_G(x) \
+static void \
+output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
+ GFC_REAL_ ## x m, char *buffer, int sign_bit, \
+ bool zero_flag, int ndigits, int edigits) \
+{ \
+ int e = f->u.real.e;\
+ int d = f->u.real.d;\
+ int w = f->u.real.w;\
+ fnode *newf;\
+ GFC_REAL_ ## x exp_d;\
+ int low, high, mid;\
+ int ubound, lbound;\
+ char *p;\
+ int save_scale_factor, nb = 0;\
+\
+ save_scale_factor = dtp->u.p.scale_factor;\
+ newf = get_mem (sizeof (fnode));\
+\
+ exp_d = calculate_exp_ ## x (d);\
+ if ((m > 0.0 && m < 0.1 - 0.05 / exp_d) || (m >= exp_d - 0.5 ) ||\
+ ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
+ { \
+ newf->format = FMT_E;\
+ newf->u.real.w = w;\
+ newf->u.real.d = d;\
+ newf->u.real.e = e;\
+ nb = 0;\
+ goto finish;\
+ }\
+\
+ mid = 0;\
+ low = 0;\
+ high = d + 1;\
+ lbound = 0;\
+ ubound = d + 1;\
+\
+ while (low <= high)\
+ { \
+ GFC_REAL_ ## x temp;\
+ mid = (low + high) / 2;\
+\
+ temp = 0.1 * calculate_exp_ ## x (mid) - 0.5\
+ * calculate_exp_ ## x (mid - d - 1);\
+\
+ if (m < temp)\
+ { \
+ ubound = mid;\
+ if (ubound == lbound + 1)\
+ break;\
+ high = mid - 1;\
+ }\
+ else if (m > temp)\
+ { \
+ lbound = mid;\
+ if (ubound == lbound + 1)\
+ { \
+ mid ++;\
+ break;\
+ }\
+ low = mid + 1;\
+ }\
+ else\
+ break;\
+ }\
+\
+ if (e < 0)\
+ nb = 4;\
+ else\
+ nb = e + 2;\
+\
+ newf->format = FMT_F;\
+ newf->u.real.w = f->u.real.w - nb;\
+\
+ if (m == 0.0)\
+ newf->u.real.d = d - 1;\
+ else\
+ newf->u.real.d = - (mid - d - 1);\
+\
+ dtp->u.p.scale_factor = 0;\
+\
+ finish:\
+ output_float (dtp, newf, buffer, sign_bit, zero_flag, ndigits, edigits);\
+ dtp->u.p.scale_factor = save_scale_factor;\
+\
+ free_mem(newf);\
+\
+ if (nb > 0)\
+ { \
+ p = write_block (dtp, nb);\
+ if (p == NULL)\
+ return;\
+ memset (p, ' ', nb);\
+ }\
+}\
+
+OUTPUT_FLOAT_FMT_G(4)
+
+OUTPUT_FLOAT_FMT_G(8)
+
+#ifdef HAVE_GFC_REAL_10
+OUTPUT_FLOAT_FMT_G(10)
+#endif
+
+#ifdef HAVE_GFC_REAL_16
+OUTPUT_FLOAT_FMT_G(16)
+#endif
+
+#undef OUTPUT_FLOAT_FMT_G
+
+/* Define a macro to build code for write_float. */
+
+#ifdef HAVE_SNPRINTF
+
+#define DTOA \
+snprintf (buffer, sizeof (buffer), "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
+ "e", ndigits - 1, tmp);
+
+#define DTOAL \
+snprintf (buffer, sizeof (buffer), "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
+ "Le", ndigits - 1, tmp);
+
+#else
+
+#define DTOA \
+sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
+ "e", ndigits - 1, tmp);
+
+#define DTOAL \
+sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
+ "Le", ndigits - 1, tmp);
+
+#endif
+
+#define WRITE_FLOAT(x,y)\
+{\
+ GFC_REAL_ ## x tmp;\
+ tmp = * (GFC_REAL_ ## x *)source;\
+ sign_bit = signbit (tmp);\
+ if (!isfinite (tmp))\
+ { \
+ write_infnan (dtp, f, isnan (tmp), sign_bit);\
+ return;\
+ }\
+ tmp = sign_bit ? -tmp : tmp;\
+ if (f->u.real.d == 0 && f->format == FMT_F)\
+ {\
+ if (tmp < 0.5)\
+ tmp = 0.0;\
+ else if (tmp < 1.0)\
+ tmp = tmp + 0.5;\
+ }\
+ zero_flag = (tmp == 0.0);\
+\
+ DTOA ## y\
+\
+ if (f->format != FMT_G)\
+ output_float (dtp, f, buffer, sign_bit, zero_flag, ndigits, edigits);\
+ else \
+ output_float_FMT_G_ ## x (dtp, f, tmp, buffer, sign_bit, zero_flag, \
+ ndigits, edigits);\
+}\
+
+/* Output a real number according to its format. */
+
+static void
+write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
+{
+
+#if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
+# define MIN_FIELD_WIDTH 46
+#else
+# define MIN_FIELD_WIDTH 31
+#endif
+#define STR(x) STR1(x)
+#define STR1(x) #x
+
+ /* This must be large enough to accurately hold any value. */
+ char buffer[MIN_FIELD_WIDTH+1];
+ int sign_bit, ndigits, edigits;
+ bool zero_flag;
+
+ /* printf pads blanks for us on the exponent so we just need it big enough
+ to handle the largest number of exponent digits expected. */
+ edigits=4;
+
+ if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
+ || ((f->format == FMT_D || f->format == FMT_E)
+ && dtp->u.p.scale_factor != 0))
+ {
+ /* Always convert at full precision to avoid double rounding. */
+ ndigits = MIN_FIELD_WIDTH - 4 - edigits;
+ }
+ else
+ {
+ /* The number of digits is known, so let printf do the rounding. */
+ if (f->format == FMT_ES)
+ ndigits = f->u.real.d + 1;
+ else
+ ndigits = f->u.real.d;
+ if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
+ ndigits = MIN_FIELD_WIDTH - 4 - edigits;
+ }
+
+ switch (len)
+ {
+ case 4:
+ WRITE_FLOAT(4,)
+ break;
+
+ case 8:
+ WRITE_FLOAT(8,)
+ break;
+
+#ifdef HAVE_GFC_REAL_10
+ case 10:
+ WRITE_FLOAT(10,L)
+ break;
+#endif
+#ifdef HAVE_GFC_REAL_16
+ case 16:
+ WRITE_FLOAT(16,L)
+ break;
+#endif
+ default:
+ internal_error (NULL, "bad real kind");
+ }
+}