GCC Bugzilla – Bug 48906

Wrong rounding results with -m32

Last modified: 2014-01-11 21:44:02 UTC

This PR is to carry forward the remaining issue identified in PR48602 comment #37 Here is a test case showing the problem. print "(ru,g15.2)", .099d0 ! > 0.10E+00 expected 0.10 print "(rc,g15.1)", .095d0 ! > 0.1E+00 expected 0.1 print "(rc,g15.2)", .0995d0 ! > 0.10E+00 expected 0.10 print "(ru,g15.3)", .0999d0 ! > 0.100E+00 expected 0.100 end $ gfc bug.f03 $ ./a.out 0.10 0.1 0.10 0.100 $ gfc -O1 bug.f03 $ ./a.out 0.10 0.1 0.10 0.100 $ gfc -m32 -O1 bug.f03 $ ./a.out 0.10E+00 0.1E+00 0.10E+00 0.100E+00 $ gfc -m32 bug.f03 $ ./a.out 0.10E+00 0.1E+00 0.10E+00 0.100E+00 The problem is rooted in excess precision and that we use sensitive floating point operations to determine the formatting parameters. I am seeking a solution that will avoid using floating point to do this.

Repeating Thomas T. suggestion: "We can treat FMT_G like FMT_E. After the rounding step, i.e. when the final value of the variable e in output_float() is known (at the label 'skip'), if e is within limits for F editing, revert the formatting to F and add the blanks at the end." This is my plan A and I have a patch started. Another idea would be to perform the floating point operations and tests with mpfr. This will avoid the -m32 issue. As a side note, mpfr is compiled with -ffloat-store on 32bit platforms by default. I will take this idea as plan B or C. ;)

(In reply to comment #1) > Another idea would be to perform the floating point operations and tests with > mpfr. This will avoid the -m32 issue. As a side note, mpfr is compiled with > -ffloat-store on 32bit platforms by default. I will take this idea as plan B > or C. ;) If the "floating point operations and tests" are part of the runtime library, then you cannot currently use mpfr because the final executable is not linked with gmp and mpfr.

Status: I have managed to eliminate the OUTPUT_FLOAT_G macro completely and all floating point operations. What remains is to adjust the trailing blanks. The patch is simple looking once you get there, though a bit subtle to sort out initially, especially when I get only 20 or 30 minute intervals to work on this every two or three days. ;) Stay tuned.

Way to go! I'll be happy to test.

Thomas, What are the correct results for this test case: implicit none integer, parameter :: RT = 8 write(*, "(rc,f11.2,4x,'<RC')") .995_RT write(*, "(rc,g15.2,'<RC')") .995_RT write(*, "(rd,f11.2,4x,'<RD')") .995_RT write(*, "(rd,g15.2,'<RD')") .995_RT write(*, "(ru,f11.2,4x,'<RU')") .995_RT write(*, "(ru,g15.2,'<RU')") .995_RT write(*,'(a)') "12345678901234567890" end I have finished re-factoring the code, but have some remaining breakages to fix. I just want to confirm what the above should give vs what I am getting at the moment. I will post a patch to test after I get all known breakage fixed. Thanks.

