BUG: float64 nanmin() returns NaN on little-endian MIPS (mipsel, mips64el)

[stefanor]

Describe the issue:

Build failures of silx on debian logs point to a problem with numpy.nanmin() for float64 on little endian MIPS platforms.

Reproduce the code example:

import numpy
data = numpy.array((float('nan'), 1.0), dtype='float64')
minimum = numpy.nanmin(data)
print(f"minimum={minimum}")
assert minimum == 1.0

Error message:

minimum=nan
Traceback (most recent call last):
  File "/tmp/testcase.py", line 6, in <module>
    assert minimum == 1.0
           ^^^^^^^^^^^^^^
AssertionError

Runtime information:

import sys, numpy; print(numpy.version); print(sys.version)
1.24.1
3.11.1 (main, Dec 31 2022, 10:23:59) [GCC 12.2.0]
print(numpy.show_runtime())
WARNING: threadpoolctl not found in system! Install it by pip install threadpoolctl. Once installed, try np.show_runtime again for more detailed build information
[{'simd_extensions': {'baseline': [], 'found': [], 'not_found': []}}]
None

Context for the issue:

This is a fundamental bug in a data type, breaking higher level verification tests in other packages.

[tylerjereddy]

Is there a GCC compile farm machine with the appropriate configuration to debug this? I checked https://cfarm.tetaneutral.net/machines/list/ and tried host gcc230, which also runs Debian, but it is Big Endian it seems:

Architecture:        mips64
Byte Order:          Big Endian
CPU(s):              4
On-line CPU(s) list: 0-3
Thread(s) per core:  1
Core(s) per socket:  4
Socket(s):           1
NUMA node(s):        1
BogoMIPS:            2000.00
L1d cache:           32K
L1i cache:           78K
L2 cache:            512K
NUMA node0 CPU(s):   0-3

[stefanor]

I'm not aware of any public machines, I'm afraid. I can test things on Debian's mips porterbox, if that helps. I haven't tried to debug this, at all, beyond pinpointing the issue in numpy.

I believe that QEMU's support for mips64el and mipsel are both usable.

[seberg]

I tried to have a bit of a look, but need to try to get emulation running. It does seem like this this should just end up using fmin in C, unless it goes down some interesting SIMD paths (but I don't think it does from np.show_runtime()). It would be surprising to fail though.

[seberg]

@stefanor any chance you could make sure that fmin itself isn't broken? My guess is, something like this should do:

#include <math.h>
#include <stdio.h>

int main() {
    volatile double f1 = 1.;
    volatile double f2 = NAN;

    printf("Results are: %f, %f\n", fmin(f1, f2), fmin(f2, f1));
    return 0;
}

Compile with gcc -lm and run. I really don't see how we will not end up with fmin.

[stefanor]

On both mips64el and mipsel:

Results are: 1.000000, 1.000000

[seberg]

Whenever I return to this, I don't really see how this can happen (sorry I didn't get to a working qemu/native setup, yet).

Do you have the NumPy build setup/flags ready or are they the defaults? Maybe there is something simple like someone being a bit overzealous on optimization and enabling -ffast-math?!
(In godbold, -ffast-math seems to replace the fmin call with __ledf2 which is not NaN aware, there may be other ways to trigger that I guess. EDIT: I doubt it, but am not sure some of the following instructions doesn't do something about it.)

[stefanor]

Looks like it's all defaults.

You can see them in a build log: https://buildd.debian.org/status/fetch.php?pkg=numpy&arch=mips64el&ver=1%3A1.24.2-1&stamp=1675927314&raw=0

[mattip]

I see gcc12 and these flags:

mips64el-linux-gnuabi64-gcc -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -O2 -ffile-prefix-map=/<<PKGBUILDDIR>>=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC

Do we have any known successful little-endian systems?

[stefanor]

Do we have any known successful little-endian systems?

We don't see the issue on any other Debian platforms: https://buildd.debian.org/status/logs.php?pkg=silx&ver=1.1.0%2Bdfsg-4 (there is a test failure on s390x, but that's unrelated)

[mattip]

Thanks. Looking for instance at the sucessful ppc64el build log, I see gcc12 but not the specific flags used. With that, I do not see any flag in the mips64el build that should cause problems.

