[tetsuo]

Metal's addmm can lose precision prematurely during the alpha * (A @ B) + beta * C computation.

Currently, the GEMV and GEMVT kernels cast the accumulator and coefficients down to the output data type before applying scaling and bias. Similarly, the split-K output transform computes the scaled bias using output precision. For BF16 and FP16, these early conversions can cause incorrect rounding or avoidable overflow.

For example:

import mlx.core as mx

a = mx.array([[1, 1 / 256]], dtype=mx.bfloat16)
b = mx.array([[1], [1]], dtype=mx.bfloat16)
c = mx.array([1 / 256], dtype=mx.bfloat16)

print(mx.addmm(c, a, b, stream=mx.gpu).item())

The current implementation returns 1.0. The corrected implementation returns the exactly representable result, 1.0078125.

This change:

• Keeps scaling and bias addition in accumulator precision, then converts the final result to the output dtype.
• Corrects the shared affine transforms. Regular dense GEMM already uses accumulator-typed transforms.
• Adds gemv.h to the precompiled kernel's build dependencies.
• Adds regression tests covering BF16, FP16, and F32, including batching, transposed inputs, strided bias, and overflow cases.
• Adds a focused addmm benchmark.

Tested and validated on Apple M4 Pro.