(In reply to comment #5) As a general rule, d specifies the number of significant digits in the result, i.e. the number of digits counting from the first non-zero digit. So in the example all numbers should be formatted with 2 significant digits, i.e. 2 decimal digits if the value rounded to 2 significant digits is smaller than 1, and 1 decimal digit if the value rounded to 2 significant digits is larger than or equal to 1. > implicit none > integer, parameter :: RT = 8 > write(*, "(rc,f11.2,4x,'<RC')") .995_RT > write(*, "(rc,g15.2,'<RC')") .995_RT > write(*, "(rd,f11.2,4x,'<RD')") .995_RT > write(*, "(rd,g15.2,'<RD')") .995_RT > write(*, "(ru,f11.2,4x,'<RU')") .995_RT > write(*, "(ru,g15.2,'<RU')") .995_RT > write(*,'(a)') "12345678901234567890" > end 0.99 <RC 0.99 <RC 0.99 <RD 0.99 <RD 1.00 <RU 1.0 <RU 12345678901234567890 and write(*, "(ru,f11.1,4x,'<RU')") .995_RT gets 1.0 <RU

It was the last case I was having trouble with. I also was getting different results from different compilers, so I needed to confirm independently. I have this part working now. Still have some bugs here.

Status: I am down to about 5 testsuite failures on the patch for this. There is a lot if interplay going on, so i will be factoring the code some as part of the cleanup. Stay tuned.

Created attachment 24364 [details] Preliminary patch I am getting two testsuite regressions with this patch. At least one of those I think is an adjustment needed to the test case. Possibly the same with the other, but I am still testing. After fixing the above, all that remains is a little code cleanup. Please begin testing this patch and report back.

The test cases from above still fail (but with a different result): print "(ru,g15.2)", .099d0 ! 1.0E-01 expected 0.10 print "(rc,g15.1)", .095d0 ! 1.E-01 expected 0.1 print "(rc,g15.2)", .0995d0 ! 1.0E-01 expected 0.10 print "(ru,g15.3)", .0999d0 ! 1.00E-01 expected 0.100 Also, G0.d tests fail in the E formatting case, the chosen format is E25.(d)E3 and not the minimal width required for G0.d: print "(rc,g0.4,'<')", .1234d-3 ! fixed width, not G0.d format print "(rc,g0.4,'<')", .1234d-2 ! fixed width, not G0.d format print "(rc,g0.4,'<')", .1234d-1 ! fixed width, not G0.d format print "(rc,g0.4,'<')", .1234d0 print "(rc,g0.4,'<')", .1234d1 print "(rc,g0.4,'<')", .1234d2 print "(rc,g0.4,'<')", .1234d3 print "(rc,g0.4,'<')", .1234d4 print "(rc,g0.4,'<')", .1234d5 ! fixed width, not G0.d format print "(rc,g0.4,'<')", .1234d6 ! fixed width, not G0.d format print "(rc,g0.4,'<')", .1234d7 ! fixed width, not G0.d format

The scale factor is applied to F editing, but it shouldn't. See Fortran 2008: NOTE 10.20 The scale factor has no effect on output unless the magnitude of the datum to be edited is outside the range that permits effective use of F editing. print "(rc,1p,g15.4e2,'<')", 0.1234 ! 1.2340E-01 expected 0.1234 print "(rc,g15.4e2,'<')", 0.1234 ! 0.1234 Additionally, there is a strange digit mixup in this case: print "(rc,-1p,g15.4e2,'<')", 0.1234 ! 1.0234 expected 0.1234

The following examples fail: print "(rc,g10.2,'<')", 99.5 ! 10. expected 0.10E+03 print "(rc,g10.2,'<')", 995. ! 1.0E+03 expected 0.10E+04 print "(rc,g10.3,'<')", 999.5 ! 100. expected 0.100E+04 print "(rc,g10.3,'<')", 9995. ! 1.0E+04. expected 0.100E+05

OK, thanks for test cases, these are very helpful.

Created attachment 24399 [details] Updated preliminary patch This updated patch resolves the issues identified so far. Please continue testing. I have regression tested and the testsuite case gfortran.dg/pr20755.f fails because with the patch I am ignoring the scale factor. Perhaps I have been too aggressive with that. Regardless, please continue testing while I study the scale factor question a bit more.

I hadn't really thought of that, but now I see what a PITA a scale factor > 0 is going to be (scale factor < 0 is simply padding with zeros and shifting digits): We have to ignore the scale factor, until we decide that we use E editing, at which point we need an additional digit in the result. For efficiency: Maybe we can "decide that we use E editing" without actually overwriting the original digits at that point yet, so we can go back when we actually switch to E editing and re-do the rounding with the extra digit.

print "(rc,g11.1)", 2. ! 0.2E+01 expected 2. print "(rc,g11.2)", 20. ! 0.20E+02 expected 20. print "(rc,g11.3)", 200. ! 0.200E+03 expected 200. print "(rc,g11.4)", 2000. ! 0.2000E+04 expected 2000. print "(rc,g11.2)", .04 ! 0.40 expected 0.40E-01 print "(rc,g11.3)", .04 ! 0.400 expected 0.400E-01

As you can tell, finding the magic formula with this new approach is not a trivial effort. ;) My fear is that there is an unlimited number of corner cases.

