I would suggest just writing:

x = pi * x

Maybe we can even continue with (or where if you like):

y = sin(x)
x[x == 0] = 1  # or where, or copyto
return y / x

That way, we only guess at 1, which seems like a smaller guess? Although, I didn't quite see through whether spelling it may be slower/require more memory.

I was briefly considering using np.nextafter(..., 1) for the filling value, which fixes the issue of custom dtypes like ml_dtypes.bfloat16.
(But, we should likely np.finfo(...).smallest_normal, even without custom dtypes, it is currently slow and actually currently not defined for all longdoubles; although, both should probably be changed.)

Now, it turns out that bfloat16 currently returns float64 in the where (I am not quite sure why, I thought I fixed that!?).
But, I think it still fixes the problem, although it will change the result for ml_dtypes.bfloat64 to be float32 rather than float64 I believe (which may be a good thing, though).

Sorry, that was silly of course, we would have to copy things twice...

I'll have to think about it, maybe we can just fix it up later, although smallest_normal has some issues that make it very slow to use right now (and it may indeed fail).
finfo().eps works for all NumPy dtypes it seems, so maybe that is the better choice (if we use it for NumPy dtypes only anyway).

@seberg, sorry, I'm not sure what the conclusion is here as to how this PR should change - do you have a concrete suggestion?

I would lean towards fixing the multiplication to remove the isdtype. Then use dtype.kind == "f" (assume that means finfo works) and go with finfo(dtype).eps because that seems to work as well, and (whether a thing that should be fixed or not) is currently vastly faster and won't have problems with longdouble on double-double architectures.

But, you need to sanity check that thought of course and maybe you have a better idea.

EDIT: If dtype != "f" you would just keep guessing at the old value for compatibility probably).