-/* Copyright (C) 2007, 2008 Free Software Foundation, Inc.
+/* Copyright (C) 2007-2015 Free Software Foundation, Inc.
Contributed by Andy Vaught
Write float code factoring to this file by Jerry DeLisle
F2003 I/O support contributed by Jerry DeLisle
-This file is part of the GNU Fortran 95 runtime library (libgfortran).
+This file is part of the GNU Fortran 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)
+the Free Software Foundation; either version 3, 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. */
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
+<http://www.gnu.org/licenses/>. */
#include "config.h"
}
+/* Determine the precision except for EN format. For G format,
+ determines an upper bound to be used for sizing the buffer. */
+
+static int
+determine_precision (st_parameter_dt * dtp, const fnode * f, int len)
+{
+ int precision = f->u.real.d;
+
+ switch (f->format)
+ {
+ case FMT_F:
+ case FMT_G:
+ precision += dtp->u.p.scale_factor;
+ break;
+ case FMT_ES:
+ /* Scale factor has no effect on output. */
+ break;
+ case FMT_E:
+ case FMT_D:
+ /* See F2008 10.7.2.3.3.6 */
+ if (dtp->u.p.scale_factor <= 0)
+ precision += dtp->u.p.scale_factor - 1;
+ break;
+ default:
+ return -1;
+ }
+
+ /* If the scale factor has a large negative value, we must do our
+ own rounding? Use ROUND='NEAREST', which should be what snprintf
+ is using as well. */
+ if (precision < 0 &&
+ (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
+ || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
+ dtp->u.p.current_unit->round_status = ROUND_NEAREST;
+
+ /* Add extra guard digits up to at least full precision when we do
+ our own rounding. */
+ if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
+ && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
+ {
+ precision += 2 * len + 4;
+ if (precision < 0)
+ precision = 0;
+ }
+
+ return precision;
+}
+
+
/* 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, size_t size,
- int sign_bit, bool zero_flag, int ndigits, int edigits)
+static bool
+output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
+ int nprinted, int precision, int sign_bit, bool zero_flag)
{
char *out;
char *digits;
- int e;
- char expchar;
+ int e, w, d, p, i;
+ char expchar, rchar;
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;
+ int ndigits, edigits;
sign_t sign;
ft = f->format;
w = f->u.real.w;
d = f->u.real.d;
+ p = dtp->u.p.scale_factor;
- nzero_real = -1;
+ rchar = '5';
/* We should always know the field width and precision. */
if (d < 0)
sign = calculate_sign (dtp, sign_bit);
- /* The following code checks the given string has punctuation in the correct
- places. Uncomment if needed for debugging.
- if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
- || buffer[ndigits + 2] != 'e'))
- internal_error (&dtp->common, "printf is broken"); */
+ /* Calculate total number of digits. */
+ if (ft == FMT_F)
+ ndigits = nprinted - 2;
+ else
+ ndigits = precision + 1;
/* Read the exponent back in. */
- e = atoi (&buffer[ndigits + 3]) + 1;
+ if (ft != FMT_F)
+ e = atoi (&buffer[ndigits + 3]) + 1;
+ else
+ e = 0;
/* 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);
-
- /* Handle special cases. */
- if (w == 0)
- w = 2;
-
- /* For this one we choose to not output a decimal point.
- F95 10.5.1.2.1 */
- if (w == 1 && ft == FMT_F)
- {
- out = write_block (dtp, w);
- if (out == NULL)
- return;
- *out = '0';
- return;
- }
-
- }
+ e = 0;
/* Normalize the fractional component. */
- buffer[2] = buffer[1];
- digits = &buffer[2];
+ if (ft != FMT_F)
+ {
+ buffer[2] = buffer[1];
+ digits = &buffer[2];
+ }
+ else
+ digits = &buffer[1];
/* Figure out where to place the decimal point. */
switch (ft)
{
case FMT_F:
- nbefore = e + dtp->u.p.scale_factor;
- if (nbefore < 0)
+ nbefore = ndigits - precision;
+ /* Make sure the decimal point is a '.'; depending on the
+ locale, this might not be the case otherwise. */
+ digits[nbefore] = '.';
+ if (p != 0)
{
- nzero = -nbefore;
- nzero_real = nzero;
- if (nzero > d)
- nzero = d;
- nafter = d - nzero;
- nbefore = 0;
+ if (p > 0)
+ {
+
+ memmove (digits + nbefore, digits + nbefore + 1, p);
+ digits[nbefore + p] = '.';
+ nbefore += p;
+ nafter = d - p;
+ if (nafter < 0)
+ nafter = 0;
+ nafter = d;
+ nzero = 0;
+ }
+ else /* p < 0 */
+ {
+ if (nbefore + p >= 0)
+ {
+ nzero = 0;
+ memmove (digits + nbefore + p + 1, digits + nbefore + p, -p);
+ nbefore += p;
+ digits[nbefore] = '.';
+ nafter = d;
+ }
+ else
+ {
+ nzero = -(nbefore + p);
+ memmove (digits + 1, digits, nbefore);
+ digits++;
+ nafter = d + nbefore;
+ nbefore = 0;
+ }
+ if (nzero > d)
+ nzero = d;
+ }
}
else
{
nzero = 0;
nafter = d;
}
+
+ while (digits[0] == '0' && nbefore > 0)
+ {
+ digits++;
+ nbefore--;
+ ndigits--;
+ }
+
expchar = 0;
+ /* If we need to do rounding ourselves, get rid of the dot by
+ moving the fractional part. */
+ if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
+ && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
+ memmove (digits + nbefore, digits + nbefore + 1, ndigits - nbefore);
break;
case FMT_E:
case FMT_D:
i = dtp->u.p.scale_factor;
- if (d <= 0 && i == 0)
+ if (d <= 0 && p == 0)
{
generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
"greater than zero in format specifier 'E' or 'D'");
- return;
+ return false;
}
- if (i <= -d || i >= d + 2)
+ if (p <= -d || p >= d + 2)
{
generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
"out of range in format specifier 'E' or 'D'");
- return;
+ return false;
}
if (!zero_flag)
- e -= i;
- if (i < 0)
+ e -= p;
+ if (p < 0)
{
nbefore = 0;
- nzero = -i;
- nafter = d + i;
+ nzero = -p;
+ nafter = d + p;
}
- else if (i > 0)
+ else if (p > 0)
{
- nbefore = i;
+ nbefore = p;
nzero = 0;
- nafter = (d - i) + 1;
+ nafter = (d - p) + 1;
}
- else /* i == 0 */
+ else /* p == 0 */
{
nbefore = 0;
nzero = 0;
internal_error (&dtp->common, "Unexpected format token");
}
- /* Round the value. */
+ if (zero_flag)
+ goto skip;
+
+ /* Round the value. The value being rounded is an unsigned magnitude. */
+ switch (dtp->u.p.current_unit->round_status)
+ {
+ /* For processor defined and unspecified rounding we use
+ snprintf to print the exact number of digits needed, and thus
+ let snprintf handle the rounding. On system claiming support
+ for IEEE 754, this ought to be round to nearest, ties to
+ even, corresponding to the Fortran ROUND='NEAREST'. */
+ case ROUND_PROCDEFINED:
+ case ROUND_UNSPECIFIED:
+ case ROUND_ZERO: /* Do nothing and truncation occurs. */
+ goto skip;
+ case ROUND_UP:
+ if (sign_bit)
+ goto skip;
+ goto updown;
+ case ROUND_DOWN:
+ if (!sign_bit)
+ goto skip;
+ goto updown;
+ case ROUND_NEAREST:
+ /* Round compatible unless there is a tie. A tie is a 5 with
+ all trailing zero's. */
+ i = nafter + nbefore;
+ if (digits[i] == '5')
+ {
+ for(i++ ; i < ndigits; i++)
+ {
+ if (digits[i] != '0')
+ goto do_rnd;
+ }
+ /* It is a tie so round to even. */
+ switch (digits[nafter + nbefore - 1])
+ {
+ case '1':
+ case '3':
+ case '5':
+ case '7':
+ case '9':
+ /* If odd, round away from zero to even. */
+ break;
+ default:
+ /* If even, skip rounding, truncate to even. */
+ goto skip;
+ }
+ }
+ /* Fall through. */
+ /* The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */
+ case ROUND_COMPATIBLE:
+ rchar = '5';
+ goto do_rnd;
+ }
+
+ updown:
+
+ rchar = '0';
+ if (ft != FMT_F && w > 0 && d == 0 && p == 0)
+ nbefore = 1;
+ /* Scan for trailing zeros to see if we really need to round it. */
+ for(i = nbefore + nafter; i < ndigits; i++)
+ {
+ if (digits[i] != '0')
+ goto do_rnd;
+ }
+ goto skip;
+
+ do_rnd:
+
if (nbefore + nafter == 0)
+ /* Handle the case Fw.0 and value < 1.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;
- }
+ if (digits[0] >= rchar)
+ {
+ /* We rounded to zero but shouldn't have */
+ nbefore = 1;
+ digits--;
+ digits[0] = '1';
+ ndigits = 1;
+ }
}
else if (nbefore + nafter < ndigits)
{
- ndigits = nbefore + nafter;
- i = ndigits;
- if (digits[i] >= '5')
+ i = ndigits = nbefore + nafter;
+ if (digits[i] >= rchar)
{
/* Propagate the carry. */
for (i--; i >= 0; i--)
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. */
+ /* 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)
}
}
+ skip:
+
/* Calculate the format of the exponent field. */
if (expchar)
{
else
edigits = 0;
- /* Zero values always output as positive, even if the value was negative
- before rounding. */
- for (i = 0; i < ndigits; i++)
+ /* Scan the digits string and count the number of zeros. If we make it
+ all the way through the loop, we know the value is zero after the
+ rounding completed above. */
+ int hasdot = 0;
+ for (i = 0; i < ndigits + hasdot; i++)
{
- if (digits[i] != '0')
+ if (digits[i] == '.')
+ hasdot = 1;
+ else if (digits[i] != '0')
break;
}
- if (i == ndigits)
+
+ /* To format properly, we need to know if the rounded result is zero and if
+ so, we set the zero_flag which may have been already set for
+ actual zero. */
+ if (i == ndigits + hasdot)
{
+ zero_flag = true;
/* 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, 0);
}
- /* Pick a field size if none was specified. */
+ /* Pick a field size if none was specified, taking into account small
+ values that may have been rounded to zero. */
if (w <= 0)
- w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+ {
+ if (zero_flag)
+ w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
+ else
+ {
+ w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+ w = w == 1 ? 2 : w;
+ }
+ }
/* Work out how much padding is needed. */
nblanks = w - (nbefore + nzero + nafter + edigits + 1);
/* Create the ouput buffer. */
out = write_block (dtp, w);
if (out == NULL)
- return;
+ return false;
/* Check the value fits in the specified field width. */
- if (nblanks < 0 || edigits == -1)
+ if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
{
+ if (unlikely (is_char4_unit (dtp)))
+ {
+ gfc_char4_t *out4 = (gfc_char4_t *) out;
+ memset4 (out4, '*', w);
+ return false;
+ }
star_fill (out, w);
- return;
+ return false;
}
/* See if we have space for a zero before the decimal point. */
else
leadzero = 0;
+ /* For internal character(kind=4) units, we duplicate the code used for
+ regular output slightly modified. This needs to be maintained
+ consistent with the regular code that follows this block. */
+ if (unlikely (is_char4_unit (dtp)))
+ {
+ gfc_char4_t *out4 = (gfc_char4_t *) out;
+ /* Pad to full field width. */
+
+ if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+ {
+ memset4 (out4, ' ', nblanks);
+ out4 += nblanks;
+ }
+
+ /* Output the initial sign (if any). */
+ if (sign == S_PLUS)
+ *(out4++) = '+';
+ else if (sign == S_MINUS)
+ *(out4++) = '-';
+
+ /* Output an optional leading zero. */
+ if (leadzero)
+ *(out4++) = '0';
+
+ /* Output the part before the decimal point, padding with zeros. */
+ if (nbefore > 0)
+ {
+ if (nbefore > ndigits)
+ {
+ i = ndigits;
+ memcpy4 (out4, digits, i);
+ ndigits = 0;
+ while (i < nbefore)
+ out4[i++] = '0';
+ }
+ else
+ {
+ i = nbefore;
+ memcpy4 (out4, digits, i);
+ ndigits -= i;
+ }
+
+ digits += i;
+ out4 += nbefore;
+ }
+
+ /* Output the decimal point. */
+ *(out4++) = dtp->u.p.current_unit->decimal_status
+ == DECIMAL_POINT ? '.' : ',';
+ if (ft == FMT_F
+ && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
+ || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
+ digits++;
+
+ /* Output leading zeros after the decimal point. */
+ if (nzero > 0)
+ {
+ for (i = 0; i < nzero; i++)
+ *(out4++) = '0';
+ }
+
+ /* Output digits after the decimal point, padding with zeros. */
+ if (nafter > 0)
+ {
+ if (nafter > ndigits)
+ i = ndigits;
+ else
+ i = nafter;
+
+ memcpy4 (out4, digits, i);
+ while (i < nafter)
+ out4[i++] = '0';
+
+ digits += i;
+ ndigits -= i;
+ out4 += nafter;
+ }
+
+ /* Output the exponent. */
+ if (expchar)
+ {
+ if (expchar != ' ')
+ {
+ *(out4++) = expchar;
+ edigits--;
+ }
+ snprintf (buffer, size, "%+0*d", edigits, e);
+ memcpy4 (out4, buffer, edigits);
+ }
+
+ if (dtp->u.p.no_leading_blank)
+ {
+ out4 += edigits;
+ memset4 (out4, ' ' , nblanks);
+ dtp->u.p.no_leading_blank = 0;
+ }
+ return true;
+ } /* End of character(kind=4) internal unit code. */
+
/* Pad to full field width. */
if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
/* Output the decimal point. */
*(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
+ if (ft == FMT_F
+ && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
+ || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
+ digits++;
/* Output leading zeros after the decimal point. */
if (nzero > 0)
*(out++) = expchar;
edigits--;
}
-#if HAVE_SNPRINTF
snprintf (buffer, size, "%+0*d", edigits, e);
-#else
- sprintf (buffer, "%+0*d", edigits, e);
-#endif
memcpy (out, buffer, edigits);
}
dtp->u.p.no_leading_blank = 0;
}
-#undef STR
-#undef STR1
-#undef MIN_FIELD_WIDTH
+ return true;
}
{
char * p, fin;
int nb = 0;
+ sign_t sign;
+ int mark;
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)
+ sign = calculate_sign (dtp, sign_bit);
+ mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
+
+ 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) || dtp->u.p.g0_no_blanks)
+ {
+ if (isnan_flag)
+ nb = 3;
+ else
+ nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
+ }
+ p = write_block (dtp, nb);
+ if (p == NULL)
+ return;
+ if (nb < 3)
+ {
+ if (unlikely (is_char4_unit (dtp)))
{
- memset (p, '*',nb);
- return;
+ gfc_char4_t *p4 = (gfc_char4_t *) p;
+ memset4 (p4, '*', nb);
}
+ else
+ memset (p, '*', nb);
+ return;
+ }
- memset(p, ' ', nb);
- if (!isnan_flag)
+ if (unlikely (is_char4_unit (dtp)))
+ {
+ gfc_char4_t *p4 = (gfc_char4_t *) p;
+ memset4 (p4, ' ', nb);
+ }
+ else
+ memset(p, ' ', nb);
+
+ if (!isnan_flag)
+ {
+ if (sign_bit)
{
- 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 the sign is negative and the width is 3, there is
+ insufficient room to output '-Inf', so output asterisks */
+ if (nb == 3)
+ {
+ if (unlikely (is_char4_unit (dtp)))
+ {
+ gfc_char4_t *p4 = (gfc_char4_t *) p;
+ memset4 (p4, '*', nb);
+ }
+ else
+ memset (p, '*', nb);
+ return;
+ }
+ /* The negative sign is mandatory */
+ fin = '-';
+ }
+ else
+ /* The positive sign is optional, but we output it for
+ consistency */
+ fin = '+';
+
+ if (unlikely (is_char4_unit (dtp)))
+ {
+ gfc_char4_t *p4 = (gfc_char4_t *) p;
- if (nb > 8)
-
- /* We have room, so output 'Infinity' */
- memcpy(p + nb - 8, "Infinity", 8);
+ if (nb > mark)
+ /* We have room, so output 'Infinity' */
+ memcpy4 (p4 + 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);
+ /* For the case of width equals mark, there is not enough room
+ for the sign and 'Infinity' so we go with 'Inf' */
+ memcpy4 (p4 + nb - 3, "Inf", 3);
+
+ if (sign == S_PLUS || sign == S_MINUS)
+ {
+ if (nb < 9 && nb > 3)
+ /* Put the sign in front of Inf */
+ p4[nb - 4] = (gfc_char4_t) fin;
+ else if (nb > 8)
+ /* Put the sign in front of Infinity */
+ p4[nb - 9] = (gfc_char4_t) fin;
+ }
+ return;
+ }
+ if (nb > mark)
+ /* 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 (sign == S_PLUS || sign == S_MINUS)
+ {
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
+ {
+ if (unlikely (is_char4_unit (dtp)))
+ {
+ gfc_char4_t *p4 = (gfc_char4_t *) p;
+ memcpy4 (p4 + nb - 3, "NaN", 3);
+ }
else
memcpy(p + nb - 3, "NaN", 3);
- return;
}
+ return;
}
+}
/* Returns the value of 10**d. */
#define CALCULATE_EXP(x) \
-inline static GFC_REAL_ ## x \
+static GFC_REAL_ ## x \
calculate_exp_ ## x (int d)\
{\
int i;\
#endif
#undef CALCULATE_EXP
+
+/* Define a macro to build code for write_float. */
+
+ /* Note: Before output_float is called, snprintf is used to print to buffer the
+ number in the format +D.DDDDe+ddd.
+
+ # 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
+
+ * prec 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. */
+
+
+#define TOKENPASTE(x, y) TOKENPASTE2(x, y)
+#define TOKENPASTE2(x, y) x ## y
+
+#define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val)
+
+#define DTOA2(prec,val) \
+snprintf (buffer, size, "%+-#.*e", (prec), (val))
+
+#define DTOA2L(prec,val) \
+snprintf (buffer, size, "%+-#.*Le", (prec), (val))
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define DTOA2Q(prec,val) \
+__qmath_(quadmath_snprintf) (buffer, size, "%+-#.*Qe", (prec), (val))
+#endif
+
+#define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val)
+
+/* For F format, we print to the buffer with f format. */
+#define FDTOA2(prec,val) \
+snprintf (buffer, size, "%+-#.*f", (prec), (val))
+
+#define FDTOA2L(prec,val) \
+snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define FDTOA2Q(prec,val) \
+__qmath_(quadmath_snprintf) (buffer, size, "%+-#.*Qf", \
+ (prec), (val))
+#endif
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define ISFINITE2Q(val) finiteq(val)
+#endif
+#define ISFINITE2(val) isfinite(val)
+#define ISFINITE2L(val) isfinite(val)
+
+#define ISFINITE(suff,val) TOKENPASTE(ISFINITE2,suff)(val)
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define SIGNBIT2Q(val) signbitq(val)
+#endif
+#define SIGNBIT2(val) signbit(val)
+#define SIGNBIT2L(val) signbit(val)
+
+#define SIGNBIT(suff,val) TOKENPASTE(SIGNBIT2,suff)(val)
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define ISNAN2Q(val) isnanq(val)
+#endif
+#define ISNAN2(val) isnan(val)
+#define ISNAN2L(val) isnan(val)
+
+#define ISNAN(suff,val) TOKENPASTE(ISNAN2,suff)(val)
+
+
+
/* 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:
m >= 10**d-0.5 Ew.d[Ee]
notes: for Gw.d , n' ' means 4 blanks
- for Gw.dEe, n' ' means e+2 blanks */
+ for Gw.dEe, n' ' means e+2 blanks
+ for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
+ the asm volatile is required for 32-bit x86 platforms. */
-#define OUTPUT_FLOAT_FMT_G(x) \
+#define OUTPUT_FLOAT_FMT_G(x,y) \
static void \
output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
GFC_REAL_ ## x m, char *buffer, size_t size, \
- int sign_bit, bool zero_flag, int ndigits, int edigits) \
+ int sign_bit, bool zero_flag, int comp_d) \
{ \
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;\
+ fnode newf;\
+ GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
int low, high, mid;\
int ubound, lbound;\
- char *p;\
+ char *p, pad = ' ';\
int save_scale_factor, nb = 0;\
+ bool result;\
+ int nprinted, precision;\
+ volatile GFC_REAL_ ## x temp;\
\
save_scale_factor = dtp->u.p.scale_factor;\
- newf = (fnode *) get_mem (sizeof (fnode));\
+\
+ switch (dtp->u.p.current_unit->round_status)\
+ {\
+ case ROUND_ZERO:\
+ r = sign_bit ? 1.0 : 0.0;\
+ break;\
+ case ROUND_UP:\
+ r = 1.0;\
+ break;\
+ case ROUND_DOWN:\
+ r = 0.0;\
+ break;\
+ default:\
+ break;\
+ }\
\
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)))\
+ r_sc = (1 - r / exp_d);\
+ temp = 0.1 * r_sc;\
+ if ((m > 0.0 && ((m < temp) || (r >= (exp_d - m))))\
+ || ((m == 0.0) && !(compile_options.allow_std\
+ & (GFC_STD_F2003 | GFC_STD_F2008)))\
+ || d == 0)\
{ \
- newf->format = FMT_E;\
- newf->u.real.w = w;\
- newf->u.real.d = d;\
- newf->u.real.e = e;\
+ newf.format = FMT_E;\
+ newf.u.real.w = w;\
+ newf.u.real.d = d - comp_d;\
+ newf.u.real.e = e;\
nb = 0;\
+ precision = determine_precision (dtp, &newf, x);\
+ nprinted = DTOA(y,precision,m); \
goto finish;\
}\
\
\
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);\
+ temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
\
if (m < temp)\
{ \
}\
}\
\
- 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);\
-\
+ nb = e <= 0 ? 4 : e + 2;\
+ nb = nb >= w ? w - 1 : nb;\
+ newf.format = FMT_F;\
+ newf.u.real.w = w - nb;\
+ newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
dtp->u.p.scale_factor = 0;\
+ precision = determine_precision (dtp, &newf, x); \
+ nprinted = FDTOA(y,precision,m); \
\
finish:\
- output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
- edigits);\
+ result = output_float (dtp, &newf, buffer, size, nprinted, precision,\
+ sign_bit, zero_flag);\
dtp->u.p.scale_factor = save_scale_factor;\
\
- free_mem(newf);\
\
if (nb > 0 && !dtp->u.p.g0_no_blanks)\
- { \
+ {\
p = write_block (dtp, nb);\
if (p == NULL)\
return;\
- memset (p, ' ', nb);\
+ if (!result)\
+ pad = '*';\
+ if (unlikely (is_char4_unit (dtp)))\
+ {\
+ gfc_char4_t *p4 = (gfc_char4_t *) p;\
+ memset4 (p4, pad, nb);\
+ }\
+ else \
+ memset (p, pad, nb);\
}\
}\
-OUTPUT_FLOAT_FMT_G(4)
+OUTPUT_FLOAT_FMT_G(4,)
-OUTPUT_FLOAT_FMT_G(8)
+OUTPUT_FLOAT_FMT_G(8,)
#ifdef HAVE_GFC_REAL_10
-OUTPUT_FLOAT_FMT_G(10)
+OUTPUT_FLOAT_FMT_G(10,L)
#endif
#ifdef HAVE_GFC_REAL_16
-OUTPUT_FLOAT_FMT_G(16)
+# ifdef GFC_REAL_16_IS_FLOAT128
+OUTPUT_FLOAT_FMT_G(16,Q)
+#else
+OUTPUT_FLOAT_FMT_G(16,L)
+#endif
#endif
#undef OUTPUT_FLOAT_FMT_G
-/* Define a macro to build code for write_float. */
-
- /* Note: Before output_float is called, sprintf is used to print to buffer 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.
-
- # 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
+/* EN format is tricky since the number of significant digits depends
+ on the magnitude. Solve it by first printing a temporary value and
+ figure out the number of significant digits from the printed
+ exponent. Values y, 0.95*10.0**e <= y <10.0**e, are rounded to
+ 10.0**e even when the final result will not be rounded to 10.0**e.
+ For these values the exponent returned by atoi has to be decremented
+ by one. The values y in the ranges
+ (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1))
+ (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2)
+ (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1)
+ are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)),
+ 100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0
+ represents d zeroes, by the lines 279 to 297. */
+
+#define EN_PREC(x,y)\
+{\
+ volatile GFC_REAL_ ## x tmp, one = 1.0;\
+ tmp = * (GFC_REAL_ ## x *)source;\
+ if (ISFINITE (y,tmp))\
+ {\
+ nprinted = DTOA(y,0,tmp);\
+ int e = atoi (&buffer[4]);\
+ if (buffer[1] == '1')\
+ {\
+ tmp = (calculate_exp_ ## x (-e)) * tmp;\
+ tmp = one - (tmp < 0 ? -tmp : tmp); \
+ if (tmp > 0)\
+ e = e - 1;\
+ }\
+ nbefore = e%3;\
+ if (nbefore < 0)\
+ nbefore = 3 + nbefore;\
+ }\
+ else\
+ nprinted = -1;\
+}\
-#define DTOA \
-snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
- "e", ndigits - 1, tmp);
+static int
+determine_en_precision (st_parameter_dt *dtp, const fnode *f,
+ const char *source, int len)
+{
+ int nprinted;
+ char buffer[10];
+ const size_t size = 10;
+ int nbefore; /* digits before decimal point - 1. */
-#define DTOAL \
-snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
- "Le", ndigits - 1, tmp);
+ switch (len)
+ {
+ case 4:
+ EN_PREC(4,)
+ break;
-#else
+ case 8:
+ EN_PREC(8,)
+ break;
-#define DTOA \
-sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
- "e", ndigits - 1, tmp);
+#ifdef HAVE_GFC_REAL_10
+ case 10:
+ EN_PREC(10,L)
+ break;
+#endif
+#ifdef HAVE_GFC_REAL_16
+ case 16:
+# ifdef GFC_REAL_16_IS_FLOAT128
+ EN_PREC(16,Q)
+# else
+ EN_PREC(16,L)
+# endif
+ break;
+#endif
+ default:
+ internal_error (NULL, "bad real kind");
+ }
-#define DTOAL \
-sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
- "Le", ndigits - 1, tmp);
+ if (nprinted == -1)
+ return -1;
-#endif
+ int prec = f->u.real.d + nbefore;
+ if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
+ && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
+ prec += 2 * len + 4;
+ return prec;
+}
+
#define WRITE_FLOAT(x,y)\
{\
GFC_REAL_ ## x tmp;\
tmp = * (GFC_REAL_ ## x *)source;\
- sign_bit = signbit (tmp);\
- if (!isfinite (tmp))\
+ sign_bit = SIGNBIT (y,tmp);\
+ if (!ISFINITE (y,tmp))\
{ \
- write_infnan (dtp, f, isnan (tmp), sign_bit);\
+ write_infnan (dtp, f, ISNAN (y,tmp), sign_bit);\
return;\
}\
tmp = sign_bit ? -tmp : tmp;\
- if (f->u.real.d == 0 && f->format == FMT_F\
- && dtp->u.p.scale_factor == 0)\
- {\
- if (tmp < 0.5)\
- tmp = 0.0;\
- else if (tmp < 1.0)\
- tmp = 1.0;\
- }\
zero_flag = (tmp == 0.0);\
-\
- DTOA ## y\
-\
- if (f->format != FMT_G)\
- output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
- edigits);\
- else \
+ if (f->format == FMT_G)\
output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
- zero_flag, ndigits, edigits);\
+ zero_flag, comp_d);\
+ else\
+ {\
+ if (f->format == FMT_F)\
+ nprinted = FDTOA(y,precision,tmp); \
+ else\
+ nprinted = DTOA(y,precision,tmp); \
+ output_float (dtp, f, buffer, size, nprinted, precision,\
+ sign_bit, zero_flag);\
+ }\
}\
/* Output a real number according to its format. */
static void
-write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
+write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
+ int len, int comp_d)
{
-
-#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;
+ int sign_bit, nprinted;
+ int precision; /* Precision for snprintf call. */
bool zero_flag;
- size_t size;
- size = MIN_FIELD_WIDTH+1;
-
- /* 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;
- }
+ if (f->format != FMT_EN)
+ precision = determine_precision (dtp, f, len);
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;
- }
+ precision = determine_en_precision (dtp, f, source, len);
+
+ /* 4932 is the maximum exponent of long double and quad precision, 3
+ extra characters for the sign, the decimal point, and the
+ trailing null, and finally some extra digits depending on the
+ requested precision. */
+ const size_t size = 4932 + 3 + precision;
+#define BUF_STACK_SZ 5000
+ char buf_stack[BUF_STACK_SZ];
+ char *buffer;
+ if (size > BUF_STACK_SZ)
+ buffer = xmalloc (size);
+ else
+ buffer = buf_stack;
switch (len)
{
#endif
#ifdef HAVE_GFC_REAL_16
case 16:
+# ifdef GFC_REAL_16_IS_FLOAT128
+ WRITE_FLOAT(16,Q)
+# else
WRITE_FLOAT(16,L)
+# endif
break;
#endif
default:
internal_error (NULL, "bad real kind");
}
+ if (size > BUF_STACK_SZ)
+ free (buffer);
}