Created attachment 24406 [details] New updated patch This updated patch takes care of Comment #16. Also, the test case, pr20755.f was originally intended to pass only with -std=legacy. This is to mimic g77 which does not ignore the scale factor. At least recently, gfortran passes this test case regardless of the -std=. I am tempted to delete the test case, since it looks like it goes against current standards. Any opinions about this?

> Also, the test case, pr20755.f was originally intended to pass only with > -std=legacy. This is to mimic g77 which does not ignore the scale factor. At > least recently, gfortran passes this test case regardless of the -std=. I am > tempted to delete the test case, since it looks like it goes against current > standards. Any opinions about this? Now you've lost me. From what I understand PR20755 was about G editing resetting the scale factor in certain cases (to 0 in the example) and not restoring it afterwards (to 1 in this case). This is a bug both in legacy and in Fortran standard modes. The bug was fixed and a testcase was created. The result is standard-compliant. Why would you want to remove the testcase? The point of having a testcase is to prevent a reappearing of the same bug in the future, so maybe we should keep it.

(In reply to comment #18) > Created attachment 24406 [details] > New updated patch > > This updated patch takes care of Comment #16. Unfortunately, now the other testcases fail again (although with still another results): print "(ru,g15.2)", .099d0 ! > 0.10E-01 expected 0.10 print "(rc,g15.1)", .095d0 ! > 0.1E-01 expected 0.1 print "(rc,g15.2)", .0995d0 ! > 0.10E-01 expected 0.10 print "(ru,g15.3)", .0999d0 ! > 0.100E-01 expected 0.100 print "(rc,g10.2)", 99.5 ! 0.10E+02 expected 0.10E+03 print "(rc,g10.2)", 995. ! 0.10E+03 expected 0.10E+04 print "(rc,g10.3)", 999.5 ! 0.100E+03 expected 0.100E+04 print "(rc,g10.3)", 9995. ! 0.100E+04 expected 0.100E+05 It seems there is still a bug that doesn't increment the exponent in the "overflow" case, i.e. when going from 0.9999...E(x) to 0.1000...E(x+1) The exponent test (decision between F and E editing) must be done after that! I can not find this in the code right now.

Regarding comment #18, the test case passes regardless 0f the -std=legacy option specified at the beginning of the source file. I am questioning why that option is ignored in current gfortran releases. Somewhere along the way, something changed. I am also asking, do we want the test case to pass under -std=legacy only, or pass all the time. Regarding, not incrementing the exponent, I missed that and will update on next iteration. thanks for reviews.

Created attachment 24414 [details] Fixed revised patch This patch fixes my omissions on the last patch. In addition, studying the test case pr20755.f, I see that it has the wrong number of significant digits in the expected result. I have revised it now so that it passes. I have also added a new test case, to capture the various errors we found here while developing the patch Feel free to continue testing and if no new errors show up, I will submit to the gfortran list for approval. /me crosses his fingers. ;)

Patch submitted to list for approval.