[seiko2plus]

Maybe disabling all kinds of optimization including disabling optimization level 3 can be used as a workaround for this issue through the following command:

python setup.py build --disable-optimization install --user 

through pip:

pip install --no-use-pep517 --global-option=build \
--global-option="--disable-optimization" ./

This issue may be related to the following manual unroll which can be disabled by build option --disable-optimization:

numpy/numpy/core/src/umath/loops_minmax.dispatch.c.src

Lines 368 to 436 in bc8cf11

	#ifndef NPY_DISABLE_OPTIMIZATION
	// scalar unrolls
	if (IS_BINARY_REDUCE) {
	// Note, 8x unroll was chosen for best results on Apple M1
	npy_intp elemPerLoop = 8;
	if((i+elemPerLoop) <= len){
	@type@ m0 = *((@type@ *)(ip2 + (i + 0) * is2));
	@type@ m1 = *((@type@ *)(ip2 + (i + 1) * is2));
	@type@ m2 = *((@type@ *)(ip2 + (i + 2) * is2));
	@type@ m3 = *((@type@ *)(ip2 + (i + 3) * is2));
	@type@ m4 = *((@type@ *)(ip2 + (i + 4) * is2));
	@type@ m5 = *((@type@ *)(ip2 + (i + 5) * is2));
	@type@ m6 = *((@type@ *)(ip2 + (i + 6) * is2));
	@type@ m7 = *((@type@ *)(ip2 + (i + 7) * is2));
	i += elemPerLoop;
	for(; (i+elemPerLoop)<=len; i+=elemPerLoop){
	@type@ v0 = *((@type@ *)(ip2 + (i + 0) * is2));
	@type@ v1 = *((@type@ *)(ip2 + (i + 1) * is2));
	@type@ v2 = *((@type@ *)(ip2 + (i + 2) * is2));
	@type@ v3 = *((@type@ *)(ip2 + (i + 3) * is2));
	@type@ v4 = *((@type@ *)(ip2 + (i + 4) * is2));
	@type@ v5 = *((@type@ *)(ip2 + (i + 5) * is2));
	@type@ v6 = *((@type@ *)(ip2 + (i + 6) * is2));
	@type@ v7 = *((@type@ *)(ip2 + (i + 7) * is2));
	m0 = SCALAR_OP(m0, v0);
	m1 = SCALAR_OP(m1, v1);
	m2 = SCALAR_OP(m2, v2);
	m3 = SCALAR_OP(m3, v3);
	m4 = SCALAR_OP(m4, v4);
	m5 = SCALAR_OP(m5, v5);
	m6 = SCALAR_OP(m6, v6);
	m7 = SCALAR_OP(m7, v7);
	}
	m0 = SCALAR_OP(m0, m1);
	m2 = SCALAR_OP(m2, m3);
	m4 = SCALAR_OP(m4, m5);
	m6 = SCALAR_OP(m6, m7);
	m0 = SCALAR_OP(m0, m2);
	m4 = SCALAR_OP(m4, m6);
	m0 = SCALAR_OP(m0, m4);
	*((@type@ *)op1) = SCALAR_OP(*((@type@ *)op1), m0);
	}
	} else{
	// Note, 4x unroll was chosen for best results on Apple M1
	npy_intp elemPerLoop = 4;
	for(; (i+elemPerLoop)<=len; i+=elemPerLoop){
	/* Note, we can't just load all, do all ops, then store all here.
	* Sometimes ufuncs are called with `accumulate`, which makes the
	* assumption that previous iterations have finished before next
	* iteration. For example, the output of iteration 2 depends on the
	* result of iteration 1.
	*/
	/**begin repeat2
	* #unroll = 0, 1, 2, 3#
	*/
	@type@ v@unroll@ = *((@type@ *)(ip1 + (i + @unroll@) * is1));
	@type@ u@unroll@ = *((@type@ *)(ip2 + (i + @unroll@) * is2));
	*((@type@ *)(op1 + (i + @unroll@) * os1)) = SCALAR_OP(v@unroll@, u@unroll@);
	/**end repeat2**/
	}
	}
	#endif // NPY_DISABLE_OPTIMIZATION