Compute float64 exp, sin, cos and erf in double

[AKnassa]

Proposed changes

Fixes #4158.

exp, sin, cos and erf in mlx/backend/cpu/simd/math.h are hand-rolled
float32 approximations, and they convert their input to Simd<float, N>
unconditionally. float64 therefore gets float32 accuracy:

max rel. error, float64	before	after
exp	3.00e-07	0.00e+00
sin	7.75e-07	0.00e+00
cos	3.42e-04	0.00e+00
erf	9.44e-04	0.00e+00

exp is worse than imprecise, because it also inherited float's saturation
points, so a float64 array returns the wrong value outright:

exp(100.0)  -> inf   (should be 2.688e+43)
exp(709.0)  -> inf   (should be 8.218e+307)
exp(-200.0) -> 0     (should be 1.384e-87)

This routes double through libm and leaves float32 untouched.

Worth flagging: the existing test_unary_ops in python/tests/test_double.py
asserts allclose(y, y_double.astype(mx.float32, mx.cpu)) — it compares the
float64 result against the float32 one, so it could not catch any of this.

Performance

This is a real cost. 2^20 elements, M2 Pro, best of 7 runs:

	float32 (unchanged)	float64 before	float64 after
exp	0.437 -> 0.438 ms	0.573 ms	5.908 ms
sin	0.804 -> 0.807 ms	1.067 ms	6.659 ms
cos	0.786 -> 0.791 ms	1.032 ms	6.743 ms
erf	0.877 -> 0.878 ms	1.160 ms	11.641 ms

float32 is unchanged, as expected from a compile-time branch. The old float64
speed was the speed of computing the wrong answer.

I checked whether that speed can be recovered. Accelerate's vectorized double
math (simd::exp / simd::erf on double4) measured no faster than scalar
libm on this machine — 2.456 vs 2.558 ms for exp, 9.653 vs 9.951 ms for erf,
with bit-identical results — so the cost looks intrinsic to evaluating these at
double precision rather than a poor choice of implementation. Happy to take
direction if you would prefer a different trade-off here.

erfinv still narrows to float. It has no libm equivalent and would need a
double precision polynomial, so I left it for a separate change rather than
half-do it.

Testing

Two new tests in python/tests/test_double.py: one pins the four ops to double
accuracy against numpy, one pins the range of exp. Both fail before this
change and pass after.

• Full Python suite: 821 tests, no new failures.
• Full C++ suite: 247 cases, 3326 assertions, all pass.

Built CPU-only (MLX_BUILD_METAL=OFF). float64 is CPU-only in MLX, so that is
where this behaviour lives, but note the GPU suite was not run on my machine.

Checklist

• I have read the CONTRIBUTING document
• I have run pre-commit run --all-files to format my code / installed pre-commit prior to committing changes
• I have added tests that prove my fix is effective or that my feature works
• I have updated the necessary documentation (if needed) — no documentation change needed, there is no API change

[zcbenz]

At the moment we don't plan to taking non-critical simd changes until we merge #3019.