BUG: Add promoter to `ldexp` for python integers to prevent overflow

[MaanasArora]

Fixes #29182.

Adds a promoter for first argument PyArray_PyLongDType to ldexp. It casts straight to double, not sure if that is more arbitrary than needed.

Thanks for reviewing!

[mattip]

What will happen with complex input?

[MaanasArora]

@mattip this only promotes integer inputs (Python scalar integer, actually), so the rest would stay the same?

[seberg]

Thanks a for looking at this, this is very cool to see this use!

I am curious if we can push this further and actually use the promoter always. That is, we:

• Set the first and last dtype to PyArray_CommonDType(first_dtype, AbstractFloatDType):
  • If the last was forced, we just set the first to that instead.
  • Yes, we ignore the second (integer) argument here considering how it is currently defined.
• Set the second dtype to PyArray_CommonDType(second_dtype, int_dtype). That just seems like it always works in practice, so I don't really feel a need to do something more complicated.

[MaanasArora]

Thanks, that extension makes a lot of sense!

I tried implementing this, but it causes a segmentation fault. It seems to be due to the fact that common_dtype is not defined on FloatPyArray_FloatAbstractDType (in the dtypemeta declaration). I'm wondering if the abstract classes support the common dtype operation yet? (I assume the python integers won't be the ones casting up to floats)

Edit: that does seem to be it; reverting the operation for the first dtype and setting the others in this way doesn't crash (though of course the promotion doesn't work properly)

[seberg]

Oh, hmmm, the pyfloat one should work ok probably, although the abstract version woukd be nicer.

[MaanasArora]

PyFloat doesn't seem to work well oddly, but Float does. I've pushed that, but yeah it's not as nice. It could be interesting to look into supporting common_dtype more fully for the abstract types if you'd like!

[seberg]

Making the abstract and PyFloat do the same thing probably makes sense, but the fact that it doesn't work here seems like a bug somewhere that I (or we) need track down to make this work (and a couple of other cases, because we should use this for other float only functions).

[MaanasArora]

To be clear it doesn't crash, it just doesn't promote further up than float16 so causes overflow.

In line 223-226 of abstractdtypes.c we do in float_common_dtype:

    else if (other == &PyArray_PyLongDType) {
        Py_INCREF(cls);
        return cls;
    }

This is equivalent to doing no promotion as we already know it's a float, but looking at NEP 50, this actually seems correct behavior - Python ints can promote to the lowest float (given we need a float) right? So for ldexp if we don't consider argument values, we would have to just use something that can fit all values?

[MaanasArora]

Yes, we promote to DoubleDType when the other dtype is any legacy integer, but do not promote (keep the same class) when it is PyLong. That seems a bit inconsistent, so it is probably buggy in at least one case? Unless I'm missing something

https://github.com/numpy/numpy/blob/583993f2f9f7cafe265b3a5505505fb8d2d376e7/numpy/_core/src/multiarray/abstractdtypes.c#L217-L226

[seberg]

Ah, that makes sense (if I am getting the right picture). We are getting the abstract version as a result. In other paths that is OK because we only need a concrete dtype instance (via default_descr), but here we need to convert the abstract one to the concrete Float64 one. I actually don't think I needed that before, it seems like a missing piece for generic promotion! We could do it via the default_descr type, which not be so bad in practice, although it's reversed...

[MaanasArora]

Ah, I see, thanks! Yes, I think that's right, we're getting the PyFloat back, which is abstract. Should common_dtype not return abstract dtypes (i.e., are dtypes concretized there)? Then yes, we could do default_descr if the dtype is abstract. I'm not sure it being reversed is a concern because we often rely on symmetry in these casts?

[seberg]

Yeah, I suppose it probably isn't a concern in practice, so the choice we have may be the pragmatic one, even if it is slightly weird/reversed if we use it here. Reverse or not, until we have a case where this fails, it may be a reasonable to path to "promote to a concrete DType".

[MaanasArora]

