[erwinzhang7]

Fix precision loss in CPU quantized matmul by accumulating in float32 instead of bfloat16 or float16.

The scalar _qmm_t path accumulates in T sum, while _qmm accumulates directly into the low-precision output buffer. This repeatedly rounds intermediate sums, causing substantial errors for large reduction dimensions.

The fix:

• Uses a float accumulator in _qmm_t.
• Uses a per-row float accumulation buffer in _qmm, converting to the output dtype once.
• Adds a deterministic regression test covering 84 combinations of dtype, layout, group size, and quantization bits.

Dequantization arithmetic and existing SIMD paths are unchanged.

Validation

Tested on M5 Max (18-core), macOS 27.0.1, CPU-only MLX builds.

• Regression test: 78/84 subtests fail on main; all pass with the fix.
• Python suite: 59 tests pass (16 skipped).
• C++ suite: 276/276 cases, 4433/4433 assertions pass.

For a 4096×4096 quantized matmul (4-bit, bf16):

Metric	Main	Fixed
Max absolute error (transposed)	0.682	0.0233
Max absolute error (non-transposed)	0.613	0.0261
Transposed latency (M=1)	72.9 ms	26.9 ms
Non-transposed latency (M=1)	14.9 ms	11.6 ms

The fix reduces maximum error by approximately 29× and improves transposed bf16 performance by 2.7×.

Remaining bf16 discrepancies come from low-precision dequantization and multiplication, which are intentionally outside this change.

• ☑️ I understand it is strictly prohibited to use AI to write PR description
• AI usage disclosure: Claude Code reviewed and worked with me to improve my fix

[zcbenz]

Instead of allocating a temporary memory with std::vector directly, can you allocate memory with MLX's allocator and register it with encoder's add_temporary and set_input_array like what GPU backends would do with intermediate buffers? It will make use of our memory cache, and will be required when we parallelize CPU scheduler.

[erwinzhang7]

Yessir done.