MinGW-w64 - for 32 and 64 bit Windows / Bugs / #466 pow() function inaccurate in 64 bits

Duncan Murdoch - 2015-03-20

I've looked at the source to pow(), exp2l() and log2l(), and can't spot any changes since version 2.0.9, which we previously used successfully. This suggests that the problem is something else in our gcc build.

Does anyone else with gcc 4.9.2 see this problem on x64?

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Avi - 2015-03-20

For what it is worth, I get the exact same results when compiling newtest.c against the MinGW personal builds of 4.9.2 v 1.7, both SJLJ and SEH. Although, the Rtools executable is about twice the size of the other two. Compiled under Windows 7 64 bit.

gcc (GCC) 4.9.2
Copyright (C) 2014 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

10^j (j=-1), double=0x1.999999999999a0p-4
10^j (j=-2), double=0x1.47ae147ae147c0p-7
10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
10^j (j=-4), double=0x1.a36e2eb1c43300p-14
10^j (j=-5), double=0x1.4f8b588e368f40p-17
10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
10^j (j=-8), double=0x1.5798ee2308c3f0p-27
10^j (j=-9), double=0x1.12e0be826d6990p-30
10^-1, double=0x1.999999999999a0p-4
10^-2, double=0x1.47ae147ae147b0p-7

gcc (x86_64-win32-sjlj-rev1, Built by MinGW-W64 project) 4.9.2
Copyright (C) 2014 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

10^j (j=-1), double=0x1.999999999999a0p-4
10^j (j=-2), double=0x1.47ae147ae147c0p-7
10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
10^j (j=-4), double=0x1.a36e2eb1c43300p-14
10^j (j=-5), double=0x1.4f8b588e368f40p-17
10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
10^j (j=-8), double=0x1.5798ee2308c3f0p-27
10^j (j=-9), double=0x1.12e0be826d6990p-30
10^-1, double=0x1.999999999999a0p-4
10^-2, double=0x1.47ae147ae147b0p-7

gcc (x86_64-win32-seh-rev1, Built by MinGW-W64 project) 4.9.2
Copyright (C) 2014 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

