Below a story, that nobody needs to read, it is a clear fix. I am looking the integer cast difficulties, but it is still a bit difficult to unravel where to best do something about it.

(My preference would be to never store an incompatible fill_value basically, which is likely easy enough, but I am not sure yet that fixes everything. Plus the easiest way to do that may be to ensure that default_fill_value() returns something very different. That seems good, but maybe not for backporting.)

The unary comparison "domains" can return masked (the others can't so there is a dissonance here that doesn't really add up).

This adds a force-cast to boolean, since the domain result only makes sense as boolean, but the masked scalar is float64...

Now, why did this change in NumPy 2? It didn't... Both before and after a TypeError was raised (incorrectly) and the context path was then skipped completely (the ufunc machinery has a try/except).

What changed is that I removed __array_prepare__ because it had basically no use. Turns out, that in this error path (and only there) it was important because the "context" path was not taken, but __array_prepare__ called __array_finalize__ so it also applied the mask.
For unary ufuncs, it thus neatly (but very round-about) fills the bad logic here.
(For non-unary ufuncs, the "domain" functions don't return masked arrays so the path cannot happen)

A revier with a keen eye may notice that the old code returned a masked array rather than the masked constant here. Note, however, that this is special the incorrectly taken no-context path, because the other path doesn't have the "return scalar" logic (which is probably also a bug, but who is counting).

Closes gh-25635