On win32, np.absolute changes floating point control word

[pv]

Consider this:

import numpy as np
import fpread
print("before:", hex(fpread._control87() & fpread._MCW_PC))
np.absolute(complex(np.inf, 0))
print("after:", hex(fpread._control87() & fpread._MCW_PC))

where _control87 returns the floating point control word. Full code in https://github.com/pv/fpread

On win32, using the current numpy 1.13.1 Python 3.6 wheels from Pypi, this prints

before: 0x10000
after: 0x0

I.e., the precision mode changes from 53-bit to 64-bit. On win64 it's 0x0 before and after.

The mode is apparently not changed unless the value passed in to np.absolute is both complex-valued and contains nan or inf.

In addition to np.absolute, also np.hypot(np.inf, 0) appears to have the same issue.

I'm not sure if this is actually a bug, but it seems strange that np.absolute changes the precision mode. This however appears to be a source for some fun as it seems as if this can break gfortran optimizations and make heisenbugs appear.

[njsmith]

Yeah, numpy shouldn't be changing the precision like this. That's very weird that it is.

It looks like np.absolute on complex just calls npy_hypot, which is selected at compile time between the platform hypot or our implementation. And our implementation starts with some checks to bail out if npy_isinf or npy_isnan return true. So it sounds like either a bug in the platform hypot, or a bug in npy_isinf/npy_isnan, which themselves have several implementations.

I guess the next step might be to look at the disassembly for CDOUBLE_absolute in your build to see what it's actually doing, and then following the trail from there?

Are you using the official numpy windows wheels? [edit: sorry, I see now that you said you are.] Is this consistent on both Python 2 and 3? (Wondering if the different MSVC versions might affect anything.)

[ghost]

Are the official binaries compiled with MSVC or mingwpy?

[njsmith]

MSVC. Build scripts, Appveyor account, sample appveyor log containing lots of calls to cl.exe

[pv]

The numpy-1.13.1-cp27-none-win32.whl wheel does not change the precision mode.
I remember seeing this also with locally compiled numpy (py36 + MSVC14).

I now have 5 different C compilers, 3 different unix environments. Now I just need a third debugger to proceed. Lovely platform.

[ghost]

I have MSVC 2017 installed, but here's my contribution:

PS C:\Users\User\Downloads> type test.c
// crt_hypot.c  
// This program prints the hypotenuse of a right triangle.  
  
#include <math.h>  
#include <stdio.h>  

#pragma fenv_access (on)
  
int main( void )  
{  
   double x = INFINITY, y = 4.0;  

   printf("FPU mode: %2.1d\n", _control87(0, 0));
   printf( "If a right triangle has sides %2.1f and %2.1f, "  
           "its hypotenuse is %2.1f\n", x, y, _hypot( x, y ) );  
   
   printf("FPU mode: %2.1d", _control87(0, 0)); 
}  
PS C:\Users\User\Downloads> .\test
FPU mode: 589855
If a right triangle has sides inf and 4.0, its hypotenuse is inf
FPU mode: 786463
PS C:\Users\User\Downloads>

[njsmith]

Odd. objdump says that the umath.cp36-win32.pyd in numpy-1.13.1-cp36-none-win32.whl is importing the hypot symbol from api-ms-win-crt-math-l1-1-0.dll. I'm not quite sure how to check if it's actually being used, though. If you have a debugger on windows then I guess you could set a breakpoint on it and see if np.absolute triggers it.

@xoviat: Try with a nan or inf?

[njsmith]

The 3.6 wheels are built with MSVC 2015, but MSVC 2015 and MSVC 2017 both use the hypot implementation from the UCRT, which is shipped as part of the operating system. So in theory MSVC 2015 and MSVC 2017 should behave the same way here. (And the 2.7 wheels are using a much older version of MSVC that doesn't use the UCRT.)

[ghost]

Jackpot:

FPU mode: 589855
If a right triangle has sides inf and 4.0, its hypotenuse is inf
FPU mode: 786463

[ghost]

@brettcannon Who should we contact about this?

[njsmith]

What an odd bug.

@brettcannon: to save you some reading: the bug is that with the 32-bit UCRT, calling hypot(inf, 1.0) apparently causes the x87 control word register to flip from regular double-precision mode – which is what it's always supposed to be set to on Windows – to "extended" precision mode – Intel's weird 80-bit float mode, which is never used on Windows.

CC @zooba as well.

[pv]

I guess Numpy should just blacklist hypot on 32-bit MSVC? If so, done in gh-9574. It seems hypot() is indeed the only thing that does this stuff, as far as the test suite can see.

[zooba]

Sounds reasonable. The implementation of hypot certainly saves and restores the control word when there is an SNaN, both arguments are zero, or there's a special value that isn't inf, SNaN or QNaN. Without digging into 20+ year old macros I wouldn't be able to figure out why it flips that bit though, and blacklisting sounds like the right approach anyway.

[pv]

Does someone want to submit a MSVC bug report (I'm not up for that.)? Do they take bugreports?

[zooba]

Yep, we do. I'll email the team with a link to this issue.

[pv]

Also cabs* seem to have this issue.