Chapel: nans do not propagate through a reduction

[npadmana]

Summary of Problem

nan's do not appear to propagate through a reduction.

Steps to Reproduce

Source Code:

var a = [1.0e150,2.0e309];
var b = [2.0e150,2.1e309];
writeln(max reduce abs(a-b));
writeln(abs(a-b));

produces

1e+150
1e+150 nan

whereas I would expect this to produce a nan.

This is on Chapel 1.19 (running on tio.run).

[vasslitvinov]

@npadmana how would you implement a max-reduction so that it returns a nan in this case?

We implement max-reduction using a set of "greater than" comparisons over the input numbers. Since nan never compares "greater than" any value (in Chapel at least), so it never raises to the top. For example:

writeln(max(NAN,1e150));

produces 1e150.

If we add an explicit check for NAN in our max reduction (edit: better yet, in our max() function), that supposedly will slow down the common case where NAN is not encountered?

[npadmana]

@vasslitvinov Agreed -- this is more subtle/complicated than I thought. It turns out that IEEE 754-2008 defines max(NAN, num)==max(num,NAN)==num. And it looks like the C99 standard follows that.

As far as I can tell, Fortran depends on which compiler you use...

Chapel seems to hedge its bets :-)

writeln(max(NAN, 10.0));
writeln(max(10.0, NAN));

produces

10.0
nan

so the order matters. Which presumably is because it is using > to implement max?

Here is a fun thread on the Julia mailing list -- they started out with the standard, and then went to NaN propagation.
https://github.com/JuliaLang/julia/issues/7866

I tend to think of NaNs as bad, so I like NaN propagation. But I'm not sure what the right answer is, or even if there is one.

[vasslitvinov]

Interesting, the C standard treats NaN as missing data(?!!).

They discuss having two different versions of min/max, the default one propagating NaN.

Note that both options require mods in our current min/max(real,real) to make the result consistent.

P.S. I personally favor NaN propagation by default.

[mppf]

I think it's reasonable to expect min and max to propagate NaN. The other operations do (+ * etc). It's just an oddity of the IEEE standard that comparison with NaN always returns false, which combines with min/max being implemented with comparison, to result in it not being propagated. It seems like something we could easily fix in the module code defining min and max.

[bradcray]

I also agree with making Chapel min and max propagate NaN. @vasslitvinov: Would this be an easy change to make?

[vasslitvinov]

Looking how Julia did it in: https://github.com/JuliaLang/julia/pull/23155

For their min/max, they switched from using < to their already-existing isless() which is:

isless(x::AbstractFloat, y::AbstractFloat) = (!isnan(x) & isnan(y)) | signless(x, y) | (x < y)
signless(x, y) = signbit(x)::Bool & !signbit(y)::Bool

I wonder whether using | instead of || (if Julia has that) allows the compiler to make the three comparisons execute concurrently on the CPU.

The change here is simpler: https://github.com/apache/incubator-mxnet/pull/14262

Given the above, I'd change our max to be:

  inline proc max(x: real(?w), y: real(w)) return if (x > y) | isnan(x) then x else y;

@npadmana what do you think?

[vasslitvinov]

How should we adjust max for imag ? Currently it is:

  inline proc max(x: imag(?w), y: imag(w)) return if x > y then x else y;
// where...
  inline proc >(a: imag(?w), b: imag(w)) return __primitive(">", a, b);

[npadmana]

I would argue that we should not implement comparisons on imaginary numbers (just as we don't implement -- I hope -- comparisons on complex numbers). Imaginary numbers are much weaker that real numbers; for instance they aren't closed on multiplication.

Of course, a user is free to convert to a real and compare if that is the appropriate thing to do.

FWIW, Python and Mathematica both fail on these; I don't believe C has a pure imaginary type.

[bradcray]

I agree that we should retire (deprecate then remove) comparison on imaginaries (and am almost certainly the one at fault for thinking that it was reasonable to do so). There's discussion of why this doesn't make sense here: https://www.reddit.com/r/math/comments/5vpgzt/ordering_of_imaginary_numbers/

We don't support comparisons on complex numbers.

[bradcray]

@vasslitvinov: Just checking... are you running with the notions being discussed here, or just discussing them? If the former would you assign the issue to yourself?

[vasslitvinov]

Yes, adding isnan to min/max on reals and deprecating inequalities on imag.