Wrong results of ma.isin() and ma.in1d() for masked arrays

For masked arrays, isin() and in1d() produce wrong results. In particular, the result is masked at wrong places. The issue was discussed initially on stackoverflow and it seems that the involved sorting does not handle the masked value correctly.

Reproducing code example:

import numpy.ma as ma
a = ma.MaskedArray([[1,2,3],[4,5,6]], [[True,False,False],[False,False,False]])
ta = ma.array([1,4,5])
ma.isin(a, ta)

The output of ma.isin() is wrong (see below). I use the following workaround:

import numpy as np
def ma_isin(array, comparison):
    return ma.MaskedArray(
        data=np.isin(array, comparison),
        mask=array.mask.copy())

Output and expected results:

These are the results from the code above

>>> a
masked_array(
  data=[[--, 2, 3],
        [4, 5, 6]],
  mask=[[ True, False, False],
        [False, False, False]],
  fill_value=999999)
>>> ta
masked_array(data=[1, 4, 5],
             mask=False,
       fill_value=999999)
>>> ma.isin(a, ta)
masked_array(
  data=[[False, False, False],
        [True, True, --]],
  mask=[[False, False, False],
        [False, False,  True]],
  fill_value=True)

The issue is mainly that the last element in the result is masked, and probably also that the first element is not masked. I would expect the following output:

>>> ma_isin(a, ta)
masked_array(
  data=[[--, False, False],
        [True, True, False]],
  mask=[[ True, False, False],
        [False, False, False]],
  fill_value=True)

in1d() shows an analogue issue.

NumPy/Python version information:

1.21.2 3.9.4 (default, Apr 9 2021, 16:34:09)
[GCC 7.3.0]

Start by running the masked-array reproducer against ma.isin() and the analogous ma.in1d() case, then trace their masked-array sorting and comparison entry points. Verify that masked input positions remain masked in the result and that unmasked values receive the correct membership results; add regression coverage for both functions.