(In reply to comment #23) > Patch submitted to list for approval. For a scale factor 0, we are done. Good work, thank you! A scale factor != 0 does not work yet, you wrote you are still working on it, is that correct? print "(-2pg12.3)", 0.02 ! 0.200E-01 expected 0.002E+01 print "(-1pg12.3)", 0.02 ! 0.200E-01 expected 0.020E+00 print "(0pg12.3)", 0.02 ! 0.200E-01 print "(1pg12.3)", 0.02 ! 0.200E-01 expected 2.000E-02 print "(2pg12.3)", 0.02 ! 0.200E-01 expected 20.00E-03 --- gcc/testsuite/gfortran.dg/pr20755.f (revision 174320) +++ gcc/testsuite/gfortran.dg/pr20755.f (working copy) @@ -5,8 +5,8 @@ character*30 s write (s,2000) 0.0, 0.02 - if (s .ne. " 0.00 2.000E-02") call abort + if (s .ne. " 0.00 0.200E-01") call abort write (s,2000) 0.01, 0.02 - if (s .ne. " 1.000E-02 2.000E-02") call abort + if (s .ne. " 0.100E-01 0.200E-01") call abort 2000 format (1PG12.3,G12.3) end I don't agree with the changes to this testcase, since the scale factor does not work yet, and the test was correct before. --- gcc/testsuite/gfortran.dg/fmt_g0_6.f08 (revision 174320) +++ gcc/testsuite/gfortran.dg/fmt_g0_6.f08 (working copy) @@ -57,7 +57,7 @@ contains do dec = d, 0, -1 lower = 10.0_RT ** (d - 1 - dec) - r * 10.0_RT ** (- dec - 1) upper = 10.0_RT ** (d - dec) - r * 10.0_RT ** (- dec) - if (lower <= mag .and. mag < upper) then + if (lower < mag .and. mag <= upper) then write(fmt_f, "('R', a, ',F', i0, '.', i0, ',', i0, 'X')") roundmode, w - n, dec, n exit end if I don't agree with this change, since it was according to the Fortran standard before (it does not actually make a difference for the values we currently check in this test, though.)

On 06/06/2011 01:38 AM, thenlich at users dot sourceforge.net wrote: > For a scale factor 0, we are done. Good work, thank you! > > A scale factor != 0 does not work yet, you wrote you are still working on it, > is that correct? I am now. ;) > > print "(-2pg12.3)", 0.02 ! 0.200E-01 expected 0.002E+01 > print "(-1pg12.3)", 0.02 ! 0.200E-01 expected 0.020E+00 > print "(0pg12.3)", 0.02 ! 0.200E-01 > print "(1pg12.3)", 0.02 ! 0.200E-01 expected 2.000E-02 > print "(2pg12.3)", 0.02 ! 0.200E-01 expected 20.00E-03 > My confusion seems to be when scale factor is to be ignored and when not, I will give the standard another read.