Makes sense! If it's a missing step, I suppose we should add it to PyArray_CommonDType unless that's too big and we'd rather have it localized? Just pushed it there! I did notice though, the python types (like PyFloat) don't seem to have the abstract flag, while tp_base is set to the abstract version. I'm just checking if default_descr exists but not sure if that's ideal.

(I realize this is an unrelated change to this PR, so completely happy to revert and open another PR if you prefer.)

[seberg]

I suppose we should add it to PyArray_CommonDType unless that's too big and we'd rather have it localized?

Yeah, I am not sure I want that function to ensure no abstract result. We could add a mini internal helper for now that ensures we have a concrete dtype, also means we can refactor it in a single place if we ever have to. We should check NPY_DT_is_abstract before doing this dance, there is a NPY_DT_CALL_default_descr macro and we ensure that it cannot be NULL (although it could fail, but that is probably OK).

(Also could check if you need the PyFloat rather than the abstract float version, since that didn't actually help initially.)

[MaanasArora]

Makes sense! I actually meant that NPY_DT_is_abstract didn't work for some reason. It doesn't seem we initiate the PyFloat with NPY_DT_ABSTRACT, though we do it for AbstractDType. If we want it to be abstract should we add the flag? (Or is it inherited from tp_base? It doesn't seem to be)

Using the AbstractDType would be a better call probably, but since common_dtype is undefined, doing that still causes a segmentation fault. (I can define the common dtype if you like, but should probably do it for all abstract dtypes then.)

[seberg]

It doesn't seem we initiate the PyFloat with NPY_DT_ABSTRACT, though we do it for AbstractDType.

Hmmmmmmm, we had some reason for changing that. FWIW, it's probably fine then, we should be using the abstract version. The common_dtype for the abstract version should be defined identically, but not sure it would be easy to not just duplicate the code (which is probably fine).

[MaanasArora]

We can just pass the float_common_dtype pointer right? There are no dt slots defined for the abstract dtypes as far as I can tell, but we could create that struct and then pass common_dtype = float_common_dtype?

Edit: unless there is a reason to not define the slots (I guess not, since the only other two are default_descr and discover_descr_from_pyobject, which would not be needed, so coincidence?)

[seberg]

We can just pass the float_common_dtype pointer right?

Basically except that may return one of the Py*DType which would be wrong when we need to return one of the Abstract*DType.

[MaanasArora]

Hmm, there doesn't seem to be a reference to the Py*DType at least in the float version! So it should probably be fine.

I've added two slots to the abstract float dtype (I can add if you want them on the other abstract dtypes as well). The internal helper, should it be only for promoters?

Edit: sorry, I suppose if it's a helper it can just be used directly in ldexp? Pushed this.

[MaanasArora]

The new test seems to reliably fail on armhf (will investigate).

[jorenham]

A rebase wouldn't hurt

[seberg]

Pushed a small fixup mostly for additional error checks. I also broadened it up to all integer inputs (including int32, etc.).

NumPy promotion should go to float64 for integers of any type by default for functions that do not support integers normally, I think. Since we add the machinery to do that, might as well do that now I think.

Let's see if tests pass now, although I am not sure what that armhf failure was... If they do, I think I'll merge it soon.

[seberg]

What will happen with complex input?

Nothing, the promoter will never be called and even if it was, it wouldn't do anything. There are no complex loops anyway.

[MaanasArora]

Thanks for the cleanup!

I somewhat suspect it may be an issue with the SIMD code, a way forward may to ignore this warning/FPE and check what the result type (and value, but if that fails just the type to see if the promotion is right).

Yes, likely some downstream bug, since it's only on a few platforms. I've added a test to check the result dtype is np.float64 and suppress the warnings (not sure if the TODO is too much, but I guess good to be explicit there).

[seberg]

Yes, likely some downstream bug

Well, this is should be in the SIMD code path I think, but let's see. So it may or may not be the system math being at fault for spuriously giving a warning.