10^j (j=-1), double=0x1.999999999999a0p-4
10^j (j=-2), double=0x1.47ae147ae147c0p-7
10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
10^j (j=-4), double=0x1.a36e2eb1c43300p-14
10^j (j=-5), double=0x1.4f8b588e368f40p-17
10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
10^j (j=-8), double=0x1.5798ee2308c3f0p-27
10^j (j=-9), double=0x1.12e0be826d6990p-30
10^-1, double=0x1.999999999999a0p-4
10^-2, double=0x1.47ae147ae147b0p-7

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Duncan Murdoch - 2015-03-20
  
  On 20/03/2015 3:51 PM, Avi wrote:
  
  For what it is worth, I get the exact same results when compiling
  newtest.c against the MinGW personal builds of 4.9.2 v 1.7, both SJLJ
  and SEH. Although, the Rtools executable is about twice the size of the
  other two. Compiled under Windows 7 64 bit.
  
  You posted 3 sets of output, all the same as far as I can see. Is one
  of them a 32 bit build? My first output is from 64 bit gcc, the second
  from 32 bit gcc, and they don't match. The 32 bit compiler gets correct
  results.
  
  Duncan Murdoch
  
  gcc (GCC) 4.9.2
  Copyright (C) 2014 Free Software Foundation, Inc.
  This is free software; see the source for copying conditions. There is NO
  warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147c0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
  10^j (j=-4), double=0x1.a36e2eb1c43300p-14
  10^j (j=-5), double=0x1.4f8b588e368f40p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
  10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
  10^j (j=-8), double=0x1.5798ee2308c3f0p-27
  10^j (j=-9), double=0x1.12e0be826d6990p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  gcc (x86_64-win32-sjlj-rev1, Built by MinGW-W64 project) 4.9.2
  Copyright (C) 2014 Free Software Foundation, Inc.
  This is free software; see the source for copying conditions. There is NO
  warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147c0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
  10^j (j=-4), double=0x1.a36e2eb1c43300p-14
  10^j (j=-5), double=0x1.4f8b588e368f40p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
  10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
  10^j (j=-8), double=0x1.5798ee2308c3f0p-27
  10^j (j=-9), double=0x1.12e0be826d6990p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  gcc (x86_64-win32-seh-rev1, Built by MinGW-W64 project) 4.9.2
  Copyright (C) 2014 Free Software Foundation, Inc.
  This is free software; see the source for copying conditions. There is NO
  warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147c0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
  10^j (j=-4), double=0x1.a36e2eb1c43300p-14
  10^j (j=-5), double=0x1.4f8b588e368f40p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
  10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
  10^j (j=-8), double=0x1.5798ee2308c3f0p-27
  10^j (j=-9), double=0x1.12e0be826d6990p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  [bugs:#466] http://sourceforge.net/p/mingw-w64/bugs/466 pow()
  function inaccurate in 64 bits
  
  Status: open
  Group: v1.0 (example)
  Created: Thu Mar 19, 2015 07:20 PM UTC by Duncan Murdoch
  Last Updated: Fri Mar 20, 2015 06:21 PM UTC
  Owner: nobody
  
  We're trying to build gcc 4.9.2 using MingW-w64 version 3.3.0, to use in
  a future R release. Our tests are finding that pow(10.0, j) is off by
  one bit in the 64 bit build for j==-2, giving (in hex notation)
  0x1.47ae147ae147c0p-7 when it should be 0x1.47ae147ae147b0p-7.
  
  When the compiler does the calculation, e.g. with pow(10.0, -2.0) coded
  literally, it gets the correct result. This leads to inconsistencies and
  problems.
  
  This seems related to a bug report here:
  https://sourceforge.net/p/tdm-gcc/bugs/219/
  https://sourceforge.net/p/tdm-gcc/bugs/219, but we're not using TDM,
  we're using our own build.
  
  I simplified the example from that report; it's attached as newtest.c.
  Here's the output I'm seeing when I compile in 64 bits using simply
  
  c:/Rtools/gcc492_64/bin/gcc newtest.c -o newtest.exe
  
  $ ./newtest
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147c0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
  10^j (j=-4), double=0x1.a36e2eb1c43300p-14
  10^j (j=-5), double=0x1.4f8b588e368f40p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
  10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
  10^j (j=-8), double=0x1.5798ee2308c3f0p-27
  10^j (j=-9), double=0x1.12e0be826d6990p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  Notice that the second line of output and the last line don't match. If
  I do the same with the 32 bit compiler, I get
  
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147b0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fc0p-10
  10^j (j=-4), double=0x1.a36e2eb1c432d0p-14
  10^j (j=-5), double=0x1.4f8b588e368f10p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed8d0p-20
  10^j (j=-7), double=0x1.ad7f29abcaf480p-24
  10^j (j=-8), double=0x1.5798ee2308c3a0p-27
  10^j (j=-9), double=0x1.12e0be826d6950p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  Notice that most of the 32 bit values differ in the last bits from the
  64 bit values.
  
  Sent from sourceforge.net because you indicated interest in
  https://sourceforge.net/p/mingw-w64/bugs/466/
  https://sourceforge.net/p/mingw-w64/bugs/466
  
  To unsubscribe from further messages, please visit
  https://sourceforge.net/auth/subscriptions/
  https://sourceforge.net/auth/subscriptions
  
  Related
  
  Bugs: #466
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Avi - 2015-03-22

Perhaps I was not clear; they are all 64 bit compilations. Yes, they are all the same—that is the point; the problem you point out with Rtools's version of GCC 4.9.2 seems to exist in MinGW64's own two versions as well. This should confirm your intuition that there is something in the 64 bit code of GCC 4.9.2 itself.

Last edit: Avi 2015-03-22

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Duncan Murdoch - 2015-03-23

I was wrong about the lack of changes -- not sure what files I was looking at.

Our older version mainly used the regular floating point instructions, with 80 bit intermediate values; the new version uses SSE2 instructions, with 64 bit intermediates. That would explain the difference.

I would still call it a bug, since the compiler is apparently using some other precision when computing immediate values.

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Kai Tietz - 2015-03-23

Yes, these differences are caused by x87 80-bit intermediate vs fixed-precision of double/float.
Well, 1 ulp is not much, nevertheless it would be good to find a way to catch it. It might be most likely a rounding-issue. Anyway even with that 1ulp difference, it would be more consistent if we would see 1 ulp difference for both cases, or for none.

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Duncan Murdoch - 2015-03-23
  
  On 23/03/2015 10:26 AM, Kai Tietz wrote:
  
  Yes, these differences are caused by x87 80-bit intermediate vs
  fixed-precision of double/float.
  Well, 1 ulp is not much, nevertheless it would be good to find a way
  to catch it. It might be most likely a rounding-issue. Anyway even
  with that 1ulp difference, it would be more consistent if we would see
  1 ulp difference for both cases, or for none.
  
  Yes, the 1 ulp error is livable, but it does cause problems that pow(10,
  -2) != pow(10, j) with j = -2. I think it would also cause us problems
  if a literal 0.01 was not equal to pow(10, -2), though we could work
  around those more easily.
  
  Duncan Murdoch
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Avi - 2015-03-24
  
  If that is the case, then shouldn't compiling it as "gcc -mfpmath=387
  -fexcess-precision=standard -ffloat-store newtest.c -o newtest.exe"
  provide equality? When using a 64bit toolchain, I still get the
  different values for j = -2. What have I missed?
  
  On Mon, Mar 23, 2015 at 10:26 AM, Kai Tietz ktietz70@users.sf.net wrote:
  
  Yes, these differences are caused by x87 80-bit intermediate vs
  fixed-precision of double/float.
  Well, 1 ulp is not much, nevertheless it would be good to find a way to
  catch it. It might be most likely a rounding-issue. Anyway even with that
  1ulp difference, it would be more consistent if we would see 1 ulp
  difference for both cases, or for none.
  
  [bugs:#466] pow() function inaccurate in 64 bits
  
  Status: open
  Group: v1.0 (example)
  Created: Thu Mar 19, 2015 07:20 PM UTC by Duncan Murdoch
  Last Updated: Mon Mar 23, 2015 02:02 PM UTC
  Owner: nobody
  
  We're trying to build gcc 4.9.2 using MingW-w64 version 3.3.0, to use in a
  future R release. Our tests are finding that pow(10.0, j) is off by one bit
  in the 64 bit build for j==-2, giving (in hex notation)
  0x1.47ae147ae147c0p-7 when it should be 0x1.47ae147ae147b0p-7.
  
  When the compiler does the calculation, e.g. with pow(10.0, -2.0) coded
  literally, it gets the correct result. This leads to inconsistencies and
  problems.
  
  This seems related to a bug report here:
  https://sourceforge.net/p/tdm-gcc/bugs/219/, but we're not using TDM, we're
  using our own build.
  
  I simplified the example from that report; it's attached as newtest.c.
  Here's the output I'm seeing when I compile in 64 bits using simply
  
  c:/Rtools/gcc492_64/bin/gcc newtest.c -o newtest.exe
  
  $ ./newtest
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147c0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
  10^j (j=-4), double=0x1.a36e2eb1c43300p-14
  10^j (j=-5), double=0x1.4f8b588e368f40p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
  10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
  10^j (j=-8), double=0x1.5798ee2308c3f0p-27
  10^j (j=-9), double=0x1.12e0be826d6990p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  Notice that the second line of output and the last line don't match. If I do
  the same with the 32 bit compiler, I get
  
  10^j (j=-1), double=0x1.999999999999a0p-4
  10^j (j=-2), double=0x1.47ae147ae147b0p-7
  10^j (j=-3), double=0x1.0624dd2f1a9fc0p-10
  10^j (j=-4), double=0x1.a36e2eb1c432d0p-14
  10^j (j=-5), double=0x1.4f8b588e368f10p-17
  10^j (j=-6), double=0x1.0c6f7a0b5ed8d0p-20
  10^j (j=-7), double=0x1.ad7f29abcaf480p-24
  10^j (j=-8), double=0x1.5798ee2308c3a0p-27
  10^j (j=-9), double=0x1.12e0be826d6950p-30
  10^-1, double=0x1.999999999999a0p-4
  10^-2, double=0x1.47ae147ae147b0p-7
  
  Notice that most of the 32 bit values differ in the last bits from the 64
  bit values.
  
  Sent from sourceforge.net because you indicated interest in
  https://sourceforge.net/p/mingw-w64/bugs/466/
  
  To unsubscribe from further messages, please visit
  https://sourceforge.net/auth/subscriptions/
  
  Related
  
  Bugs: #466
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
  - Duncan Murdoch - 2015-03-24
    
    On 24/03/2015 12:36 AM, Avi wrote:
    
    If that is the case, then shouldn't compiling it as "gcc -mfpmath=387
    -fexcess-precision=standard -ffloat-store newtest.c -o newtest.exe"
    provide equality? When using a 64bit toolchain, I still get the
    different values for j = -2. What have I missed?
    
    Try it, and let me know if it works.
    
    Duncan Murdoch
    
    On Mon, Mar 23, 2015 at 10:26 AM, Kai Tietz ktietz70@users.sf.net
    ktietz70@users.sf.net wrote:
    
    Yes, these differences are caused by x87 80-bit intermediate vs fixed-precision of double/float. Well, 1 ulp is not much, nevertheless it would be good to find a way to catch it. It might be most likely a rounding-issue. Anyway even with that 1ulp difference, it would be more consistent if we would see 1 ulp difference for both cases, or for none. ------------------------------------------------------------------------ [bugs:#466] <http://sourceforge.net/p/mingw-w64/bugs/466> pow() function inaccurate in 64 bits Status: open Group: v1.0 (example) Created: Thu Mar 19, 2015 07:20 PM UTC by Duncan Murdoch Last Updated: Mon Mar 23, 2015 02:02 PM UTC Owner: nobody We're trying to build gcc 4.9.2 using MingW-w64 version 3.3.0, to use in a future R release. Our tests are finding that pow(10.0, j) is off by one bit in the 64 bit build for j==-2, giving (in hex notation) 0x1.47ae147ae147c0p-7 when it should be 0x1.47ae147ae147b0p-7. When the compiler does the calculation, e.g. with pow(10.0, -2.0) coded literally, it gets the correct result. This leads to inconsistencies and problems. This seems related to a bug report here: https://sourceforge.net/p/tdm-gcc/bugs/219/ <https://sourceforge.net/p/tdm-gcc/bugs/219>, but we're not using TDM, we're using our own build. I simplified the example from that report; it's attached as newtest.c. Here's the output I'm seeing when I compile in 64 bits using simply c:/Rtools/gcc492_64/bin/gcc newtest.c -o newtest.exe $ ./newtest 10^j (j=-1), double=0x1.999999999999a0p-4 10^j (j=-2), double=0x1.47ae147ae147c0p-7 10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10 10^j (j=-4), double=0x1.a36e2eb1c43300p-14 10^j (j=-5), double=0x1.4f8b588e368f40p-17 10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20 10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24 10^j (j=-8), double=0x1.5798ee2308c3f0p-27 10^j (j=-9), double=0x1.12e0be826d6990p-30 10^-1, double=0x1.999999999999a0p-4 10^-2, double=0x1.47ae147ae147b0p-7 Notice that the second line of output and the last line don't match. If I do the same with the 32 bit compiler, I get 10^j (j=-1), double=0x1.999999999999a0p-4 10^j (j=-2), double=0x1.47ae147ae147b0p-7 10^j (j=-3), double=0x1.0624dd2f1a9fc0p-10 10^j (j=-4), double=0x1.a36e2eb1c432d0p-14 10^j (j=-5), double=0x1.4f8b588e368f10p-17 10^j (j=-6), double=0x1.0c6f7a0b5ed8d0p-20 10^j (j=-7), double=0x1.ad7f29abcaf480p-24 10^j (j=-8), double=0x1.5798ee2308c3a0p-27 10^j (j=-9), double=0x1.12e0be826d6950p-30 10^-1, double=0x1.999999999999a0p-4 10^-2, double=0x1.47ae147ae147b0p-7 Notice that most of the 32 bit values differ in the last bits from the 64 bit values. ------------------------------------------------------------------------ Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/mingw-w64/bugs/466/ <https://sourceforge.net/p/mingw-w64/bugs/466> To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/ <https://sourceforge.net/auth/subscriptions>
    
    [bugs:#466] http://sourceforge.net/p/mingw-w64/bugs/466 pow()
    function inaccurate in 64 bits
    
    Status: open
    Group: v1.0 (example)
    Created: Thu Mar 19, 2015 07:20 PM UTC by Duncan Murdoch
    Last Updated: Mon Mar 23, 2015 04:08 PM UTC
    Owner: nobody
    
    We're trying to build gcc 4.9.2 using MingW-w64 version 3.3.0, to use in
    a future R release. Our tests are finding that pow(10.0, j) is off by
    one bit in the 64 bit build for j==-2, giving (in hex notation)
    0x1.47ae147ae147c0p-7 when it should be 0x1.47ae147ae147b0p-7.
    
    When the compiler does the calculation, e.g. with pow(10.0, -2.0) coded
    literally, it gets the correct result. This leads to inconsistencies and
    problems.
    
    This seems related to a bug report here:
    https://sourceforge.net/p/tdm-gcc/bugs/219/
    https://sourceforge.net/p/tdm-gcc/bugs/219, but we're not using TDM,
    we're using our own build.
    
    I simplified the example from that report; it's attached as newtest.c.
    Here's the output I'm seeing when I compile in 64 bits using simply
    
    c:/Rtools/gcc492_64/bin/gcc newtest.c -o newtest.exe
    
    $ ./newtest
    10^j (j=-1), double=0x1.999999999999a0p-4
    10^j (j=-2), double=0x1.47ae147ae147c0p-7
    10^j (j=-3), double=0x1.0624dd2f1a9fd0p-10
    10^j (j=-4), double=0x1.a36e2eb1c43300p-14
    10^j (j=-5), double=0x1.4f8b588e368f40p-17
    10^j (j=-6), double=0x1.0c6f7a0b5ed900p-20
    10^j (j=-7), double=0x1.ad7f29abcaf4e0p-24
    10^j (j=-8), double=0x1.5798ee2308c3f0p-27
    10^j (j=-9), double=0x1.12e0be826d6990p-30
    10^-1, double=0x1.999999999999a0p-4
    10^-2, double=0x1.47ae147ae147b0p-7
    
    Notice that the second line of output and the last line don't match. If
    I do the same with the 32 bit compiler, I get
    
    10^j (j=-1), double=0x1.999999999999a0p-4
    10^j (j=-2), double=0x1.47ae147ae147b0p-7
    10^j (j=-3), double=0x1.0624dd2f1a9fc0p-10
    10^j (j=-4), double=0x1.a36e2eb1c432d0p-14
    10^j (j=-5), double=0x1.4f8b588e368f10p-17
    10^j (j=-6), double=0x1.0c6f7a0b5ed8d0p-20
    10^j (j=-7), double=0x1.ad7f29abcaf480p-24
    10^j (j=-8), double=0x1.5798ee2308c3a0p-27
    10^j (j=-9), double=0x1.12e0be826d6950p-30
    10^-1, double=0x1.999999999999a0p-4
    10^-2, double=0x1.47ae147ae147b0p-7
    
    Notice that most of the 32 bit values differ in the last bits from the
    64 bit values.
    
    Sent from sourceforge.net because you indicated interest in
    https://sourceforge.net/p/mingw-w64/bugs/466/
    https://sourceforge.net/p/mingw-w64/bugs/466
    
    To unsubscribe from further messages, please visit
    https://sourceforge.net/auth/subscriptions/
    https://sourceforge.net/auth/subscriptions
    
    Related
    
    Bugs: #466
    
    If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Avi - 2015-03-23

Is it possible that the difference is that per https://gcc.gnu.org/onlinedocs/gcc-4.9.2/gcc/i386-and-x86-64-Options.html#i386-and-x86-64-Options the GCC default for i386 is -mfpmath=387 but for x86-64 the default is -mfpmath=sse?

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

I did, it doesn't help, still the c/b difference (I did say "I still get the different values for j = -2.").

Last time I used assembler was over 20 years ago and for the 80386, so I don't know what the differences mean, but below is the differences between the assembler for the case where "-mfpmath=387 -fexcess-precision=standard -ffloat-store" was thrown, and where it wasn't (> = flags thrown; < flags not thrown):

19,20c19,20
<   subq    $64, %rsp
<   .seh_stackalloc 64
---
>   subq    $80, %rsp
>   .seh_stackalloc 80
26,27c26,28
<   pxor    %xmm0, %xmm0
<   cvtsi2sd    -4(%rbp), %xmm0
---
>   fildl   -4(%rbp)
>   fstpl   -16(%rbp)
>   movq    -16(%rbp), %rdx
29,31c30,33
<   movapd  %xmm0, %xmm1
<   movq    %rax, -24(%rbp)
<   movsd   -24(%rbp), %xmm0
---
>   movq    %rdx, -40(%rbp)
>   movsd   -40(%rbp), %xmm1
>   movq    %rax, -40(%rbp)
>   movsd   -40(%rbp), %xmm0
34,36c36,38
<   movq    %rax, -16(%rbp)
<   movq    -16(%rbp), %rcx
<   movq    -16(%rbp), %rax
---
>   movq    %rax, -24(%rbp)
>   movq    -24(%rbp), %rcx
>   movq    -24(%rbp), %rax
38,39c40,41
<   movq    %rcx, -24(%rbp)
<   movsd   -24(%rbp), %xmm2
---
>   movq    %rcx, -40(%rbp)
>   movsd   -40(%rbp), %xmm2
48,52c50,54
<   movq    %rax, -16(%rbp)
<   movq    -16(%rbp), %rdx
<   movq    -16(%rbp), %rax
<   movq    %rdx, -24(%rbp)
<   movsd   -24(%rbp), %xmm1
---
>   movq    %rax, -24(%rbp)
>   movq    -24(%rbp), %rdx
>   movq    -24(%rbp), %rax
>   movq    %rdx, -40(%rbp)
>   movsd   -40(%rbp), %xmm1
57,61c59,63
<   movq    %rax, -16(%rbp)
<   movq    -16(%rbp), %rdx
<   movq    -16(%rbp), %rax
<   movq    %rdx, -24(%rbp)
<   movsd   -24(%rbp), %xmm1
---
>   movq    %rax, -24(%rbp)
>   movq    -24(%rbp), %rdx
>   movq    -24(%rbp), %rax
>   movq    %rdx, -40(%rbp)
>   movsd   -40(%rbp), %xmm1
66c68
<   addq    $64, %rsp
---
>   addq    $80, %rsp

Last edit: Avi 2015-03-24

Dominick Samperi - 2015-04-18

This bug report suggests that the rounding problem is specific to GCC 4.9.2, but I see the same problem (1-bit error when j == -2) with GCC 4.6.3 (from Rtools 3.3.0.1958). The attached modification of newtest.c contains a specialization of pow() for integral powers named pown() that is used instead of pow(). It solves the 1-bit error problem. When used with C++ (i.e., g++) where overloading is supported, pown() can be renamed to pow().

newtestwithpown.c

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

jeroen - 2015-09-16

KK has found that powl() gives the correct result on win64. So this is an alternative workaround.

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

pow() function inaccurate in 64 bits

A complete runtime environment for gcc

Group

Searches

Help

#466 pow() function inaccurate in 64 bits

Related

Discussion

Related

Related

Related