fix(ailego): accumulate FP16 dense distances in FP32 on AVX-512 FP16 CPUs

[JasonXuDeveloper]

Problem

Closes #786.

The ailego AVX512-FP16 kernels for FP16 inner product and squared L2 (single-pair and batch) accumulate in __m512h, so every partial sum is rounded to FP16. On 1024-d unit vectors this adds ~3e-4 of error to each score. On Sapphire Rapids, FP16 COSINE/IP recall@10 falls from ~0.997 to ~0.965, and the same index returns different neighbours on CPUs without AVX512-FP16. The numbers are in the issue.

Change

• Remove the AVX512_FP16 branches from the FP16 dense dispatchers (InnerProductMatrix, MinusInnerProductMatrix, SquaredEuclideanDistanceMatrix, InnerProductDistanceBatchImpl, SquaredEuclideanDistanceBatchImpl). These calls now fall through to the existing AVX-512F kernels, which widen with vcvtph2ps and accumulate in FP32.
• Delete the FP16-accumulating kernels and the helpers that only they used (FMA_FP16_AVX512FP16, HorizontalAdd/Max_FP16_V512).

Rewriting the AVX512-FP16 kernels to use _mm512_cvtph_ps + _mm512_fmadd_ps would duplicate the AVX-512F kernels, so this PR routes to those instead. This matches turbo, which leaves its own AVX512-FP16 kernels unregistered for the same reason (turbo.cc:359-362).

Not touched:

• The sparse InnerProductSparseInSegmentFp16AVX512FP16.
• The ARM FEAT_FP16 branch of ACCUM_FP16_1X1_NEON, which the default build doesn't compile because the TUs get -march=armv8-a.
• turbo's neon_fp16 cosine.

These could be follow-ups if you want them.

Performance

Not benchmarked: I don't have an AVX512-FP16 machine. What to expect on one:

• Per 32 elements per vector, the AVX-512F loop issues one load, one vextracti64x4, two vcvtph2ps and two FP32 FMAs. The removed loop issued one load and one vfmadd231ph. When the data is cache-resident and the loop is compute-bound, this will be slower.
• In HNSW search at large dimensions (a 1024-d FP16 vector is 32 cache lines, fetched from random locations), memory access should dominate, so the end-to-end cost should be small. turbo.cc:359-362 says the same about turbo's FP16 kernels: they only helped for cache-resident dims <= 256.
• CPUs without AVX512-FP16 (Zen 4, Ice Lake, Skylake-SP) already run exactly this path, so it isn't new code.

Existing coverage is thin. Only the 1x1 baseline timing in DISABLED_InnerProduct_Benchmark (inner_product_matrix_fp16_test, dims 64 and 256) goes through the changed single-pair dispatch, and nothing benchmarks the batch kernel or large dims. The most useful number would be HNSW search latency on a 1024-d VECTOR_FP16 COSINE collection on Sapphire Rapids, before and after this change. If the small-dimension cost matters, a PR comment below proposes an alternative.

Test

DistanceMatrix.InnerProduct_Fp32Accumulation compares single-pair and batch (12 + tail) FP16 inner products on skewed 1024-d unit vectors against a double-precision reference over the same FP16 values, with a bound of 2e-5. FP32 accumulation is off by less than 1e-6. FP16 accumulation is off by about 1.6e-4.

• macOS arm64: inner_product_matrix_fp16_test, euclidean_distance_matrix_fp16_test and cosine_distance_matrix_fp16_test pass.
• Negative control: building the ailego NEON TUs with -march=armv8.2-a+fp16, which turns on the FP16-accumulating NEON path, makes the new test fail (single-pair error 2.0e-4, batch 6.9e-5).
• The changed x86 TUs compile cleanly with -march=sapphirerapids. I don't have AVX512-FP16 hardware to run the test on x86. On such a machine it should fail before this change and pass after it.

🤖 Generated with Claude Code

[JasonXuDeveloper]

A question for reviewers before this goes further: do you want to keep a native FP16 FMA path on AVX512-FP16 CPUs?

This PR takes the simplest fix and routes these calls to the existing AVX-512F kernels (convert, then FP32 FMA). That fixes precision, but in compute-bound, cache-resident cases it will be slower than the removed vfmadd231ph loop (see "Performance" in the description; not benchmarked). If that matters to you, here is an alternative that keeps most of the FP16 FMA throughput.

Proposal: FP16 FMA with periodic widening to FP32

Accumulate in __m512h only for a short block of kFlushSteps iterations, then convert the partial sums to FP32 and add them into FP32 accumulators:

template <size_t kFlushSteps>
float InnerProductFp16Flush(const Float16 *lhs, const Float16 *rhs, size_t size) {
  __m512 acc_lo = _mm512_setzero_ps(), acc_hi = _mm512_setzero_ps();
  size_t i = 0;
  while (i + 32 <= size) {
    __m512h partial = _mm512_setzero_ph();
    for (size_t s = 0; s < kFlushSteps && i + 32 <= size; ++s, i += 32) {
      partial = _mm512_fmadd_ph(_mm512_loadu_ph(lhs + i),
                                _mm512_loadu_ph(rhs + i), partial);
    }
    __m512i bits = _mm512_castph_si512(partial);
    acc_lo = _mm512_add_ps(acc_lo, _mm512_cvtph_ps(_mm512_castsi512_si256(bits)));
    acc_hi = _mm512_add_ps(acc_hi, _mm512_cvtph_ps(_mm512_extracti64x4_epi64(bits, 1)));
  }
  float sum = _mm512_reduce_add_ps(_mm512_add_ps(acc_lo, acc_hi));
  // masked tail as in the current kernels
  return sum;
}

• Cost: with kFlushSteps = 4, each 128 elements take 4 vfmadd231ph + 1 vextracti64x4 + 2 vcvtph2ps + 2 vaddps. The FP32 path in this PR takes 4 extracts + 8 converts + 8 FP32 FMAs for the same 128 elements. The batch kernel would do the same per row, sharing the query loads as it does now.
• Precision: each FP16 partial sum now covers at most kFlushSteps products instead of the whole vector, so rounding stays near the magnitude of a few products rather than of the running total. It is still worse than FP32 accumulation, because every vfmadd231ph (including the first in a block) rounds its result to FP16. So kFlushSteps is a precision/speed knob, and the precision test in this PR would need to set its bound for it.
• Verified so far: the sketch above compiles for -march=sapphirerapids and emits the expected instruction mix. It has not been run or benchmarked, since I don't have AVX512-FP16 hardware.

If you prefer this direction, I can rework the PR to use it for the IP (single and batch) and squared L2 kernels. I would pick kFlushSteps by measuring recall on the repro from #786 and bounding the error in the test. That still needs someone with a Sapphire Rapids machine to confirm it is actually faster than the FP32 path. Otherwise, I'd keep this PR as is.