MAINT/BUG: improve gradient dtype handling - #9411

Let's give it a shot, thanks Eric.

Has there been an issue regarding float16? Be interesting for float16 applications if that case matters.

Has there been an issue regarding float16?

You mean here on github? I don't think so, it would be hard to spot - previously the result of gradient(some_float16) would just be promoted to float64, which is a higher precision than gradient(some_float32) gives!

The reason I asked is that float16 is not a type commonly used for computation. That has changed somewhat in recent years, where it is now used for some types of machine learning when range rather than precision is important. That has led to some requests for computation done in float16 rather than higher precision types.

Thinking it over, this change probably needs a compatibility note.

Note that the actually type used for computation has not changed - the old code was equivalent to doing a .astype(np.float64) right at the end. I suppose that still changes the behaviour of code doing np.gradient(some_f16) + some_other_f16