User account creation filtered due to spam.

Bug 35862 - [F2003] Implement new rounding modes for run time
Summary: [F2003] Implement new rounding modes for run time
Status: RESOLVED FIXED
Alias: None
Product: gcc
Classification: Unclassified
Component: libfortran (show other bugs)
Version: 4.4.0
: P3 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2008-04-07 22:33 UTC by Jerry DeLisle
Modified: 2013-07-21 12:00 UTC (History)
3 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed: 2008-04-09 05:15:10


Attachments
Proof of concept (.c file) (355 bytes, text/plain)
2013-07-04 16:03 UTC, Tobias Burnus
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Jerry DeLisle 2008-04-07 22:33:41 UTC
The front-end parsing and translation are completed.  We now need to determine how to implement UP, DOWN, COMPATIBLE, NEAREST, etc. in output_float.
Comment 1 Tobias Burnus 2008-04-08 08:49:57 UTC
> The front-end parsing and translation are completed.
Except for the r* edit descriptors.
Comment 2 Jerry DeLisle 2008-04-09 05:15:10 UTC
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.
Comment 3 Francois-Xavier Coudert 2008-04-15 14:02:38 UTC
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?
Comment 4 Tobias Burnus 2008-04-15 19:08:20 UTC
> 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.
Comment 5 Francois-Xavier Coudert 2008-04-15 20:59:24 UTC
(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.)
Comment 6 Jerry DeLisle 2009-09-29 02:44:56 UTC
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

Comment 7 Jerry DeLisle 2009-09-29 02:48:07 UTC
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

Comment 8 Jerry DeLisle 2009-09-29 02:51:01 UTC
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

Comment 9 Jerry DeLisle 2009-10-01 01:37:39 UTC
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?
Comment 10 Tobias Burnus 2009-10-01 07:35:37 UTC
(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.
Comment 11 Dominique d'Humieres 2009-10-01 20:18:28 UTC
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.
Comment 12 Jerry DeLisle 2009-10-02 02:32:03 UTC
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.
Comment 13 Jerry DeLisle 2009-10-03 15:10:56 UTC
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\
 

Comment 14 Jerry DeLisle 2009-10-03 16:17:22 UTC
Correction to comment #13,  Shifting bytes right not left.  Also, I have completed testing using memmove instead of the explicit for loop.
Comment 15 Jerry DeLisle 2009-10-06 03:08:33 UTC
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

Comment 16 Jerry DeLisle 2009-10-06 03:12:32 UTC
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

Comment 17 John David Anglin 2009-10-10 15:45:32 UTC
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
Comment 18 Jerry DeLisle 2009-10-10 16:03:06 UTC
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.

Comment 19 Jerry DeLisle 2009-10-10 17:55:24 UTC
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. 
Comment 20 Jerry DeLisle 2009-10-10 18:57:47 UTC
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

Comment 21 Jerry DeLisle 2009-10-10 23:02:26 UTC
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

Comment 22 Tobias Burnus 2010-07-12 08:07:00 UTC
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
Comment 23 Tobias Burnus 2013-07-04 10:08:47 UTC
(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-*.
Comment 24 Tobias Burnus 2013-07-04 16:03:58 UTC
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
Comment 25 Tobias Burnus 2013-07-17 09:04:14 UTC
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
Comment 26 Uroš Bizjak 2013-07-17 09:37:44 UTC
(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.
Comment 27 Tobias Burnus 2013-07-21 11:54:54 UTC
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
Comment 28 Tobias Burnus 2013-07-21 12:00:14 UTC
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.)