From patchwork Tue Jun 14 04:51:27 2011 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Jerry DeLisle X-Patchwork-Id: 100270 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) by ozlabs.org (Postfix) with SMTP id 94595B6F8F for ; Tue, 14 Jun 2011 14:51:41 +1000 (EST) Received: (qmail 24642 invoked by alias); 14 Jun 2011 04:51:34 -0000 Received: (qmail 24625 invoked by uid 22791); 14 Jun 2011 04:51:29 -0000 X-SWARE-Spam-Status: No, hits=-0.9 required=5.0 tests=AWL, BAYES_05, RCVD_IN_DNSWL_NONE, RFC_ABUSE_POST, TW_CP, TW_XP, T_RP_MATCHES_RCVD X-Spam-Check-By: sourceware.org Received: from mta21.charter.net (HELO mta21.charter.net) (216.33.127.81) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 14 Jun 2011 04:51:12 +0000 Received: from imp09 ([10.20.200.9]) by mta21.charter.net (InterMail vM.7.09.02.04 201-2219-117-106-20090629) with ESMTP id <20110614045111.TEQV11595.mta21.charter.net@imp09>; Tue, 14 Jun 2011 00:51:11 -0400 Received: from quava.localdomain ([68.186.96.128]) by imp09 with smtp.charter.net id vgr91g0072mBQ6U05gr9XF; Tue, 14 Jun 2011 00:51:11 -0400 X-Authority-Analysis: v=1.1 cv=1b2X7W/SifksZeClH/haT1SUt4udqxFGF00pZw2/jJk= c=1 sm=1 a=qd9eegJWHRQA:10 a=tCKlK9Mz1gsA:10 a=yUnIBFQkZM0A:10 a=W3wzxSix39fL1hz8i2isiQ==:17 a=ddUUp38s0aoqahrZGw4A:9 a=wPNLvfGTeEIA:10 a=YHLvdqciOmIWQTWpOEIA:9 a=ZIYK2u4v1T0hRNmCPnIA:7 a=m4AhyEkc1QjQaWJ-:21 a=L7wT3qObOsuBuUDD:21 a=W3wzxSix39fL1hz8i2isiQ==:117 Message-ID: <4DF6E8CF.4070401@charter.net> Date: Mon, 13 Jun 2011 21:51:27 -0700 From: jerry DeLisle User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv:1.9.2.17) Gecko/20110428 Fedora/3.1.10-1.fc15 Thunderbird/3.1.10 MIME-Version: 1.0 To: gfortran CC: Thomas Henlich , Janne Blomqvist , gcc patches Subject: Re: [patch, libgfortran] PR48906 Wrong rounding results with -m32 References: <4DE8D8D6.5030506@charter.net> <4DF25409.1080509@charter.net> In-Reply-To: Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Delivered-To: mailing list gcc-patches@gcc.gnu.org On 06/11/2011 12:23 AM, Thomas Henlich wrote: >> I don't agree with this; with the patch we now output 10 significant >> digits, whereas 9 is sufficient for a binary->ascii->binary roundtrip. >> So please retain the "reduce d by one when E editing is used" thing >> for list format and G0. This is just a side effect of using 1PGw.d >> format for list format and G0 in order to avoid duplicating code, but >> we don't need to follow this particular craziness that is due to how >> the scale factor is specified in the standard. > > I hadn't noticed this, but I agree with Janne. > > It should be easy to implement: > > After the switch between F and E editing, we just need to shift the > decimal point and decrement the exponent. No new rounding is required, > because we keep the number of significant digits. > OK, after a little bit of experimentation, I have arrived at the updated patch attached. This has been regression tested and passes all test cases I am aware of. I also have included a new test case gcc/testsuite/gfortran.dg/fmt_g.f90. OK for trunk? Jerry Index: gcc/testsuite/gfortran.dg/fmt_g.f90 =================================================================== --- gcc/testsuite/gfortran.dg/fmt_g.f90 (revision 0) +++ gcc/testsuite/gfortran.dg/fmt_g.f90 (revision 0) @@ -0,0 +1,88 @@ +! { dg-do run } +! PR48906 Wrong rounding results with -m32 +! This test case presents various corner cases bumped into while reworking +! handling of G formatting of reals without using floating point operations. + implicit none + integer, parameter :: RT = 8 + character(28) :: buf + write(buf, "(rc,f11.2,4x,'u.p.scale_factor; dtp->u.p.scale_factor = 1; set_fnode_default (dtp, &f, length); - write_float (dtp, &f, source , length, 1); + write_float (dtp, &f, source , length); dtp->u.p.scale_factor = org_scale; } @@ -1487,19 +1487,11 @@ void write_real_g0 (st_parameter_dt *dtp, const char *source, int length, int d) { fnode f; - int comp_d; set_fnode_default (dtp, &f, length); if (d > 0) f.u.real.d = d; - - /* Compensate for extra digits when using scale factor, d is not - specified, and the magnitude is such that E editing is used. */ - if (dtp->u.p.scale_factor > 0 && d == 0) - comp_d = 1; - else - comp_d = 0; dtp->u.p.g0_no_blanks = 1; - write_float (dtp, &f, source , length, comp_d); + write_float (dtp, &f, source , length); dtp->u.p.g0_no_blanks = 0; } Index: libgfortran/io/write_float.def =================================================================== --- libgfortran/io/write_float.def (revision 174673) +++ libgfortran/io/write_float.def (working copy) @@ -59,8 +59,22 @@ calculate_sign (st_parameter_dt *dtp, int negative } -/* Output a real number according to its format which is FMT_G free. */ +/* Output a real number according to its format. + What this does: + + - Calculates the number of digits to emit before and after the + decimal point. + - Performs the required rounding. + - Calculates the number of leading blanks to emit. + - Calculates the number of trailing blanks to emit. + - Calculates the exponent format if any. + - Determines if the sign should be displayed. + - Allocates the write buffer according to the width. + - Emits the formatted number according to the above. + + buffer contains the digit string when we enter this function. */ + static try output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, int sign_bit, bool zero_flag, int ndigits, int edigits) @@ -78,8 +92,10 @@ output_float (st_parameter_dt *dtp, const fnode *f int nafter; /* Number of zeros after the decimal point, whatever the precision. */ int nzero_real; + /* Flag for optional leading zero, if there is room. */ int leadzero; - int nblanks; + /* Number of leading and trailing blanks. */ + int lblanks, tblanks; sign_t sign; ft = f->format; @@ -87,7 +103,6 @@ output_float (st_parameter_dt *dtp, const fnode *f d = f->u.real.d; p = dtp->u.p.scale_factor; - rchar = '5'; nzero_real = -1; /* We should always know the field width and precision. */ @@ -142,19 +157,57 @@ output_float (st_parameter_dt *dtp, const fnode *f expchar = 0; break; + case FMT_G: + nbefore = 0; + nzero = 0; + nafter = d; + expchar = 0; + + if (e == d) + { + if (e > 0) + { + nbefore = e; + nafter = d > nbefore ? d - nbefore : 0; + } + else + expchar = 'E'; + break; + } + else if (e < d && e >= 0) + { + nbefore = e; + nafter = d > nbefore ? d - nbefore : 0; + break; + } + else if (e >= -d && e < 0) + { + if (f->u.real.e == -1) + nafter = d; + else if (p == 1) + { + p = 0; + e--; + nbefore++; + nafter--; + expchar = 'E'; + break; + } + } + /* Fall through. */ 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'"); + "greater than zero in format specifier 'D', 'E', " + "or 'G'"); return FAILURE; } if (p <= -d || p >= d + 2) { generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor " - "out of range in format specifier 'E' or 'D'"); + "out of range in format specifier 'E', 'D', or 'G'"); return FAILURE; } @@ -179,10 +232,11 @@ output_float (st_parameter_dt *dtp, const fnode *f nafter = d; } - if (ft == FMT_E) + if (ft == FMT_E || ft == FMT_G) expchar = 'E'; else expchar = 'D'; + break; case FMT_EN: @@ -224,6 +278,7 @@ output_float (st_parameter_dt *dtp, const fnode *f /* Round the value. The value being rounded is an unsigned magnitude. The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */ + rchar = '5'; switch (dtp->u.p.current_unit->round_status) { case ROUND_ZERO: /* Do nothing and truncation occurs. */ @@ -275,6 +330,7 @@ output_float (st_parameter_dt *dtp, const fnode *f rchar = '0'; if (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++) { @@ -284,7 +340,7 @@ output_float (st_parameter_dt *dtp, const fnode *f goto skip; do_rnd: - + if (nbefore + nafter == 0) { ndigits = 0; @@ -321,6 +377,7 @@ output_float (st_parameter_dt *dtp, const fnode *f zero. */ digits--; digits[0] = '1'; + if (ft == FMT_F) { if (nzero > 0) @@ -331,6 +388,27 @@ output_float (st_parameter_dt *dtp, const fnode *f else nbefore++; } + else if (ft == FMT_G) + { + if (e < d && e >= 0) + { + nbefore++; + nafter--; + } + else if (e == d) + { + nbefore -= d; + nafter += d; + e++; + expchar = 'E'; + } + else + { + e++; + if (e == 0) + expchar = 0; + } + } else if (ft == FMT_EN) { nbefore++; @@ -346,11 +424,68 @@ output_float (st_parameter_dt *dtp, const fnode *f } } + /* 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. 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. */ + for (i = 0; i < ndigits; i++) + { + if (digits[i] != '0') + break; + } + if (i == ndigits) + { + 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); + } + skip: + /* 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; + } + } + + /* Figure out the leading and trailing blanks. */ + lblanks = tblanks = 0; + + /* G formatting is a special case of several things. */ + if (ft == FMT_G) + { + if (zero_flag && nbefore == 0 && nafter > 0) + { + nbefore++; + nafter--; + } + tblanks = edigits = f->u.real.e == -1 ? 4 : f->u.real.e + 2; + lblanks = w - (nbefore + nzero + nafter + 1 + tblanks + + (sign != S_NONE ? 1 : 0)); + } + else + { + lblanks = w - (nbefore + nzero + nafter + 1 + e); + if (sign != S_NONE) + lblanks--; + } + /* Calculate the format of the exponent field. */ if (expchar) { + tblanks = 0; edigits = 1; for (i = abs (e); i >= 10; i /= 10) edigits++; @@ -379,76 +514,57 @@ output_float (st_parameter_dt *dtp, const fnode *f 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. */ - for (i = 0; i < ndigits; i++) + if (ft == FMT_G) { - 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) - { - 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 + if (dtp->u.p.g0_no_blanks) { - w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1); - w = w == 1 ? 2 : w; + w = w - tblanks - lblanks; + lblanks = tblanks = 0; } + else if (expchar == 0) + { + lblanks = w - (nbefore + nzero + nafter + 1 + tblanks + + (sign != S_NONE ? 1 : 0)); + if (dtp->u.p.no_leading_blank) + { + tblanks += lblanks; + lblanks = 0; + } + } } - - /* Work out how much padding is needed. */ - nblanks = w - (nbefore + nzero + nafter + edigits + 1); - if (sign != S_NONE) - nblanks--; + else + lblanks = w - (nbefore + nzero + nafter + 1 + edigits + + (sign != S_NONE ? 1 : 0)); - if (dtp->u.p.g0_no_blanks) - { - w -= nblanks; - nblanks = 0; - } - /* Create the ouput buffer. */ - out = write_block (dtp, w); + if ((ft == FMT_G && dtp->u.p.g0_no_blanks) + || (ft == FMT_F && f->u.real.w == 0)) + out = write_block (dtp, w); + else + out = write_block (dtp, f->u.real.w); + if (out == NULL) return FAILURE; /* Check the value fits in the specified field width. */ - if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE)) + if (lblanks < 0 || tblanks < 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); + memset4 (out4, '*', f->u.real.w); return FAILURE; } - star_fill (out, w); + star_fill (out, f->u.real.w); return FAILURE; } /* See if we have space for a zero before the decimal point. */ - if (nbefore == 0 && nblanks > 0) + if (nbefore == 0 && lblanks > 0) { leadzero = 1; - nblanks--; + lblanks--; } else leadzero = 0; @@ -461,10 +577,10 @@ output_float (st_parameter_dt *dtp, const fnode *f gfc_char4_t *out4 = (gfc_char4_t *) out; /* Pad to full field width. */ - if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank) + if ((lblanks > 0) && !dtp->u.p.no_leading_blank) { - memset4 (out4, ' ', nblanks); - out4 += nblanks; + memset4 (out4, ' ', lblanks); + out4 += lblanks; } /* Output the initial sign (if any). */ @@ -539,21 +655,17 @@ output_float (st_parameter_dt *dtp, const fnode *f memcpy4 (out4, buffer, edigits); } - if (dtp->u.p.no_leading_blank) - { - out4 += edigits; - memset4 (out4, ' ' , nblanks); - dtp->u.p.no_leading_blank = 0; - } + /* Emit trailing blanks if needed. */ + memset4 (out4, ' ' , tblanks); return SUCCESS; } /* End of character(kind=4) internal unit code. */ - /* Pad to full field width. */ + /* Pad leading to full field width. */ - if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank) + if (lblanks > 0 && !dtp->u.p.no_leading_blank) { - memset (out, ' ', nblanks); - out += nblanks; + memset (out, ' ', lblanks); + out += lblanks; } /* Output the initial sign (if any). */ @@ -627,13 +739,10 @@ output_float (st_parameter_dt *dtp, const fnode *f memcpy (out, buffer, edigits); } - if (dtp->u.p.no_leading_blank) - { - out += edigits; - memset( out , ' ' , nblanks ); - dtp->u.p.no_leading_blank = 0; - } - + /* Emit trailing blanks if needed. */ + if (tblanks) + memset( out , ' ' , tblanks ); + #undef STR #undef STR1 #undef MIN_FIELD_WIDTH @@ -770,185 +879,6 @@ write_infnan (st_parameter_dt *dtp, const fnode *f } } - -/* 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 - 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) \ -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 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 rexp_d, r = 0.5;\ - int low, high, mid;\ - int ubound, lbound;\ - char *p, pad = ' ';\ - int save_scale_factor, nb = 0;\ - try result;\ -\ - 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;\ - }\ -\ - rexp_d = calculate_exp_ ## x (-d);\ - if ((m > 0.0 && ((m < 0.1 - 0.1 * r * rexp_d) || (rexp_d * (m + r) >= 1.0)))\ - || ((m == 0.0) && !(compile_options.allow_std\ - & (GFC_STD_F2003 | GFC_STD_F2008))))\ - { \ - newf->format = FMT_E;\ - newf->u.real.w = w;\ - newf->u.real.d = d - comp_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)\ - { \ - volatile GFC_REAL_ ## x temp;\ - mid = (low + high) / 2;\ -\ - temp = (calculate_exp_ ## x (mid - 1) * (1 - r * rexp_d));\ -\ - 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;\ - }\ - }\ -\ - 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;\ -\ - finish:\ - result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \ - ndigits, edigits);\ - dtp->u.p.scale_factor = save_scale_factor;\ -\ - free (newf);\ -\ - if (nb > 0 && !dtp->u.p.g0_no_blanks)\ - {\ - p = write_block (dtp, nb);\ - if (p == NULL)\ - return;\ - if (result == FAILURE)\ - 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(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. */ /* Note: Before output_float is called, snprintf is used to print to buffer the @@ -1003,19 +933,15 @@ __qmath_(quadmath_snprintf) (buffer, sizeof buffer \ DTOA ## y\ \ - if (f->format != FMT_G)\ - output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \ - edigits);\ - else \ - output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \ - zero_flag, ndigits, edigits, comp_d);\ + output_float (dtp, f, buffer, size, 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, int comp_d) + int len) { #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18