[Bug fortran/40070] New: Some math expressions containing exponents fail on a Windows 64 build

[Bug fortran/40070] New: Some math expressions containing exponents fail on a Windows 64 build

[jonathan dot d dot wetherbee at saic dot com]

Some equations using exponents on Windows XP Professional 64-bit evaluate
incorrectly, but can be computed correctly by splitting up the equation into
multiple statements.

I apologize in advance if this turns out to be just a MinGW bug, but thought
I'd start with gcc first.

I am using MinGW GFortran with the following built-in specs.
Target: x86_64-pc-mingw32
Configured with: ../gcc/configure --target=x86_64-pc-mingw32
--prefix=/c/buildbot/vista64-mingw32/mingw-x86/build/root
--with-sysroot=/c/buildbot/vista64-mingw32/mingw-x86/build/root
--disable-multilib --enable-languages=all,obj-c++
Thread model: win32
gcc version 4.4.0 20090110 (experimental) (GCC)


Here is some sample code that demonstrates the problem:
myfunc.f
------------------------------------------------------------
      Subroutine myfunc(a,b,c,d,e,res1,res2)

      Real a,b,c,d,e
      Parameter (g = 273.15)


      res1 = a*(b/a)**((c-d-g)/(e-d))

      res2 = (b/a)**((c-d-g)/(e-d))
      res2 = res2 * a

      Return
      End


      Subroutine myfunc2(a,b,c,d,res1,res2)

      Real a,b,c,d,res1,res2

      res1 = a*(b/c)**d
      res2 = (b/c)**d
      res2 = a * res2

      Return
      End

      Subroutine myfunc3(a,b,c,d,e,f,g,h,i,res1,res2)

      Real a,b,c,d,e,f,g,h,i,res1,res2

      res1 = a*((b/(b+c*(d-e)))**((f*g)/(h*i)))*100.
      res2 = ((b/(b+c*(d-e)))**((f*g)/(h*i)))*100.
      res2 = a*res2

      Return
      End




main.c
------------------------------------------------------

#include <stdio.h>


extern void myfunc(float *a, float *b, float *c, float *d, float *e, float
*res1, float *res2); 
extern void myfunc2(float *a, float *b, float *c, float *d, float *res1, float
*res2); 
extern void myfunc3(float *a, float *b, float *c, float *d, float *e, float *f,
float *g, float *h, float *i, float *res1, float *res2); 