(In reply to comment #25) > My confusion seems to be when scale factor is to be ignored and when not, I > will give the standard another read. As it happens, you're not the only one to find (this part of) the Fortran standard confusing and buggy, so if all goes as planned there will be a rewrite: --------------------------------------------------------------------- 11-xxx To: J3 From: John Reid and Thomas Henlich Subject: Interp: G editing for reals Date: 2011 June 5 --------------------------------------------------------------------- NUMBER: F08/xxxx TITLE: G editing for reals KEYWORDS: format, G editing DEFECT TYPE: Erratum STATUS: J3 consideration in progress QUESTION: 1. Gw.d editing for a real value that is in the range (0.1,10**d) and is not near an integer power of 10 uses F editing to produce exactly the same value as 0PEw.d editing would produce. For values in this range that are near an integer power of 10, is it intended that F editing be used to produce exactly the same value as 0PEw.d editing would produce? The rules in 10.7.5.2.2 usually have this effect, but the following examples illustrate exceptions for rounding UP and to ZERO. i. print "(ru,g11.2)", -9.95 or print "(rz,g11.2)", -9.95 Here, 0PE11.2 editing would produce the value -0.99E+01, which can be represented as -9.9. However, the standard requires F7.0 to be used, which gives the result -9. Note that the standard requires print "(rd,g11.2)", 9.95 and print "(rz,g11.2)", 9.95 to give the result 9.9. ii. print "(ru,0p,g11.2)", -99.5 The standard requires 0PE11.2 editing to be used, which gives -0.99E+02. This is representable as -99. iii. print "(ru,0p,g11.2)", 99. The standard requires 0PE11.2 editing to be used, which gives 0.99E+02. This is representable as 99. 2. COMPATIBLE and NEAREST modes of rounding differ only when the two nearest representable values are equidistant from the given value. The similarity appears not to be represented in the second table. What is meant by "if the higher value is even"? 3. Why is no account taken of the effects when PROCESSOR_DEFINED rounding is in effect? ANSWER: 1. Yes, this was the intention and it would be clearer for the standard to state this directly. It would also be easier for implementers to implement. 2. If the standard is rewritten as proposed in the first question, these further problems would be resolved. 3. If the standard is rewritten as proposed in the first question, PROCESSOR_DEFINED rounding would be covered. EDITS to 10-007r1: [258:14-20] In 10.7.5.2.2, replace paragraph 4 by "Otherwise, the method of representation in the output field depends on the internal value being edited. Let N be the value that would be output with 0PEw.dEe editing and let s be its exponent part unless N is identically 0 in which case let s be 1. Let k be the scale factor (10.8.5). Let b be a blank. For Gw.d editing, if s lies in the range 0 <= s <= d, F(w-4).(d-s),4('b') editing is used to represent N in the output field; otherwise, kPEw.d editing is used to represent the internal value. For Gw.dEe editing, if s lies in the range 0 <= s <= d, F(w-n).(d-s),n('b') editing where n=e+2 is used to represent N in the output field; otherwise, kPEw.dEe editing is used to represent the internal value." SUBMITTED BY: John Reid and Thomas Henlich HISTORY: 11-xxx m195 Submitted ------------------------------------------------------------------------

> > print "(-2pg12.3)", 0.02 ! 0.200E-01 expected 0.002E+01 > print "(-1pg12.3)", 0.02 ! 0.200E-01 expected 0.020E+00 > print "(0pg12.3)", 0.02 ! 0.200E-01 > print "(1pg12.3)", 0.02 ! 0.200E-01 expected 2.000E-02 < Too many significant digits? > print "(2pg12.3)", 0.02 ! 0.200E-01 expected 20.00E-03 < Too many significant digits? Should these last two cases by 2.00E-02 and 20.0E-03 ? Otherwise we seem to be adding an extra significant digit. Help me understand this. Jerry

I had a look at the code and what I think we should do is the following: If G editing and a scale factor p != 0 is in effect, split the rounding step into 2 steps: First, check if overflow occurs (.9999.. => 1.0000...) but do not actually overwrite the digits yet. With this information, we determine the final value of e. If it turns out we need F editing: Repeat the rounding as above (or remember if overflow occured the first time), and this time actually overwrite the digits. If it turns out we need E editing: Set the new value for the number of digits according to p and repeat the rounding. It may help to refactor some of the code into a separate function, since we need to call it twice.

(In reply to comment #27) > > > > print "(-2pg12.3)", 0.02 ! 0.200E-01 expected 0.002E+01 > > print "(-1pg12.3)", 0.02 ! 0.200E-01 expected 0.020E+00 > > print "(0pg12.3)", 0.02 ! 0.200E-01 > > print "(1pg12.3)", 0.02 ! 0.200E-01 expected 2.000E-02 < Too many significant digits? > > print "(2pg12.3)", 0.02 ! 0.200E-01 expected 20.00E-03 < Too many significant digits? > > Should these last two cases by 2.00E-02 and 20.0E-03 ? Otherwise we seem to be > adding an extra significant digit. > > Help me understand this. Don't worry, it's not your fault... Say goodbye to the world of common sense, welcome to the world of the Fortran standard! ;-) "The scale factor k controls the decimal normalization (10.3.2, 10.8.5). If −d < k <= 0, the output field contains exactly |k| leading zeros and d − |k| significant digits after the decimal symbol. If 0 < k < d + 2, the output field contains exactly k significant digits to the left of the decimal symbol and d − k + 1 significant digits to the right of the decimal symbol. Other values of k are not permitted."

Created attachment 24454 [details] Patch for scale factor. PUBLIC DOMAIN

(In reply to comment #30) > Created attachment 24454 [details] > Patch for scale factor. PUBLIC DOMAIN A had a go at this. Feel free to improve.

I managed to get pr20755.f fixed last night with only a minor change and I found a problem elsewhere for fmt_g0_6.f08 which I am still working on. I will have a look at your patch in the next few days.

The last test case I am working is fmt_g0_6.f08. The apparent failing case is: print "(rc,g15.2)", 0.995000_8 Which is resulting in 0.99 and we expect it to be 1.0. However, the raw internal representation of 0.995 is not exact and is: 99499999999999999555 This places the value less then the threshold given in the table in the Standard, giving the correct F format as: print "(f11.2,4x)", 0.995000_8 ! 0.99 Instead of print "(f11.1,4x)", 0.995000_8 ! 1.0 Our F formatting in the test case is using the f11.1 and the logic does not look at the third digit and dutifully rounds the value as it should. The new G formatting is using exact integer logic and dutifully does look at the third digit. If we do: print "(rc,g15.2)", 0.9950001_8 The internal value does cross the 0.995 threshold and we do get the 1.0 results. I think this issues is one of those splitting hairs things. We could artificially round the digits string early before we process it, but this would effectively be rounding twice to get there. I would prefer to revise the test case like so: call check_all(0.995000000001_RT, 15, 2, 0) and be done with this.

Additional note: The standard states: "Let N be the magnitude of the internal value" The internal value is to be used to determine the conversion to F formatting. I think this adds to my point. I wonder if the standards committee knew the thresholds selected do not have an exact binary representation so that the internal values will always be above or below it.

(In reply to comment #33) > The last test case I am working is fmt_g0_6.f08. > > The apparent failing case is: > > print "(rc,g15.2)", 0.995000_8 > > Which is resulting in 0.99 and we expect it to be 1.0. > > However, the raw internal representation of 0.995 is not exact and is: > > 99499999999999999555 > > This places the value less then the threshold given in the table in the > Standard, giving the correct F format as: > > print "(f11.2,4x)", 0.995000_8 ! 0.99 > > Instead of > > print "(f11.1,4x)", 0.995000_8 ! 1.0 > > > Our F formatting in the test case is using the f11.1 and the logic does not > look at the third digit and dutifully rounds the value as it should. The new G > formatting is using exact integer logic and dutifully does look at the third > digit. > > If we do: > > print "(rc,g15.2)", 0.9950001_8 > > The internal value does cross the 0.995 threshold and we do get the 1.0 > results. > > I think this issues is one of those splitting hairs things. We could > artificially round the digits string early before we process it, but this would > effectively be rounding twice to get there. I would prefer to revise the test > case like so: > > call check_all(0.995000000001_RT, 15, 2, 0) > > and be done with this. The result is correct, but our expectation was wrong. 0.995000_8 which is internally 0.99499999999999999555... is smaller than 0.995 (exact value) so it must output 0.99, not 1.0. Please don't waste time with this test case (fmt_g0_6.f08), it is broken. I will fix this some time.

(In reply to comment #34) > Additional note: The standard states: > > "Let N be the magnitude of the internal value" > > The internal value is to be used to determine the conversion to F formatting. I > think this adds to my point. > > I wonder if the standards committee knew the thresholds selected do not have an > exact binary representation so that the internal values will always be above or > below it. Not always (e.g. 99.5 is the exact internal value), but it many cases, yes. This and other problems with G editing will be fixed in the standard soon, if all goes well. I helped to prepare a request for this, which has now been submitted: http://j3-fortran.org/doc/year/2011/11-174.txt

Updated patch posted for approval: http://gcc.gnu.org/ml/fortran/2011-06/msg00097.html Thomas, thanks for working the Standard issues and your testing.

The Fortran standards committee has voted to edit the standard: http://j3-fortran.org/doc/meeting/195/11-174r2.txt This makes our current approach standard compliant (with the corrigendum to be published) This also means an additional todo item for this bug is "If d is zero, kPEw.0 or kPEw.0Ee editing is used for Gw.0 editing or Gw.0Ee editing respectively."

Very Good, I will work toward this end.

This PR seems to have been fixed from 4.7.3 up to trunk. Could it be closed as fixed or is there some issue I am missing?

Closing this long PR and opening a new one for the cleanup mentioned in comment #38