1090 lines
25 KiB
Modula-2
1090 lines
25 KiB
Modula-2
/* Copyright (C) 2007-2018 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 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 3, or (at your option)
|
|
any later version.
|
|
|
|
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.
|
|
|
|
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"
|
|
|
|
typedef enum
|
|
{ S_NONE, S_MINUS, S_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 = S_NONE;
|
|
|
|
if (negative_flag)
|
|
s = S_MINUS;
|
|
else
|
|
switch (dtp->u.p.sign_status)
|
|
{
|
|
case SIGN_SP: /* Show sign. */
|
|
s = S_PLUS;
|
|
break;
|
|
case SIGN_SS: /* Suppress sign. */
|
|
s = S_NONE;
|
|
break;
|
|
case SIGN_S: /* Processor defined. */
|
|
case SIGN_UNSPECIFIED:
|
|
s = options.optional_plus ? S_PLUS : S_NONE;
|
|
break;
|
|
}
|
|
|
|
return s;
|
|
}
|
|
|
|
|
|
/* 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;
|
|
}
|
|
|
|
|
|
/* Build a real number according to its format which is FMT_G free. */
|
|
|
|
static void
|
|
build_float_string (st_parameter_dt *dtp, const fnode *f, char *buffer,
|
|
size_t size, int nprinted, int precision, int sign_bit,
|
|
bool zero_flag, int npad, char *result, size_t *len)
|
|
{
|
|
char *put;
|
|
char *digits;
|
|
int e, w, d, p, i;
|
|
char expchar, rchar;
|
|
format_token ft;
|
|
/* 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;
|
|
int leadzero;
|
|
int nblanks;
|
|
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;
|
|
|
|
rchar = '5';
|
|
|
|
/* We should always know the field width and precision. */
|
|
if (d < 0)
|
|
internal_error (&dtp->common, "Unspecified precision");
|
|
|
|
sign = calculate_sign (dtp, sign_bit);
|
|
|
|
/* Calculate total number of digits. */
|
|
if (ft == FMT_F)
|
|
ndigits = nprinted - 2;
|
|
else
|
|
ndigits = precision + 1;
|
|
|
|
/* Read the exponent back in. */
|
|
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;
|
|
|
|
/* Normalize the fractional component. */
|
|
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 = ndigits - precision;
|
|
if ((w > 0) && (nbefore > (int) size))
|
|
{
|
|
*len = w;
|
|
star_fill (result, w);
|
|
result[w] = '\0';
|
|
return;
|
|
}
|
|
/* Make sure the decimal point is a '.'; depending on the
|
|
locale, this might not be the case otherwise. */
|
|
digits[nbefore] = '.';
|
|
if (p != 0)
|
|
{
|
|
if (p > 0)
|
|
{
|
|
memmove (digits + nbefore, digits + nbefore + 1, p);
|
|
digits[nbefore + p] = '.';
|
|
nbefore += p;
|
|
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);
|
|
nafter = d - nzero;
|
|
if (nafter == 0 && d > 0)
|
|
{
|
|
/* This is needed to get the correct rounding. */
|
|
memmove (digits + 1, digits, ndigits - 1);
|
|
digits[1] = '0';
|
|
nafter = 1;
|
|
nzero = d - 1;
|
|
}
|
|
else if (nafter < 0)
|
|
{
|
|
/* Reset digits to 0 in order to get correct rounding
|
|
towards infinity. */
|
|
for (i = 0; i < ndigits; i++)
|
|
digits[i] = '0';
|
|
digits[ndigits - 1] = '1';
|
|
nafter = d;
|
|
nzero = 0;
|
|
}
|
|
nbefore = 0;
|
|
}
|
|
}
|
|
}
|
|
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 && p == 0)
|
|
{
|
|
generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
|
|
"greater than zero in format specifier 'E' or 'D'");
|
|
return;
|
|
}
|
|
if (p <= -d || p >= d + 2)
|
|
{
|
|
generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
|
|
"out of range in format specifier 'E' or 'D'");
|
|
return;
|
|
}
|
|
|
|
if (!zero_flag)
|
|
e -= p;
|
|
if (p < 0)
|
|
{
|
|
nbefore = 0;
|
|
nzero = -p;
|
|
nafter = d + p;
|
|
}
|
|
else if (p > 0)
|
|
{
|
|
nbefore = p;
|
|
nzero = 0;
|
|
nafter = (d - p) + 1;
|
|
}
|
|
else /* p == 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");
|
|
}
|
|
|
|
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 (digits[0] >= rchar)
|
|
{
|
|
/* We rounded to zero but shouldn't have */
|
|
nbefore = 1;
|
|
digits--;
|
|
digits[0] = '1';
|
|
ndigits = 1;
|
|
}
|
|
}
|
|
else if (nbefore + nafter < ndigits)
|
|
{
|
|
i = ndigits = nbefore + nafter;
|
|
if (digits[i] >= rchar)
|
|
{
|
|
/* 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++;
|
|
}
|
|
}
|
|
}
|
|
|
|
skip:
|
|
|
|
/* Calculate the format of the exponent field. */
|
|
if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
|
|
{
|
|
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;
|
|
|
|
/* 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] == '.')
|
|
hasdot = 1;
|
|
else if (digits[i] != '0')
|
|
break;
|
|
}
|
|
|
|
/* 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, sign_bit);
|
|
else
|
|
sign = calculate_sign (dtp, 0);
|
|
}
|
|
|
|
/* Pick a field size if none was specified, taking into account small
|
|
values that may have been rounded to zero. */
|
|
if (w <= 0)
|
|
{
|
|
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);
|
|
if (sign != S_NONE)
|
|
nblanks--;
|
|
|
|
/* See if we have space for a zero before the decimal point. */
|
|
if (nbefore == 0 && nblanks > 0)
|
|
{
|
|
leadzero = 1;
|
|
nblanks--;
|
|
}
|
|
else
|
|
leadzero = 0;
|
|
|
|
if (dtp->u.p.g0_no_blanks)
|
|
{
|
|
w -= nblanks;
|
|
nblanks = 0;
|
|
}
|
|
|
|
/* Create the final float string. */
|
|
*len = w + npad;
|
|
put = result;
|
|
|
|
/* Check the value fits in the specified field width. */
|
|
if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
|
|
{
|
|
star_fill (put, *len);
|
|
return;
|
|
}
|
|
|
|
/* Pad to full field width. */
|
|
if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
|
|
{
|
|
memset (put, ' ', nblanks);
|
|
put += nblanks;
|
|
}
|
|
|
|
/* Set the initial sign (if any). */
|
|
if (sign == S_PLUS)
|
|
*(put++) = '+';
|
|
else if (sign == S_MINUS)
|
|
*(put++) = '-';
|
|
|
|
/* Set an optional leading zero. */
|
|
if (leadzero)
|
|
*(put++) = '0';
|
|
|
|
/* Set the part before the decimal point, padding with zeros. */
|
|
if (nbefore > 0)
|
|
{
|
|
if (nbefore > ndigits)
|
|
{
|
|
i = ndigits;
|
|
memcpy (put, digits, i);
|
|
ndigits = 0;
|
|
while (i < nbefore)
|
|
put[i++] = '0';
|
|
}
|
|
else
|
|
{
|
|
i = nbefore;
|
|
memcpy (put, digits, i);
|
|
ndigits -= i;
|
|
}
|
|
|
|
digits += i;
|
|
put += nbefore;
|
|
}
|
|
|
|
/* Set the decimal point. */
|
|
*(put++) = 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++;
|
|
|
|
/* Set leading zeros after the decimal point. */
|
|
if (nzero > 0)
|
|
{
|
|
for (i = 0; i < nzero; i++)
|
|
*(put++) = '0';
|
|
}
|
|
|
|
/* Set digits after the decimal point, padding with zeros. */
|
|
if (nafter > 0)
|
|
{
|
|
if (nafter > ndigits)
|
|
i = ndigits;
|
|
else
|
|
i = nafter;
|
|
|
|
memcpy (put, digits, i);
|
|
while (i < nafter)
|
|
put[i++] = '0';
|
|
|
|
digits += i;
|
|
ndigits -= i;
|
|
put += nafter;
|
|
}
|
|
|
|
/* Set the exponent. */
|
|
if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
|
|
{
|
|
if (expchar != ' ')
|
|
{
|
|
*(put++) = expchar;
|
|
edigits--;
|
|
}
|
|
snprintf (buffer, size, "%+0*d", edigits, e);
|
|
memcpy (put, buffer, edigits);
|
|
put += edigits;
|
|
}
|
|
|
|
if (dtp->u.p.no_leading_blank)
|
|
{
|
|
memset (put , ' ' , nblanks);
|
|
dtp->u.p.no_leading_blank = 0;
|
|
put += nblanks;
|
|
}
|
|
|
|
if (npad > 0 && !dtp->u.p.g0_no_blanks)
|
|
{
|
|
memset (put , ' ' , npad);
|
|
put += npad;
|
|
}
|
|
|
|
/* NULL terminate the string. */
|
|
*put = '\0';
|
|
|
|
return;
|
|
}
|
|
|
|
|
|
/* Write "Infinite" or "Nan" as appropriate for the given format. */
|
|
|
|
static void
|
|
build_infnan_string (st_parameter_dt *dtp, const fnode *f, int isnan_flag,
|
|
int sign_bit, char *p, size_t *len)
|
|
{
|
|
char fin;
|
|
int nb = 0;
|
|
sign_t sign;
|
|
int mark;
|
|
|
|
if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
|
|
{
|
|
sign = calculate_sign (dtp, sign_bit);
|
|
mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
|
|
|
|
nb = f->u.real.w;
|
|
*len = nb;
|
|
|
|
/* 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;
|
|
*len = nb;
|
|
}
|
|
|
|
p[*len] = '\0';
|
|
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 > 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
|
|
memcpy(p + nb - 3, "NaN", 3);
|
|
}
|
|
}
|
|
|
|
|
|
/* Returns the value of 10**d. */
|
|
|
|
#define CALCULATE_EXP(x) \
|
|
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
|
|
|
|
|
|
/* Define macros to build code for format_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) \
|
|
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) \
|
|
quadmath_snprintf (buffer, size, "%+-#.*Qf", \
|
|
(prec), (val))
|
|
#endif
|
|
|
|
|
|
/* 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 (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;\
|
|
}\
|
|
|
|
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. */
|
|
|
|
switch (len)
|
|
{
|
|
case 4:
|
|
EN_PREC(4,)
|
|
break;
|
|
|
|
case 8:
|
|
EN_PREC(8,)
|
|
break;
|
|
|
|
#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");
|
|
}
|
|
|
|
if (nprinted == -1)
|
|
return -1;
|
|
|
|
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;
|
|
}
|
|
|
|
|
|
/* Generate corresponding I/O format. 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
|
|
for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
|
|
the asm volatile is required for 32-bit x86 platforms. */
|
|
#define FORMAT_FLOAT(x,y)\
|
|
{\
|
|
int npad = 0;\
|
|
GFC_REAL_ ## x m;\
|
|
m = * (GFC_REAL_ ## x *)source;\
|
|
sign_bit = signbit (m);\
|
|
if (!isfinite (m))\
|
|
{ \
|
|
build_infnan_string (dtp, f, isnan (m), sign_bit, result, res_len);\
|
|
return;\
|
|
}\
|
|
m = sign_bit ? -m : m;\
|
|
zero_flag = (m == 0.0);\
|
|
if (f->format == FMT_G)\
|
|
{\
|
|
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, r = 0.5, r_sc;\
|
|
int low, high, mid;\
|
|
int ubound, lbound;\
|
|
int save_scale_factor;\
|
|
volatile GFC_REAL_ ## x temp;\
|
|
save_scale_factor = dtp->u.p.scale_factor;\
|
|
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);\
|
|
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 - comp_d;\
|
|
newf.u.real.e = e;\
|
|
npad = 0;\
|
|
precision = determine_precision (dtp, &newf, x);\
|
|
nprinted = DTOA(y,precision,m);\
|
|
}\
|
|
else \
|
|
{\
|
|
mid = 0;\
|
|
low = 0;\
|
|
high = d + 1;\
|
|
lbound = 0;\
|
|
ubound = d + 1;\
|
|
while (low <= high)\
|
|
{\
|
|
mid = (low + high) / 2;\
|
|
temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
|
|
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\
|
|
{\
|
|
mid++;\
|
|
break;\
|
|
}\
|
|
}\
|
|
npad = e <= 0 ? 4 : e + 2;\
|
|
npad = npad >= w ? w - 1 : npad;\
|
|
npad = dtp->u.p.g0_no_blanks ? 0 : npad;\
|
|
newf.format = FMT_F;\
|
|
newf.u.real.w = w - npad;\
|
|
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);\
|
|
}\
|
|
build_float_string (dtp, &newf, buffer, size, nprinted, precision,\
|
|
sign_bit, zero_flag, npad, result, res_len);\
|
|
dtp->u.p.scale_factor = save_scale_factor;\
|
|
}\
|
|
else\
|
|
{\
|
|
if (f->format == FMT_F)\
|
|
nprinted = FDTOA(y,precision,m);\
|
|
else\
|
|
nprinted = DTOA(y,precision,m);\
|
|
build_float_string (dtp, f, buffer, size, nprinted, precision,\
|
|
sign_bit, zero_flag, npad, result, res_len);\
|
|
}\
|
|
}\
|
|
|
|
/* Output a real number according to its format. */
|
|
|
|
|
|
static void
|
|
get_float_string (st_parameter_dt *dtp, const fnode *f, const char *source,
|
|
int kind, int comp_d, char *buffer, int precision,
|
|
size_t size, char *result, size_t *res_len)
|
|
{
|
|
int sign_bit, nprinted;
|
|
bool zero_flag;
|
|
|
|
switch (kind)
|
|
{
|
|
case 4:
|
|
FORMAT_FLOAT(4,)
|
|
break;
|
|
|
|
case 8:
|
|
FORMAT_FLOAT(8,)
|
|
break;
|
|
|
|
#ifdef HAVE_GFC_REAL_10
|
|
case 10:
|
|
FORMAT_FLOAT(10,L)
|
|
break;
|
|
#endif
|
|
#ifdef HAVE_GFC_REAL_16
|
|
case 16:
|
|
# ifdef GFC_REAL_16_IS_FLOAT128
|
|
FORMAT_FLOAT(16,Q)
|
|
# else
|
|
FORMAT_FLOAT(16,L)
|
|
# endif
|
|
break;
|
|
#endif
|
|
default:
|
|
internal_error (NULL, "bad real kind");
|
|
}
|
|
return;
|
|
}
|