int main( int argc, char *argv[] )
{
    float a,b,c,d,e,f,g,h,i,res1,res2;

    a=1000.0;
    b=850.0;
    c=143.204;
    d=30.5280;
    e=8.2351;
    res1 = 0;
    res2 = 0;

    myfunc(&a,&b,&c,&d,&e,&res1,&res2);

    printf("Function 1:\n--------------------\nresult 1: %f\nresult 2:
%f\n\n",res1,res2);


    a=103780.0;
    b=253.13005;
    c=231.676270;
    d=-16.904803;

    myfunc2(&a,&b,&c,&d,&res1,&res2);

    printf("Function 2:\n--------------------\nresult 1: %f\nresult 2:
%f\n\n",res1,res2);


    a=1513.25;
    b=238.149994;
    c=-16.500040;
    d=3.2;
    e=0;
    f=9.806650;
    g=20.964399;
    h=7614.320313;
    i=-0.1065;

    myfunc3(&a,&b,&c,&d,&e,&f,&g,&h,&i,&res1,&res2);

    printf("Function 3:\n--------------------\nresult 1: %f\nresult 2:
%f\n\n",res1,res2);

    return 0;

}



--------------------------------------------

Compiled with the following options:
-fno-uncerscoring
-fno-range-check
-W
-std=f95


-- 
           Summary: Some math expressions containing exponents fail on a
                    Windows 64 build
           Product: gcc
           Version: 4.4.0
            Status: UNCONFIRMED
          Severity: blocker
          Priority: P3
         Component: fortran
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: jonathan dot d dot wetherbee at saic dot com


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[dominiq at lps dot ens dot fr]

------- Comment #1 from dominiq at lps dot ens dot fr  2009-05-08 14:55 -------
Why are you using a code in C for the tests? If I use the following as the main
program

      real res1, res2
      call myfunc(1000.0,850.0,143.204,30.5280,8.2351,res1,res2)
      print *, res1, res1-res2

      call myfunc2(103780.0,253.13005,231.676270,-16.904803,res1,res2)
      print *, res1, res1-res2

      call
myfunc3(1513.25,238.149994,-16.500040,3.2,0.0,9.806650,20.964399,7614.320313,-0.1065,res1,res2)
      print *, res1, res1-res2, spacing(res1)

I get

   310.40363       0.0000000    
   23223.025       0.0000000    
   142007.84      1.56250000E-02  1.56250000E-02

So the difference between the res1 and res2 is either 0  or one LSB. What did
you expect?

Tests on intel-darwin with or without -m64.


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jonathan dot d dot wetherbee at saic dot com]

------- Comment #2 from jonathan dot d dot wetherbee at saic dot com  2009-05-08 15:03 -------
I used C code to call the function because the fortran code is being compiled
into a C based application.

When I run the sample main function, the values are obviously different.

Function 1:
result 1: 0.000
result 2: 310.403625

Function 2:
result 1: 0.000
result 2: 23223.025391

Function 3:
result 1: -41718457312547043000000000000000000.0000
result 2: 142007.828125


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[kargl at gcc dot gnu dot org]

------- Comment #3 from kargl at gcc dot gnu dot org  2009-05-08 15:07 -------
Fortran bug reports never have a Severity of Blocker unless
a GCC Release Manager sets that value.  Reset to normal.

You forgot to include the actual output you got from
running your code.

You forgot to include the output you expected to get
from the code.


-- 

kargl at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Severity|blocker                     |normal
           Priority|P3                          |P4


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jonathan dot d dot wetherbee at saic dot com]

------- Comment #4 from jonathan dot d dot wetherbee at saic dot com  2009-05-08 15:27 -------
When I run Dominique's main function written in fortran, I get:

res1        res2       res1-res2
 6.09E-24   310.40363  -310.40363
 8.92E-23   23223.025  -23223.025
-4.01E-06   142007.59  -142007.59

Oddly enough, just about every other line in this 20,000 line fortran code
works fine.  It is starting to look like it might just be my build of gcc.

I'm running it on an Intel Xeon E5450.


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[dfranke at gcc dot gnu dot org]

------- Comment #5 from dfranke at gcc dot gnu dot org  2009-05-08 15:30 -------
You use 'float' in the C routines and 'REAL' in the fortran code. What their
respective size? Is float always 4 byte in C? I think default-real-kind on
win64 could be 8, not 4? Could also explain why it worked for Dominique, using
the same kind everywhere?!

You could try and see if using

  USE ISO_C_BINDING
  REAL(C_FLOAT)

instead of plain REAL makes a difference?


-- 

dfranke at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |dfranke at gcc dot gnu dot
                   |                            |org


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[dominiq at lps dot ens dot fr]

------- Comment #6 from dominiq at lps dot ens dot fr  2009-05-08 15:36 -------
Adding the print *, ... in the subroutines and mixing C and Fortran, I get:

   310.40363       0.0000000    
Function 1:
--------------------
result 1: 310.403625
result 2: 310.403625

   23223.025       0.0000000    
Function 2:
--------------------
result 1: 23223.025391
result 2: 23223.025391

   142007.84      1.56250000E-02  1.56250000E-02
Function 3:
--------------------
result 1: 142007.843750
result 2: 142007.828125

> I think default-real-kind on win64 could be 8

It would as weird as it was with Crays. default-real-kind is usually the
shortest float available, i.e. real(4) on intel procs.


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jonathan dot d dot wetherbee at saic dot com]

------- Comment #7 from jonathan dot d dot wetherbee at saic dot com  2009-05-08 15:40 -------
Checked the size of a float with sizeof(float) and got 4 bytes.  How can I
check the size of REAL in fortran?

When I try the USE ISO_C_BINDING I get the same result.


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[kargl at gcc dot gnu dot org]

------- Comment #8 from kargl at gcc dot gnu dot org  2009-05-08 15:42 -------
(In reply to comment #0)

>       Subroutine myfunc(a,b,c,d,e,res1,res2) 
>       Real a,b,c,d,e
>       Parameter (g = 273.15)
>       res1 = a*(b/a)**((c-d-g)/(e-d))
>       res2 = (b/a)**((c-d-g)/(e-d))
>       res2 = res2 * a
>       Return
>       End

Can you compile the above subroutine with -fdump-tree-original
and post the *.original file?

>     a=1513.25;
>     b=238.149994;
>     c=-16.500040;
>     d=3.2;
>     e=0;
>     f=9.806650;
>     g=20.964399;
>     h=7614.320313;
>     i=-0.1065;

Have you read Goldberg's paper?

What does 

   program a
   print *, precision(1.e0), digits(1.e0)
   end program a

produce?


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[dominiq at lps dot ens dot fr]

------- Comment #9 from dominiq at lps dot ens dot fr  2009-05-08 15:49 -------
Couldn't it be an alignement problem? 


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jonathan dot d dot wetherbee at saic dot com]

