[mycomputercrashsalot]

• ☑️ I understand it is strictly prohibited to use AI to write PR description
• AI usage disclosure: I used AI to help code and diagnose the bug with my guidance

Found while investigating #4614.

On CPU, exp, sin, cos and erf with float64 inputs went through the
float32 SIMD paths in mlx/backend/cpu/simd/math.h. The input was truncated
to float32, giving about 1e-8 error, and the float32 range clamps applied:
exp(-94.8) returned 0.0 instead of 6.7e-42.

Doubles now take the scalar std:: overloads per lane. The float32 and
complex paths are unchanged. erfinv still computes in float32 for double
inputs since there is no portable std::erfinv, that can be a follow-up.

test_unary_ops_keep_double_precision in test_double.py fails before the
change and passes after. test_double.py and test_ops.py pass (187 tests,
597 subtests).

Cost: float64 exp on CPU goes from 0.50 to 4.84 ms per 1M elements,
float32 is unchanged at 0.34 ms.