This PR adds NEP50 logic (https://numpy.org/neps/nep-0050-scalar-promotion.html) to masked arrays, which was overlooked in all binary operations.

For in-place operations specifically, this resulted in errors, see #27679. But in general, wrong dtypes were returned for all binary operations.

Tagging @mattip @seberg as related to the above issue and NEP50 that you are familiar with.

Resolves #27679

Short example of the bug

For an ndarray, this passes:

import numpy as np
res = np.array([1], dtype=np.uint8) + 2
assert res.dtype == np.uint8

The masked array equivalent fails:

res = np.ma.masked_array([1], dtype=np.uint8) + 2
assert res.dtype == np.uint8

Logic of the fix

Binary ops not in-place

We modify getdata() with a return_pyscalar argument to optionally return the scalar itself (instead of always converting to an array), and update the getdata() calls in _MaskedBinaryOperation and _DomainedBinaryOperation to use the scalar directly in their call to the corresponding np.umath function.

This fixes the issue for all not in-place binary ops.

Binary ops in-place

The problem:
For in-place operations, the problem is more delicate. As masked arrays vow to preserve the .data under the mask during an operation (see first lines here: https://numpy.org/doc/stable/reference/maskedarray.generic.html#operations-on-masked-arrays), they need a specific in-place logic.
In case of binary ops, the resulting .data is preserved under the final mask with respect to the left term, so self.data for in-place. To perform this without adding more memory usage, the existing logic in the ma replaces data under the future mask in other_data with the invariant value for the operation (0 for add/subtract, 1 for power/division), before passing it to the corresponding np.umath.
Therefore, we cannot rely on the np.umath function to take care of the NEP50 scalar logic, as the scalar other_data has to be converted to an array to substitute some cells to these invariant values before calling the np.umath.

Proposed fix:
We still use the return_pyscalar argument of getdata(), as for not in-place above, to get other_data as a scalar. Then, in the case where it is a scalar, we run an equivalent of result_dtype() on the scalar+array for this operation. However, result_dtype itself is not sufficient here. We need a function that can take as input (scalar, array_type, operator) and can both return the proper result_dtype for NEP50, and flag if we need to raise an OverflowError.

Does such a function exist internally somewhere in NumPy? (I couldn't find any)

For now, as a placeholder operation simulating the same behaviour, I've used a dubious call to a single value numpy array + the scalar for this operation.
For example for __iadd__:

if isinstance(other_data, (bool, int, float)):
    placehold_arr = np.array([0], dtype=self._data.dtype)  # Placeholder data
    placehold_arr.__iadd__(other_data)   # Raises appropriate errors
    other_dtype = placehold_arr.dtype   # Returns appropriate dtype to create a other_data as array

This fixes the issue for all in-place binary ops.

New tests

First, I've mirrored all the relevant tests in https://github.com/numpy/numpy/blob/main/numpy/_core/tests/test_nep50_promotions.py for masked arrays.

However, as those don't cover in-place operations, I've also added a test_nep50_inplace_ma that checks the behaviour for all combinations of in-place ops, scalars and dtypes.