[Bug libfortran/47945] New: REAL(8) output conversion error on MinGW32

[Bug libfortran/47945] New: REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

           Summary: REAL(8) output conversion error on MinGW32
           Product: gcc
           Version: 4.6.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libfortran
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: thenlich@users.sourceforge.net


The last digits of REAL(8) numbers are converted inaccurately on output.

For instance, the output of:

real(8) :: r8
real(16) :: r16

r8 = .14285714285714286d0
r16 = r8
write(*, '(f35.32)') r8
write(*, '(f35.32)') r16
end

is:
 0.14285714285714284921875000000000
 0.14285714285714284921269268124888

but the expected output is:
 0.14285714285714284921269268124888
 0.14285714285714284921269268124888

because the internal 64-bit representation of .14285714285714286 is
approximately equal to 0.14285714285714284921269268124888...

Arguably, an output of 0.14285714285714285000000000000000 would equally be
acceptable, because it is the shortest decimal representation of the same
internal value. (This is in fact what the GCC C Compiler on mingw32 does)

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #1 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-01 14:01:32 UTC ---
Created attachment 23501
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23501
Test case

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #2 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-01 14:01:56 UTC ---
Created attachment 23502
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23502
C test case

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[schwab@linux-m68k.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #3 from Andreas Schwab <schwab@linux-m68k.org> 2011-03-01 14:27:54 UTC ---
0.14285714285714284921875000000000 is still within the accuracy of IEEE double.
 All numbers map to the same IEEE double.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #4 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-01 14:44:54 UTC ---
(In reply to comment #3)
> 0.14285714285714284921875000000000 is still within the accuracy of IEEE double.
>  All numbers map to the same IEEE double.

This is technically correct; however a program should also observe portability
and usability rules and not just pick any number from a valid range, but the
number which is the most useful to the program user: either the one most common
across platforms and other compilers (0.14285714285714284921269268124888) or
the "shortest" one (0.14285714285714285000000000000000).

It is very hard to compare long numeric result listings obtained on different
platforms/compilers if some digits (however insignificant) never match.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[jvdelisle at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

Jerry DeLisle <jvdelisle at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jvdelisle at gcc dot
                   |                            |gnu.org

--- Comment #5 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> 2011-03-02 00:22:41 UTC ---
For real 8, it is meaningless to carry the conversion on adinfinatum just
because someone asked for extra width to be printed.  So the zeros shown are
perfectly OK.

On Cygwin I get with real 8:

0.14285714285714284921269000000000

and with nearest(r8, -1.0)

0.14285714285714282145712000000000

and nearest(r8, +1.0)

0.14285714285714287696827000000000

I would be curious about the last three non-zero digits you see on MingW

0.14285714285714284921875000000000

but I see no justification for the digits to match at all after the
...571428 sequence with two different implementations.  You just can not get
more precision then is really there, a finite number of bits. (IMHO)

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

Thomas Henlich <thenlich at users dot sourceforge.net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Severity|normal                      |minor

--- Comment #6 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-02 06:45:01 UTC ---
I still think it makes sense to demand that the same binary value always
converts to the same decimal number, regardless of compiler vendor, or
platform.

It should even be the same for the same internal value, when the variables are
of different real kinds (which now it isn't, on mingw). Consider:

real(8) :: r8
real(16) :: r16

r8 = .14285714285714286d0
r16 = r8
write(*, '(f35.32)') r8
write(*, '(f35.32)') r16
end

output:
 0.14285714285714284921875000000000
 0.14285714285714284921269268124888

Sometimes it is necessary to specify more decimal digits than necessary in the
edit descriptor because the magnitude of the result may vary. Now, if the same
binary value would always result in the same decimal number, that would make it
easier to compare results obtained from the same program, obtained with two
different compilers or on different platforms (a method used for program
verification).

The Fortran 2008 standard even demands this behaviour (in my interpretation):

===
10.7.2.3.7 I/O rounding mode

2. In what follows, the term "decimal value" means the exact decimal number as
given by the character string, while the term "internal value" means the number
actually stored in the processor.  For example, in dealing with the decimal
constant 0.1, the decimal value is the mathematical quantity 1/10, which has no
exact representation in binary form.  Formatted output of real data involves
conversion from an internal value to a decimal value; formatted input involves
conversion from a decimal value to an internal value.

3.  When the I/O rounding mode is UP, the value resulting from conversion shall
be the smallest representable value that is greater than or equal to the
original value. When the I/O rounding mode is DOWN, the value resulting
from conversion shall be the largest representable value that is less than or
equal to the original value. [etc]
===

In the example, 0.14285714285714284921269268124888 is the largest representable
(with 32 decimal digits) value that is less than the original value (binary
1.001001001001001001001001001001001001001001001001001 * 2^-3 = decimal
0.1428571428571428492126926812488818...).

The point is, what the standard demands is the closest approximation that can
be represented with the specified output precision.

The standard does not state that an implementation may truncate the result
after an arbitrary number of decimal digits (even if that number is higher than
the possible precision of the internal bit width used).

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[kargl at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |kargl at gcc dot gnu.org

--- Comment #7 from kargl at gcc dot gnu.org 2011-03-02 07:13:40 UTC ---
(In reply to comment #6)
> I still think it makes sense to demand that the same binary value always
> converts to the same decimal number, regardless of compiler vendor, or
> platform.

It does if one follows IEEE 754.  A number with 53-bits
of precision will give at most 17 decimal digits.


> It should even be the same for the same internal value, when the variables are
> of different real kinds (which now it isn't, on mingw). Consider:
> 
> real(8) :: r8
> real(16) :: r16
> 
> r8 = .14285714285714286d0
> r16 = r8
> write(*, '(f35.32)') r8
> write(*, '(f35.32)') r16
> end
> 
> output:
>  0.14285571471428284921875000000000
>  0.14285714285714284921269268124888

r8 has 53 bit of precision.
r16 has 113 bits of precision. 

r8 has 15 decimal digits.  r16 has 35
decimal digits.  You're getting well
beyond 15 digits in your r8 output.

> Sometimes it is necessary to specify more decimal digits than necessary in the
> edit descriptor because the magnitude of the result may vary.

Ah no!  With IEEE binary64 (ie double precision), you can expect
at most 15 decimal digits.  Nothing more and quite often quite
less.

laptop:kargl[203] cat > a.f90
program a
  print *, precision(1.d0), digits(1.d0)
end program a
laptop:kargl[204] 
laptop:kargl[204] gfc4x -o z a.f90
laptop:kargl[205] ./z
          15          53

If you expect more than 15 digits, then I think it may 
be profitable to re-read Goldberg's paper.

> 
> The Fortran 2008 standard even demands this behaviour (in my interpretation):
>
> ===
> 10.7.2.3.7 I/O rounding mode


gfortran is not a F2008 compiler and it currently does not support
I/O rounding modes.  I'll note that your program isn't setting a
rounding mode, so what happens in rounding UP and DOWN is irrelevant.
gfortran by default will use round to nearest.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

Tobias Burnus <burnus at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |burnus at gcc dot gnu.org

--- Comment #8 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-03-02 07:24:57 UTC ---
It is not completely clear to me what you get for your C test case. My feeling
is that the differences are due to the use of different "printf" routines.

For "double" or "REAL(8)", I think Microsoft's snprintf is used - while for
__float128 or REAL(16) a GLIBC-based snprintf is used in libquadmath.

Not surprisingly, I get on Linux the r16 result - which makes sense as all
snprintf routines are based on GLIBC (program from comment 6)
 0.14285714285714284921269268124888
 0.14285714285714284921269268124888

As both numbers are IEEE identical, one cannot really say that Microsoft libc's
implementation or libquadmath's is wrong.

Note: You could use libquadmath with __float128 also from C, cf.
http://gcc.gnu.org/onlinedocs/libquadmath/quadmath_005fsnprintf.html


(In reply to comment #7)
> > 10.7.2.3.7 I/O rounding mode
> 
> gfortran is not a F2008 compiler and it currently does not support
> I/O rounding modes.

To my knowledge gfortran supports I/O rounding (for output, not for input)
since 4.5.0. However, it neither supports setting the rounding for run-time
arithmetic via F2003's IEEE modules -- nor does it allow to set the rounding
for compile-time simplification (folding).

Besides, I am not really sure whether rounding is/should work if one ask for
more digits than the number actually has...

(I think Steve talked about the compile-time folding when he wrote: gfortran by
default will use round to nearest.)

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #9 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-02 13:49:15 UTC ---
It seems that MinGW has its own implementation of snprintf called
__mingw_snprintf which can be activated by defining __USE_MINGW_ANSI_STDIO

__mingw_snprintf has the desired behaviour:

#include <stdio.h>
main() {
    double d = .14285714285714286;
    char str[200];
    __mingw_snprintf(str, sizeof(str), "MinGW:  %35.32f", d);
    puts(str);
    _snprintf(str, sizeof(str), "MSVCRT: %35.32f", d);
    puts(str);
    return;
}

MinGW:   0.14285714285714284921269268124888
MSVCRT:  0.14285714285714285000000000000000

Is there a way to use this in Fortran too?

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #10 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-03-02 16:45:26 UTC ---
(In reply to comment #9)
> It seems that MinGW has its own implementation of snprintf called
> __mingw_snprintf [...] __mingw_snprintf has the desired behaviour

>From libgfortran/libgfortran.h:

#if HAVE_MINGW_SNPRINTF
extern int __mingw_snprintf (char *, size_t, const char *, ...)
     __attribute__ ((format (gnu_printf, 3, 4)));
#undef snprintf
#define snprintf __mingw_snprintf
#endif

That should take care of it - but seemingly it does not.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #11 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-03-02 17:03:47 UTC ---
(In reply to comment #9)
> It seems that MinGW has its own implementation of snprintf called
> __mingw_snprintf which can be activated by defining __USE_MINGW_ANSI_STDIO

In MinGW64, __USE_MINGW_ANSI_STDIO is set when _POSIX is specified. In MinGW
that does not seem to be the case and one needs to specify it either directly
or use another define such as _POSIX_SOURCE. Cf.
http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/include/_mingw.h?rev=1.61&cvsroot=src

>From libgfortran/libgfortran.h:

#if defined __MINGW32__
#  define _POSIX 1
#  define gfc_printf gnu_printf
#else

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #12 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-02 17:54:34 UTC ---
(In reply to comment #11)

> From libgfortran/libgfortran.h:
> 
> #if defined __MINGW32__
> #  define _POSIX 1
> #  define gfc_printf gnu_printf
> #else

Since the comment above that reads:

/* Ensure that ANSI conform stdio is used. This needs to be set before
   any system header file is included. */

could it be that it was the intention to set __USE_MINGW_ANSI_STDIO in effect?

So is this a bug in libgfortran.h and it was intended to put
# define _POSIX_SOURCE 1 
there?

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[burnus at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #13 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-03-02 18:02:15 UTC ---
(In reply to comment #12)
> could it be that it was the intention to set __USE_MINGW_ANSI_STDIO in effect?

Yes - and that is what does happen for _POSIX=1 on MinGW64. But it does not
happen for _POSIX=1 on the 32bit MinGW.

Jerry et al.: Do you think the following patch is OK?

--- a/libgfortran/libgfortran.h
+++ b/libgfortran/libgfortran.h
@@ -33,6 +33,7 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If
not, see
    any system header file is included.  */
 #if defined __MINGW32__
 #  define _POSIX 1
+#  define _POSIX_SOURCE 1
 #  define gfc_printf gnu_printf
 #else
 #  define gfc_printf __printf__

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[sgk at troutmask dot apl.washington.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #14 from Steve Kargl <sgk at troutmask dot apl.washington.edu> 2011-03-02 18:17:30 UTC ---
On Wed, Mar 02, 2011 at 06:02:22PM +0000, burnus at gcc dot gnu.org wrote:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945
> 
> --- Comment #13 from Tobias Burnus <burnus at gcc dot gnu.org> 2011-03-02 18:02:15 UTC ---
> (In reply to comment #12)
> > could it be that it was the intention to set __USE_MINGW_ANSI_STDIO in effect?
> 
> Yes - and that is what does happen for _POSIX=1 on MinGW64. But it does not
> happen for _POSIX=1 on the 32bit MinGW.
> 
> Jerry et al.: Do you think the following patch is OK?
> 
> --- a/libgfortran/libgfortran.h
> +++ b/libgfortran/libgfortran.h
> @@ -33,6 +33,7 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If
> not, see
>     any system header file is included.  */
>  #if defined __MINGW32__
>  #  define _POSIX 1
> +#  define _POSIX_SOURCE 1
>  #  define gfc_printf gnu_printf
>  #else
>  #  define gfc_printf __printf__
> 

I though a leading underscore in C denoted the namespace
of the C compiler, and a userspace program should not 
set such a macro.  Guess I'll pull out n1256.pdf later
to read up on the convention.  You may want to ask Joseph.  

It may also be prudent to add

#if define _POSIX_SOURCE
#undef _POSIX_SOURCE
#endif 

prior to (re)setting it to 1.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[jb at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

Janne Blomqvist <jb at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jb at gcc dot gnu.org

--- Comment #15 from Janne Blomqvist <jb at gcc dot gnu.org> 2011-03-03 14:03:40 UTC ---
(In reply to comment #13)
> (In reply to comment #12)
> > could it be that it was the intention to set __USE_MINGW_ANSI_STDIO in effect?
> 
> Yes - and that is what does happen for _POSIX=1 on MinGW64. But it does not
> happen for _POSIX=1 on the 32bit MinGW.
> 
> Jerry et al.: Do you think the following patch is OK?
> 
> --- a/libgfortran/libgfortran.h
> +++ b/libgfortran/libgfortran.h
> @@ -33,6 +33,7 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If
> not, see
>     any system header file is included.  */
>  #if defined __MINGW32__
>  #  define _POSIX 1
> +#  define _POSIX_SOURCE 1
>  #  define gfc_printf gnu_printf
>  #else
>  #  define gfc_printf __printf__

Well, first a disclaimer: I don't use Windows myself, and I have little
experience with win32 programming, and finally, I have no way of testing mingw
specific functionality.

That being said, I'm not sure this patch accomplishes anything. Looking at the
_mingw.h you linked to in another post, __USE_MINGW_ANSI_STDIO is set also if
_GNU_SOURCE is set. Which AFAICS is always set when compiling libgfortran
(there's even a bit of belt-and-suspenders approach here, in that _GNU_SOURCE
is set from both AC_USE_SYSTEM_EXTENSIONS as well as explicitly added to
AM_CFLAGS in Makefile.am; probably the second one can be removed as a
janitorial patch once 4.7 opens).

Secondly, IIRC mingw64 also sets __MINGW32__, so there should be no difference
between mingw32 and ming64 here.

Third, I'm somewhat wary of playing with preprocessor defines this close to a
release. As this issue doesn't seem particularly critical, I'd vote for
postponing it until 4.7. That is, unless Kai or someone else who actually
understands mingw pops in and says it's completely safe and will fix the issue.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #16 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-03 14:58:35 UTC ---
My _mingw.h has the following:

#if defined(_POSIX) && !defined(__USE_MINGW_ANSI_STDIO)
/* Enable __USE_MINGW_ANSI_STDIO if _POSIX defined
 * and If user did _not_ specify it explicitly... */
#  define __USE_MINGW_ANSI_STDIO            1
#endif

So _POSIX_SOURCE really seems to be the wrong way and I wonder what the values
of _POSIX and __USE_MINGW_ANSI_STDIO actually are.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

Thomas Henlich <thenlich at users dot sourceforge.net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
         Resolution|                            |INVALID

--- Comment #17 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-04 06:48:37 UTC ---
The problem has been confirmed to be an issue with the custom snprintf
implementation used in the equation.com distribution and not related to
gfortran/mingw.

In the new snapshot distributed on their website they have fixed this issue,
and are apparently using __mingw_snprintf now, which solves the problem for the
reporter.

Thank you all for your input on this.

Closing.

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #18 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-04 06:53:33 UTC ---
Created attachment 23537
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23537
Test 1 for pedantic IO rounding

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #19 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-04 06:54:03 UTC ---
Created attachment 23538
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23538
Test 2 for pedantic IO rounding

[Bug libfortran/47945] REAL(8) output conversion error on MinGW32

[thenlich at users dot sourceforge.net]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=47945

--- Comment #20 from Thomas Henlich <thenlich at users dot sourceforge.net> 2011-03-04 06:56:56 UTC ---
Added two tests just in case this comes up again.