------- Comment #10 from jonathan dot d dot wetherbee at saic dot com  2009-05-08 15:51 -------
When I run the fortran code: print *, precision(1.e0), digits(1.e0)
I get: 6,24


*.original file yields:

myfunc (real(kind=4) & a, real(kind=4) & b, real(kind=4) & c, real(kind=4) & d,
real(kind=4) & e, real(kind=4) & res1, real(kind=4) & res2)
{
  {
    real(kind=4) D.1529;

    D.1529 = ((*b / *a));
    *res1 = *a * __builtin_powf (D.1529, (((((*c - *d) -
2.73149993896484375e+2)) / ((*e - *d)))));
  }
  {
    real(kind=4) D.1530;

    D.1530 = ((*b / *a));
    *res2 = __builtin_powf (D.1530, (((((*c - *d) - 2.73149993896484375e+2)) /
((*e - *d)))));
  }
  *res2 = *res2 * *a;
  goto __return_myfunc;
  __return_myfunc:;
}


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[kargl at gcc dot gnu dot org]

------- Comment #11 from kargl at gcc dot gnu dot org  2009-05-08 16:18 -------
(In reply to comment #10)
> When I run the fortran code: print *, precision(1.e0), digits(1.e0)
> I get: 6,24

Read Goldberg's paper.

I suspect you have either a rounding issue
(hint, see Goldberg, the number 6 and your
input list) or whatever __builtin_powf maps
to on your system is broken?  I lay my money
on the former.


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jb at gcc dot gnu dot org]

------- Comment #12 from jb at gcc dot gnu dot org  2009-05-20 08:37 -------
Both the mixed C/Fortran and the pure Fortran version by Dominique works as
expected for me on x86_64-linux. I.e. I get the same results as reported by
Dominique in comments #1 and #6.

This looks like a target bug.


-- 

jb at gcc dot gnu dot org changed:

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


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[dominiq at lps dot ens dot fr]

------- Comment #13 from dominiq at lps dot ens dot fr  2009-05-20 08:55 -------
> This looks like a target bug.

This looks also like a problem with the way arguments and results are handled
(this is why I have suggested an alignement problem in comment #9).

The result in comment #2:

Function 3:
result 1: -41718457312547043000000000000000000.0000
result 2: 142007.828125

seems also to point in this direction.


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[ktietz at gcc dot gnu dot org]

------- Comment #14 from ktietz at gcc dot gnu dot org  2009-09-05 09:17 -------
(In reply to comment #13)
> > This looks like a target bug.
> 
> This looks also like a problem with the way arguments and results are handled
> (this is why I have suggested an alignement problem in comment #9).
> 
> The result in comment #2:
> 
> Function 3:
> result 1: -41718457312547043000000000000000000.0000
> result 2: 142007.828125
> 
> seems also to point in this direction.
> 

We had a strict-aliasing bug in mingw's runtime. Could you retest the bug with
newer runtime, if it is solved.

Cheers,
Kai


-- 

ktietz at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |ktietz at gcc dot gnu dot
                   |                            |org
 GCC target triplet|                            |*-*-mingw32


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[dfranke at gcc dot gnu dot org]

-- 

dfranke at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |WAITING


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jvdelisle at gcc dot gnu dot org]

------- Comment #15 from jvdelisle at gcc dot gnu dot org  2010-02-04 01:12 -------
Any further word on this?


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[ktietz at gcc dot gnu dot org]

------- Comment #16 from ktietz at gcc dot gnu dot org  2010-02-04 08:47 -------
(In reply to comment #15)
> Any further word on this?

As I said in comment #14, we fixed a strict-aliasing bug in our C runtime
related to POSIX printf. As I tested it with current runtime, result looks ok
to me. It isn't an alignment issue, and neiter it was a gcc-bug.

Kai


-- 


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

[Bug fortran/40070] Some math expressions containing exponents fail on a Windows 64 build

[jvdelisle at gcc dot gnu dot org]

------- Comment #17 from jvdelisle at gcc dot gnu dot org  2010-02-05 03:04 -------
Closing. Not a gcc/gfortran bug.


-- 

jvdelisle at gcc dot gnu dot org changed:

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


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