BUG: masked array division broken with np.seterr(under="raise")

[tacaswell]

Describe the issue:

When np.seterr(under="raise") is set, division does not work with masked arrays due to an internal multiplication underflow.

Reproduce the code example:

np.seterr(under="raise")
x=np.arange(0, 3, 0.1)
X = np.ma.array(x)
x2 = x / 2.0  # <- works
X2 = X / 2.0  # <- fails

Error message:

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:4275, in MaskedArray.__truediv__(self, other)
   4273 if self._delegate_binop(other):
   4274     return NotImplemented
-> 4275 return true_divide(self, other)

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:1171, in _DomainedBinaryOperation.__call__(self, a, b, *args, **kwargs)
   1169 domain = ufunc_domain.get(self.f, None)
   1170 if domain is not None:
-> 1171     m |= domain(da, db)
   1172 # Take care of the scalar case first
   1173 if not m.ndim:

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:858, in _DomainSafeDivide.__call__(self, a, b)
    856 a, b = np.asarray(a), np.asarray(b)
    857 with np.errstate(invalid='ignore'):
--> 858     return umath.absolute(a) * self.tolerance >= umath.absolute(b)

FloatingPointError: underflow encountered in multiply

Python and NumPy Versions:

Python 3.11, numpy 1.26.3

Runtime Environment:

[{'numpy_version': '1.26.3',
  'python': '3.11.6 (main, Nov 14 2023, 09:36:21) [GCC 13.2.1 20230801]',
  'uname': uname_result(system='Linux', node='belanna', release='6.7.4-arch1-1', version='#1 SMP PREEMPT_DYNAMIC Mon, 05 Feb 2024 22:07:49 +0000', machine='x86_64')},
 {'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
                      'found': ['SSSE3',
                                'SSE41',
                                'POPCNT',
                                'SSE42',
                                'AVX',
                                'F16C',
                                'FMA3',
                                'AVX2'],
                      'not_found': ['AVX512F',
                                    'AVX512CD',
                                    'AVX512_KNL',
                                    'AVX512_KNM',
                                    'AVX512_SKX',
                                    'AVX512_CLX',
                                    'AVX512_CNL',
                                    'AVX512_ICL',
                                    'AVX512_SPR']}}]

Context for the issue:

This was reported against Matplotlib in matplotlib/matplotlib#27770

@2sn will have to provide any additional context.

[tacaswell]

This can be worked around by using X2 = X*0.5 instead.

[seberg]

This domain stuff is fatally flawed, but I guess np.errstate(all="ignore") makes sense in any domain related calculation especially.

[carlosilva10260]

I'm working on this issue