BUG: overflow encountered in multiply when multiplying masked arrays of float32

Describe the issue:

When doing _MaskedBinaryOperation on two masked arrays the actual operation is done on underlying data arrays ignoring mask. The relevant code is https://github.com/numpy/numpy/blob/main/numpy/ma/core.py#L1008-L1013

        # Get the data, as ndarray
        (da, db) = (getdata(a), getdata(b))
        # Get the result
        with np.errstate():
            np.seterr(divide='ignore', invalid='ignore')
            result = self.f(da, db, *args, **kwargs)

Division by zero and invalid operations are ignored for this operation, but not over- or underflow. Also this code silently swallows errors occurred for nonmasked elements.

Reproduce the code example:

import numpy as np

np.seterr(all='raise')

x = np.array([[1, 1e20], [1e20, 2]], dtype=np.float32)
y = np.array([[3, 1e20], [1e20, 4]], dtype=np.float32)

undef = np.float32(1e20)

z = np.ma.masked_equal(x, undef)
w = np.ma.masked_equal(y, undef)

z*w

Error message:

Traceback (most recent call last):
  File "test.py", line 13, in <module>
    z*w
  File "<snip>/miniconda3/lib/python3.7/site-packages/numpy/ma/core.py", line 4169, in __mul__
    return multiply(self, other)
  File "<snip>/miniconda3/lib/python3.7/site-packages/numpy/ma/core.py", line 1021, in __call__
    result = self.f(da, db, *args, **kwargs)
FloatingPointError: overflow encountered in multiply

NumPy/Python version information:

1.19.2 3.7.5 (default, Oct 25 2019, 15:51:11)
[GCC 7.3.0]

Start with the masked-array multiplication path in numpy/ma/core.py around lines 1008-1013, then run the reproduction against the current NumPy behavior. Check how masked and unmasked elements are handled under np.seterr(all='raise'). Done means masked overflow does not raise while errors from unmasked elements are not silently suppressed, with regression coverage for the reported example.