The front-end parsing and translation are completed. We now need to determine how to implement UP, DOWN, COMPATIBLE, NEAREST, etc. in output_float.
> The front-end parsing and translation are completed. Except for the r* edit descriptors.
After studying the F2003 standard and our code in output_float, I believe what gfortran does now is ROUND="compatible". We round to the nearest and when there is a tie, we round away from zero. I think I can manage the other modes.
I agree with you for output: we currently do COMPATIBLE rounding and it shouldn't be too hard to change. The only case where we might not have control is for the last significant digit output by snprintf(), where we might rely somehow on the system library's implementation. For input, however, I think we don't really have control on what we do (we rely on the system strtod) and it is probably not that easy to change. I'll try to think about it some more. PS: for input, I have come to the conclusion that NEAREST and COMPATIBLE are really the same, as it is not possible for a decimal value to be exactly halfway between two machine-representable floating point values. Can someone confirm/infirm this?
> I agree with you for output: we currently do COMPATIBLE rounding I wanted to write this when I approved the code, but I found one *printf function and checking the POSIX standard one finds: "f, F The double argument shall be converted to decimal notation [...] The low-order digit shall be rounded in an implementation-defined manner." Therefore, I think COMPATIBLE is a good guess, but I'm not sure whether COMPATIBLE is always true. (In reply to comment #3) > I agree with you for output: we currently do COMPATIBLE rounding and it > shouldn't be too hard to change. The only case where we might not have control > is for the last significant digit output by snprintf(), where we might rely > somehow on the system library's implementation. Well, I would argue that this digit is the most important for rounding ... Does this mean we have to roll out our own implementation which replaces snprintf? > For input, however, I think we don't really have control on what we do > (we rely on the system strtod) and it is probably not that easy to change. I probably miss something, but where on input do you need to round? > PS: for input, I have come to the conclusion that NEAREST and COMPATIBLE are > really the same, as it is not possible for a decimal value to be exactly > halfway between two machine-representable floating point values. Can someone > confirm/infirm this? I wanted to mention BOZ I/O, but there you don't have a decimal point. But still I'm not sure that you cannot dream up a decimal number which has one digit too much for the machine-representible FP number. But in any case NEAREST is contained in COMPATIBLE as NEAREST does not describe what to to if the numbers are equidistant.
(In reply to comment #4) > Well, I would argue that this digit is the most important for rounding... Only for list-directed output and if people use really large formats (which is a lot of cases). Otherwise, not. > Does this mean we have to roll out our own implementation which replaces > snprintf? Might be. > I probably miss something, but where on input do you need to round? Rounding is done for all I/O, including input. For example, if you have: read ("0.1", *, rounding="up") x read ("0.1", *, rounding="down") y then you will have x == nearest(y, -1.)
Subject: Bug 35862 Author: jvdelisle Date: Tue Sep 29 02:44:38 2009 New Revision: 152262 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152262 Log: 2009-09-28 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR fortran/35862 * io.c (format_token): Add enumerators for rounding format specifiers. (format_lex): Tokenize the rounding format specifiers. (gfc_match_open): Enable rounding modes in OPEN statement. Modified: trunk/gcc/fortran/ChangeLog trunk/gcc/fortran/io.c
Subject: Bug 35862 Author: jvdelisle Date: Tue Sep 29 02:47:54 2009 New Revision: 152263 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152263 Log: 2009-09-28 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR libgfortran/35862 * io.h (gfc_unit): Add round_status. (format_token): Add enumerators for rounding format specifiers. * transfer.c (round_opt): New options table. (formatted_transfer_scalar_read): Add set round_status for each rounding format token. (formatted_transfer_scalar_write): Likewise. * format.c (format_lex): Tokenize the rounding format specifiers. (parse_format_list): Parse the rounding format specifiers. * write_float.def (outout_float): Modify rounding code to use new variable rchar to set the appropriate rounding. Fix some whitespace. * unit.c (get_internal_unit): Initialize rounding mode for internal units. (init_units): Likewise. Modified: trunk/libgfortran/ChangeLog trunk/libgfortran/io/format.c trunk/libgfortran/io/io.h trunk/libgfortran/io/transfer.c trunk/libgfortran/io/unit.c trunk/libgfortran/io/write_float.def
Subject: Bug 35862 Author: jvdelisle Date: Tue Sep 29 02:50:48 2009 New Revision: 152264 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152264 Log: 2009-09-28 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR libgfortran/35862 * gfortran.dg/round_1.f03: New test. * gfortran.dg/f2003_io_3.f03: Update test. Added: trunk/gcc/testsuite/gfortran.dg/round_1.f03 Modified: trunk/gcc/testsuite/ChangeLog trunk/gcc/testsuite/gfortran.dg/f2003_io_3.f03
Rounding modes are now implemented for formatted WRITE operations. I do not have a clear enough idea of what the rounding modes really mean for READ operations. For example: From Comment #4 "Rounding is done for all I/O, including input. For example, if you have: read ("0.1", *, rounding="up") x read ("0.1", *, rounding="down") y then you will have x == nearest(y, -1.)" Is this really what is meant by rounding up or down? The internal representation of .1 may already be rounded up or down from a binary representation that carried more bits in the conversion process. If we do not have access to more bits, how are we to know whether it is not already "half way in between" Similarly with nearest. By definition the result of the intrinsic nearest is already the next bit up or down from the internally stored value. So what does "rounding" ,mean here? Does it mean if the LSB is a 1 turn it into a 0 and if it is a 0, leave it a zero? This all seems absurd to my practical side. The Fortran standard seems very vague to me and yet it does imply rounding on input. Suggestions anyone?
(In reply to comment #9) > Rounding modes are now implemented for formatted WRITE operations. I do not > have a clear enough idea of what the rounding modes really mean for READ > operations. For example: > > "Rounding is done for all I/O, including input. For example, if you have: > read ("0.1", *, rounding="up") x > read ("0.1", *, rounding="down") y > then you will have x == nearest(y, -1.)" For REAL(8) I would expect for 0.1_8 with rounding up and down: 0.10000000000000001 9.99999999999999917E-002 or written as binary: 11111110111001100110011001100110011001100110011001100110011010 (= 0.1 rnd up) 11111110111001100110011001100110011001100110011001100110011001 (= 0.1 rnd down) > Is this really what is meant by rounding up or down? I think so. > The internal representation of .1 may already be rounded up or down > from a binary representation that carried more bits in the conversion > process. If we do not have access to more bits, how are we to know > whether it is not already "half way in between" Well, the rounding has to be based on the string (= decimal) number; if you already have a binary number, it is too late. > Similarly with nearest. Using "nearest" with "-1." in this example only works if you know that "read" automatically rounds up for 0.1. * * * Handling this efficiently and reliably is presumably nontrivial. See for instance: http://netlib.org/fp/ (esp. http://netlib.org/fp/gdtoa.tgz) http://dx.doi.org/10.1145/989393.989430 http://www.informatik.uni-trier.de/Reports/TR-08-2004/rnc6_10_hack.pdf and references therein.
There is probably a bug with round to nearest for values below 1: print '(RN, 4F10.3)', 0.0625, 0.1875 print '(RN, 4F10.2)', 0.125, 0.375, 1.125, 1.375 print '(RN, 4F10.1)', 0.25, 0.75, 1.25, 1.75 print '(RN, 4F10.0)', 0.5, 1.5, 2.5, 3.5 end gives 0.063 0.188 0.13 0.38 1.12 1.38 0.3 0.8 1.2 1.8 1. 2. 2. 4. ifort gives (what I was expecting): 0.062 0.188 0.12 0.38 1.12 1.38 0.2 0.8 1.2 1.8 0. 2. 2. 4.
Reply to comment #10: Thanks for the input and references. The perspective of rounding before conversion to binary representation is interesting. I will think on that a bit. Reply to comment #11: Thanks for test case. I will try to fix this.
There are two unrelated bugs in output rounding. The following patch fixes both. The smaller hunks were obvious. The larger gets rid of a hack I did to handle the special case with d=0. The rounding logic always saw a value of 1.0. Turns out with this special case, we have no leading digit from the sprintf routines so the rounding logic is looking in the wrong places. Fixed by shifting digits left and adjusting the exponent. Regression tested and passes the test case from Dominique in comment #11. Index: write_float.def =================================================================== --- write_float.def (revision 152422) +++ write_float.def (working copy) @@ -141,6 +141,14 @@ output_float (st_parameter_dt *dtp, const fnode *f switch (ft) { case FMT_F: + if (d == 0 && e <= 0 && dtp->u.p.scale_factor == 0) + { + for (i = ndigits - 1; i >= 0; i--) + digits[i] = digits[i - 1]; + digits[0] = '0'; + e++; + } + nbefore = e + dtp->u.p.scale_factor; if (nbefore < 0) { @@ -255,7 +263,7 @@ output_float (st_parameter_dt *dtp, const fnode *f case ROUND_NEAREST: /* Round compatible unless there is a tie. A tie is a 5 with all trailing zero's. */ - i = nafter + 1; + i = nafter + nbefore; if (digits[i] == '5') { for(i++ ; i < ndigits; i++) @@ -264,7 +272,7 @@ output_float (st_parameter_dt *dtp, const fnode *f goto do_rnd; } /* It is a tie so round to even. */ - switch (digits[nafter]) + switch (digits[nafter + nbefore - 1]) { case '1': case '3': @@ -818,14 +826,6 @@ sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" 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\
Correction to comment #13, Shifting bytes right not left. Also, I have completed testing using memmove instead of the explicit for loop.
Subject: Bug 35862 Author: jvdelisle Date: Tue Oct 6 03:08:20 2009 New Revision: 152483 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152483 Log: 2009-10-05 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR libgfortran/35862 * write_float.def (outout_float): Fix handling of special case where no digits after the decimal point and values less than 1.0. Adjust index into digits string. (WRITE_FLOAT): Remove special case code from macro. Modified: trunk/libgfortran/ChangeLog trunk/libgfortran/io/write_float.def
Subject: Bug 35862 Author: jvdelisle Date: Tue Oct 6 03:12:21 2009 New Revision: 152484 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152484 Log: 2009-10-05 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR libgfortran/35862 * gfortran.dg/round_2.f03: New test. Added: trunk/gcc/testsuite/gfortran.dg/round_2.f03 Modified: trunk/gcc/testsuite/ChangeLog
round_2.f03 fails on hppa-unknown-linux-gnu: Excess errors: /home/dave/gnu/gcc-4.5/gcc/gcc/testsuite/gfortran.dg/round_2.f03:6: Error: Missing kind-parameter at (1) /home/dave/gnu/gcc-4.5/gcc/gcc/testsuite/gfortran.dg/round_2.f03:9: Error: Missing kind-parameter at (1) /home/dave/gnu/gcc-4.5/gcc/gcc/testsuite/gfortran.dg/round_2.f03:12: Error: Missing kind-parameter at (1) /home/dave/gnu/gcc-4.5/gcc/gcc/testsuite/gfortran.dg/round_2.f03:15: Error: Missing kind-parameter at (1) UNRESOLVED: gfortran.dg/round_2.f03 -O0 compilation failed to produce executable
Subject: Re: [F2003] Implement new rounding modes for run time On Sat, 2009-10-10 at 15:45 +0000, danglin at gcc dot gnu dot org wrote: > > ------- Comment #17 from danglin at gcc dot gnu dot org 2009-10-10 15:45 ------- > round_2.f03 fails on hppa-unknown-linux-gnu: I was waiting for this to show up. We cab either XFAIL it or adjust the kind selection to use kind=8 if greater than 8 does not exist on the test system.
Reoly to Comment #17: See PR41612. I committed a fix to round_2.f03 that tests for the kind parameter being valid and also adds kind=8 checking.
Subject: Bug 35862 Author: jvdelisle Date: Sat Oct 10 18:57:35 2009 New Revision: 152627 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152627 Log: 2009-10-10 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR libgfortran/35862 * gfortran.dg/round_2.f03: Eliminate possible compile error. Modified: trunk/gcc/testsuite/ChangeLog trunk/gcc/testsuite/gfortran.dg/round_2.f03
Subject: Bug 35862 Author: jvdelisle Date: Sat Oct 10 23:02:11 2009 New Revision: 152632 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=152632 Log: 2009-10-10 Jerry DeLisle <jvdelisle@gcc.gnu.org> PR libgfortran/35862 * gfortran.dg/round_2.f03: Eliminate possible compile error. Use max function correctly. Modified: trunk/gcc/testsuite/ChangeLog trunk/gcc/testsuite/gfortran.dg/round_2.f03
Implemented: Rounding on output. TODO: Rounding on input. Current result for the program below is: 0.10000000149011611938E+00 0.10000000149011611938E+00 0.10000000000000000555E+00 0.10000000000000000555E+00 Expected: 0.10000000149011611938E+00 0.99999882280826568604E-01 0.10000000000000000555E+00 0.99999999999999783507E-01 (i.e. '"round up" == "round down" + epsilon(r)') (Side remark: None of my tested compiles handles this; not gfortran, not ifort nor crayftn.) Test program: character(len=5) :: str real(4) :: r4_1, r4_2 real(8) :: r8_1, r8_2 str = "0.1" read(str,'(ru,g3.1)') r4_1 read(str,'(rd,g3.1)') r4_2 read(str,'(ru,g3.1)') r8_1 read(str,'(rd,g3.1)') r8_2 print '(e30.20)', r4_1 print '(e30.20)', r4_2 print '(e30.20)', r8_1 print '(e30.20)', r8_2 end
(In reply to Tobias Burnus from comment #22) > Implemented: Rounding on output. It seems as if one handles this best via: #ifdef HAVE_FENV_H && (defined(FE_DOWNWARD) || defined(FE_TONEAREST) \ || defined(FE_TOWARDZERO) || defined(FE_UPWARD)) int mode = fegetround (); if (mode != wished_mode) fesetround(round_dir); // Ignore result value; should be == 0 #endif ... #ifdef HAVE_FENV_H && (defined(FE_DOWNWARD) || defined(FE_TONEAREST) \ || defined(FE_TOWARDZERO) || defined(FE_UPWARD)) if (mode != wished_mode) fesetround(round_dir); // Ignore result value #endif Or, possibly better, via some wrapper function in libgfortran/config/fpu-*.
Created attachment 30458 [details] Proof of concept (.c file) Proof of concept test case (in C, tested on x86-64-gnu linux with libc 2.17). [Somehow printf doesn't like my long double/__float128 example and prints 0.0 (long double) or garbage numbers (__float128).] Default 0.1: 0.1000000015 / 0.100000000000000006 -0.1: -0.1000000015 / -0.100000000000000006 FE_DOWNWARD 0.1: 0.0999999940 / 0.099999999999999991 -0.1: -0.1000000015 / -0.100000000000000006 FE_TONEAREST 0.1: 0.1000000015 / 0.100000000000000006 -0.1: -0.1000000015 / -0.100000000000000006 FE_TOWARDZERO 0.1: 0.0999999940 / 0.099999999999999991 -0.1: -0.0999999940 / -0.099999999999999991 FE_UPWARD 0.1: 0.1000000015 / 0.100000000000000006 -0.1: -0.0999999940 / -0.099999999999999991
Rounding on input. Patch, which sets the CPU rounding mode and relies on strtof to honour it: http://gcc.gnu.org/ml/fortran/2013-07/msg00042.html
(In reply to Tobias Burnus from comment #24) > [Somehow printf doesn't like my long double/__float128 example and prints > 0.0 (long double) or garbage numbers (__float128).] You should use %25.22Lf instead of %25.22lf.
Author: burnus Date: Sun Jul 21 11:54:27 2013 New Revision: 201093 URL: http://gcc.gnu.org/viewcvs?rev=201093&root=gcc&view=rev Log: 2013-07-21 Tobias Burnus <burnus@net-b.de> PR fortran/35862 * libgfortran.h (GFC_FPE_DOWNWARD, GFC_FPE_TONEAREST, GFC_FPE_TOWARDZERO, GFC_FPE_UPWARD): New defines. 2013-07-21 Tobias Burnus <burnus@net-b.de> PR fortran/35862 * libgfortran.h (set_fpu_rounding_mode, get_fpu_rounding_mode): New prototypes. * config/fpu-387.h (set_fpu_rounding_mode, get_fpu_rounding_mode): New functions. * config/fpu-aix.h (set_fpu_rounding_mode, get_fpu_rounding_mode): Ditto. * config/fpu-generic.h (set_fpu_rounding_mode, get_fpu_rounding_mode): Ditto. * config/fpu-glibc.h (set_fpu_rounding_mode, get_fpu_rounding_mode): Ditto. * config/fpu-sysv.h (set_fpu_rounding_mode, get_fpu_rounding_mode): Ditto. * configure.ac: Check for fp_rnd and fp_rnd_t. * io/io.h (enum unit_round): Use GFC_FPE_* for the value. * io/read.c (convert_real): Set FP ronding mode. * Makefile.in: Regenerate. * aclocal.m4: Regenerate. * config.h.in: Regenerate. * configure: Regenerate. 2013-07-21 Tobias Burnus <burnus@net-b.de> PR fortran/35862 * gfortran.dg/round_4.f90: New. Added: trunk/gcc/testsuite/gfortran.dg/round_4.f90 Modified: trunk/gcc/fortran/ChangeLog trunk/gcc/fortran/libgfortran.h trunk/gcc/testsuite/ChangeLog trunk/libgfortran/ChangeLog trunk/libgfortran/Makefile.in trunk/libgfortran/aclocal.m4 trunk/libgfortran/config.h.in trunk/libgfortran/config/fpu-387.h trunk/libgfortran/config/fpu-aix.h trunk/libgfortran/config/fpu-generic.h trunk/libgfortran/config/fpu-glibc.h trunk/libgfortran/config/fpu-sysv.h trunk/libgfortran/configure trunk/libgfortran/configure.ac trunk/libgfortran/io/io.h trunk/libgfortran/io/read.c trunk/libgfortran/libgfortran.h
Close as FIXED (on the trunk [= GCC 4.9]) - as now also the input rounding (READ) works. Note: - COMPATIBLE rounding is not supported for READ. (Compatible = like nearest, except it should always round away from zero for a tie.) - But I think this shouldn't matter much for input rounding. - Not on all systems, rounding on input works. It relies on rounding support of strtof/strtod/strtold. At least GLIBC handles it. (Setting the rounding mode of the CPU works at least for GLIBC, SysV, AIX and on x86/x86-64 systems.)