floats

README

floats

SIMD-accelerated reductions over []float32 and []float64 — the dot product, sum, min, max, Euclidean (L2) distance and cosine similarity that sit at the heart of vector search, embeddings and numerics, plus the float64 Axpy/Xpby/AxpyDot/Scale updates behind iterative solvers. Pure Go: CGO_ENABLED=0, stable Go, no GOEXPERIMENT. The hot reduction loops are real FMA SIMD on all six 64-bit Go targets, assembly generated by go-asmgen — on a plain go build.

import "github.com/go-simd/floats"

sim := floats.CosineSimilarity(query, doc)        // []float64 relevance score
d   := floats.Float32Distance(embA, embB)         // []float32 L2 distance
dp  := floats.Dot(a, b)                            // a·b

API

For []float64 (top-level) and []float32 (the Float32* functions):

float64	float32	meaning
Dot(a, b)	Float32Dot(a, b)	a·b = Σ aᵢbᵢ
Sum(a)	Float32Sum(a)	Σ aᵢ
Min(a) / Max(a)	Float32Min(a) / Float32Max(a)	min / max element (NaN-propagating)
Distance(a, b)	Float32Distance(a, b)	‖a−b‖₂ = √Σ(aᵢ−bᵢ)²
CosineSimilarity(a, b)	Float32CosineSimilarity(a, b)	a·b / (‖a‖₂‖b‖₂)

Dot, Distance and CosineSimilarity require len(a) == len(b) and panic otherwise; Min/Max panic on an empty slice; Sum of an empty slice is 0. The float32 dot/sum/squared-difference are accumulated in float32; cosine and distance take the final square root in float64 for accuracy.

Elementwise updates (BLAS level 1, float64)

The in-place updates that iterative solvers (CG, GMRES), Householder QR and LU elimination spend their time in:

function	computes
Axpy(alpha, x, y)	y[i] = math.FMA(alpha, x[i], y[i]) (y += αx)
Xpby(x, beta, y)	y[i] = math.FMA(beta, y[i], x[i]) (y = x + βy)
AxpyDot(alpha, x, y, w)	Axpy, then returns Σ y[i]*w[i] over the updated y, in one pass
Scale(alpha, x)	x[i] *= alpha

// one conjugate-gradient step
floats.Axpy(alpha, p, x)              // x += αp
rr := floats.AxpyDot(-alpha, ap, r, r) // r -= αAp, rr = r·r
floats.Xpby(r, rr/rrOld, p)            // p = r + βp

They operate on the first len(x) elements and panic if y (or w) is shorter. y may be the same slice as x or w; partially overlapping slices are not supported.

Bit-identical on every architecture. Each updated element is rounded once, exactly as math.FMA rounds it: VFMADD on amd64 (AVX2+FMA), VFMLA/FMADDD on arm64, the math.FMA loop everywhere else (one hardware FMA instruction on ppc64le, s390x, riscv64 and loong64). The updated slice never depends on the architecture or the CPU features detected at run time (NaN payloads aside). This is not the unfused y[i] + alpha*x[i] that amd64 Go code computes, so replacing such a loop can change the last bit; the fused result is the more accurate one. AxpyDot's returned sum is a reduction and follows the lane-order rules of Floating-point determinism.

arch	elementwise kernel
amd64	AVX2+FMA, ×4 unrolled (16 elements/iteration), gated on HasAVX2 && HasFMA
arm64	NEON, ×4 unrolled (8 elements/iteration, four-register VLD1/VST1)
ppc64le, s390x, riscv64, loong64	portable math.FMA loop: no kernel until one is measured on the real silicon, as for the reductions

Why this library

The Go SIMD-float landscape, compared:

	dot/dist/cosine	float32	amd64	arm64	ppc64le / s390x / riscv64 / loong64
gonum/floats	dot, L2 only	✗ (float64 only)	SSE	✗	✗
viterin/vek	✓	✓	AVX2+FMA	✗	✗
kelindar/simd	✗ (sum/min/max only)	✓	AVX2	NEON	✗
go-simd/floats	✓	✓	AVX2+FMA	NEON	✓ (all four)

This is the only pure-Go float reduction library that ships the full vector-search API (dot / Euclidean / cosine) for both float32 and float64 on all six of Go's 64-bit SIMD architectures.

How the SIMD is built

Every kernel is an FMA-based vectorised reduction: a multi-lane accumulator folded with a fused multiply-add at the vector stride, then a horizontal reduce and a scalar tail.

arch	ISA	FMA op	notes
amd64	AVX2 + FMA3	VFMADD231P{S,D}	runtime-gated on cpu.X86.HasAVX2 && HasFMA; scalar fallback otherwise
arm64	NEON	VFMLA	4 accumulators, ×4 unroll for ILP; baseline, no gate
riscv64	RVV	VFMUL+VFREDOSUMVS	length-agnostic VSETVLI stripmining; gated on cpu.RISCV64.HasV
s390x	vector facility (big-endian)	VFMADB	Sum/SumSqDiff vectorised; float32 uses the scalar reference (no .SB ops in the assembler); z13 baseline, no gate. float64 Dot has no VX kernel — on real z15 the VFMADB dot kernel measured 0.56× the gc-autovectorized scalar loop (11.6 vs 20.6 GB/s), so Dot routes to that loop instead (see below)
ppc64le	VSX (float32 only)	XVMADDASP	the XV* ops are WORD-encoded (not in the released assembler); gated on cpu.PPC64.IsPOWER9 (the SP ops are ISA-3.0). float64 has no VSX kernel — on real POWER9 the gc-autovectorized scalar loop beats it (see below), so float64 routes to that loop
loong64	LSX (float32 only)	vfmadd.s	vfmadd.s is WORD-encoded (no Go mnemonic yet); add/sub/xor use VADDF/VSUBF/VXORV. float64 routes to the gc-autovectorized loop. LA464 baseline, no gate

The .s files are committed; regenerate with go run *_gen.go (the generators are build-tagged //go:build ignore). All six are validated in CI — amd64/arm64 natively, the other four under QEMU (tonistiigi/binfmt) — at 100 % test coverage, with differential tests and FuzzDot/FuzzDistance/FuzzAxpy against the scalar reference on every architecture. The elementwise kernels are also checked bit for bit against math.FMA for every length 0..257 at four element offsets, and run on slices placed right against PROT_NONE guard pages so that a vector load past either end would fault.

Floating-point determinism

A SIMD reduction does not add its terms left-to-right: it keeps several independent accumulator lanes (for instruction-level parallelism) and folds them at the end, so its rounding differs from a naïve scalar loop. For dot / distance / cosine, bit-exactness against a naïve loop is neither meaningful nor required — correctness to a few ULP is.

So the reference the tests and fuzzers check against (dotLanes, sumSqDiffLanes, …) reduces in the same lane-blocked order the kernels use, and each kernel is required to match it to a tight, documented relative tolerance (1e-12 for float64, 1e-5 for float32) on every architecture. The differential tests and fuzzers enforce this, making the results deterministic and consistent across arches.

Performance

Apple M-series (arm64, native), go test -bench, ns/op (lower is better); simd is this library, naive is the textbook left-to-right Go loop (which the Go compiler already turns into a pipelined scalar FMADD loop):

op	n=8	n=64	n=512	n=4096
Dot float64 simd	7.6	14.3	71.7	539
Dot float64 naive	2.4	18.5	152	1075
Float32Dot simd	9.5	11.7	39.4	276

The SIMD dot product is ~2× the compiler's scalar loop on medium/large vectors (≈57 GB/s float64, ≈59 GB/s float32 at n=512), and loses only on tiny inputs (n≤8) where the unroll/setup overhead dominates — the regime where you would not reach for SIMD anyway.

Elementwise updates

Median ns/op of 3 runs (-cpu 1 -benchtime 200ms, lower is better). naive is the loop callers write by hand (y[i] += alpha*x[i], …), single core:

op	n	Apple M4 Max simd	naive	speedup	AMD EPYC 7763 (CI) simd	naive	speedup
Axpy	1e3	185	381	2.06×	121	634	5.25×
Axpy	1e5	22 992	40 087	1.74×	20 258	63 152	3.12×
Axpy	1e6	270 906	449 693	1.66×	201 142	631 576	3.14×
Xpby	1e3	184	368	2.01×	121	633	5.25×
Xpby	1e5	22 888	36 138	1.58×	20 373	63 207	3.10×
Xpby	1e6	234 768	387 107	1.65×	201 311	631 323	3.14×
AxpyDot	1e3	231	1 223	5.30×	160	933	5.85×
AxpyDot	1e5	33 085	160 507	4.85×	31 208	94 524	3.03×
AxpyDot	1e6	439 445	1 279 334	2.91×	301 463	946 170	3.14×
Scale	1e3	101	362	3.57×	83	633	7.67×
Scale	1e5	15 753	35 058	2.23×	11 967	62 810	5.25×
Scale	1e6	158 718	351 821	2.22×	108 396	629 221	5.80×

The M4 Max runs were on a shared workstation (1-minute load ≈ 6), the EPYC ones on a shared GitHub Actions runner (ubuntu-latest), so treat the absolute numbers as indicative. Every row is a win. The ubuntu-24.04-arm runner (Neoverse) measured 1.34–2.83× over the same benchmarks. ppc64le, s390x, riscv64 and loong64 run the scalar math.FMA loop and were not measured.

Honesty note on the other arches. amd64's AVX2+FMA kernel is validated for correctness on real x86 (an AVX2/FMA-capable VM) but the native-hardware throughput numbers are pending (the CI/dev x86 runner used here is TCG-emulated, so its timings are not representative). riscv64 and loong64 are QEMU-validated for correctness; native-hardware performance numbers are pending access to that hardware. s390x is now measured on real z15 (below).

ppc64le — measured on real POWER9 (GCC Compile Farm cfarm433, go1.26.4, 2026-06-27, Dot throughput MB/s, higher is better):

n	float64 (autovectorized)	float64 if it used VSX	float32 VSX	float32 naive
8	5700	2520	1755	2100
512	8700	6200	6340	4100
4096	9030	6200	6650	4150

For float64, the gc compiler already autovectorizes the plain reduction loop, and on POWER9 that beats the VSX kernel (VSX was ~0.82× at n≥512), so per this library's dispatch principle — never run a SIMD kernel where it loses to the scalar/autovectorized path — there is no float64 VSX kernel on ppc64le; the reductions route to the autovectorizable loop instead. For float32, VSX wins (~1.55–1.61× over naive for n≥64), so the float32 kernel is kept (POWER9-gated). Correctness on ppc64le is validated under QEMU and on real POWER9 silicon.

s390x (IBM z15, VXE2, measured 2026-07-03) — linux1, go1.26.4, -count=6, Dot throughput GB/s, higher is better. For float64 Dot, the gc compiler already autovectorizes the plain reduction loop, and on z15 that beats the VFMADB VX kernel (VX 11.6 GB/s vs the autovectorized loop's 20.6 GB/s = 0.56×). Per this library's dispatch principle — never run a SIMD kernel where it loses to the scalar/autovectorized path — there is no float64 Dot VX kernel on s390x; Dot routes to the autovectorizable loop instead, so it now runs at 1.00× parity (20.6 GB/s) rather than the pre-fix 0.56× loss. Sum and SumSqDiff keep their VFMADB VX kernels (not measured as losses on z15), and Distance/CosineSimilarity build on those and are unchanged; float32 uses the lane-blocked scalar reference as before. Correctness on s390x is validated under QEMU and on real z15 silicon.

License

BSD-3-Clause. © The go-simd authors.

Documentation

Overview

Package floats provides SIMD-accelerated reductions over float32 and float64 slices — the dot product, sum, min, max, Euclidean (L2) distance and cosine similarity that sit at the heart of vector search, embeddings and numerics. It is pure Go (CGO_ENABLED=0, stable Go, no GOEXPERIMENT): the hot reduction loops are SIMD assembly generated by go-asmgen on all six of Go's 64-bit SIMD targets — amd64 (AVX2+FMA), arm64 (NEON FMLA), ppc64le (VSX), s390x (vector facility), riscv64 (RVV) and loong64 (LSX) — with a portable scalar fallback everywhere else, all on a plain `go build`.

Operations

For both []float32 (the Float32 functions) and []float64:

Dot(a, b)              a·b = Σ aᵢbᵢ
Sum(a)                 Σ aᵢ
Min(a) / Max(a)        the minimum / maximum element
Distance(a, b)         ‖a−b‖₂ = √Σ(aᵢ−bᵢ)²
CosineSimilarity(a, b) a·b / (‖a‖₂‖b‖₂)

Dot, Distance and CosineSimilarity require len(a) == len(b) and panic otherwise. Min and Max panic on an empty slice. Sum of an empty slice is 0.

For []float64, the in-place elementwise updates (BLAS level 1) behind iterative solvers, QR and LU:

Axpy(α, x, y)          yᵢ = fma(α, xᵢ, yᵢ)          (y += αx)
Xpby(x, β, y)          yᵢ = fma(β, yᵢ, xᵢ)          (y = x + βy)
AxpyDot(α, x, y, w)    y += αx, then returns Σ yᵢwᵢ  (one pass)
Scale(α, x)            xᵢ = α·xᵢ

They operate on the first len(x) elements and panic if y (or w) is shorter. Unlike the reductions, the updated slices are bit-identical on every architecture: each element is rounded once, exactly as math.FMA rounds it (AVX2 FMA on amd64, NEON FMLA on arm64, the math.FMA loop elsewhere). Only AxpyDot's returned sum is a reduction and follows the rules below.

Floating-point determinism

A SIMD reduction does not add its terms left-to-right: it keeps several independent accumulator lanes (for instruction-level parallelism) and folds them together at the end, so the rounding differs from a naïve scalar loop. To stay deterministic and architecture-independent, the scalar reference used internally and by the tests (sumLanes / dotLanes / …) reduces in the *same* lane-blocked order the kernels use (laneCount lanes), and the SIMD result is required to agree with it to a tight, documented ULP bound across every architecture (the differential tests and FuzzDot/FuzzDistance enforce this). For the similarity use case (cosine / Euclidean distance) bit-exact equality with a naïve scalar loop is neither required nor meaningful; correctness to a few ULP is.

Index

• func Axpy(alpha float64, x, y []float64)
• func AxpyDot(alpha float64, x, y, w []float64) float64
• func CosineSimilarity(a, b []float64) float64
• func Distance(a, b []float64) float64
• func Dot(a, b []float64) float64
• func Float32CosineSimilarity(a, b []float32) float32
• func Float32Distance(a, b []float32) float32
• func Float32Dot(a, b []float32) float32
• func Float32Max(a []float32) float32
• func Float32Min(a []float32) float32
• func Float32Sum(a []float32) float32
• func Max(a []float64) float64
• func Min(a []float64) float64
• func Scale(alpha float64, x []float64)
• func Sum(a []float64) float64
• func Xpby(x []float64, beta float64, y []float64)

Functions

func Axpy

func Axpy(alpha float64, x, y []float64)

Axpy computes y[i] += alpha*x[i] for i < len(x), fused: y[i] = math.FMA(alpha, x[i], y[i]). Elements of y past len(x) are left untouched. It panics if len(y) < len(x).

func AxpyDot

func AxpyDot(alpha float64, x, y, w []float64) float64

AxpyDot computes y[i] += alpha*x[i] for i < len(x) exactly as Axpy does and returns Σ y[i]*w[i] over the updated elements, in one pass over the three slices — the residual update and its norm in conjugate gradient (r -= αAp; rr = r·r is AxpyDot(-α, Ap, r, r)). It panics if len(y) < len(x) or len(w) < len(x).

func CosineSimilarity

func CosineSimilarity(a, b []float64) float64

CosineSimilarity returns a·b / (‖a‖₂‖b‖₂), a value in [-1, 1] measuring the angle between a and b — the standard relevance score for embedding vectors. It panics if len(a) != len(b). If either vector has zero magnitude the result is NaN (0/0), matching the mathematical definition.

func Distance

func Distance(a, b []float64) float64

Distance returns the Euclidean (L2) distance ‖a−b‖₂ = √Σ(aᵢ−bᵢ)². It panics if len(a) != len(b).

func Dot

func Dot(a, b []float64) float64

Dot returns the dot product a·b = Σ aᵢbᵢ. It panics if len(a) != len(b).

func Float32CosineSimilarity

func Float32CosineSimilarity(a, b []float32) float32

Float32CosineSimilarity returns a·b / (‖a‖₂‖b‖₂) over float32 slices — the relevance score for float32 embedding vectors. It panics if len(a) != len(b). If either vector has zero magnitude the result is NaN.

func Float32Distance

func Float32Distance(a, b []float32) float32

Float32Distance returns the Euclidean (L2) distance ‖a−b‖₂ over float32 slices. The squared-difference sum is accumulated in float32 and the square root taken in float64 for accuracy. It panics if len(a) != len(b).

func Float32Dot

func Float32Dot(a, b []float32) float32

Float32Dot returns the dot product a·b = Σ aᵢbᵢ over float32 slices, accumulated in float32. It panics if len(a) != len(b).

func Float32Max

func Float32Max(a []float32) float32

Float32Max returns the largest element of a. It panics if a is empty. If a contains a NaN the result is NaN.

func Float32Min

func Float32Min(a []float32) float32

Float32Min returns the smallest element of a. It panics if a is empty. If a contains a NaN the result is NaN.

func Float32Sum

func Float32Sum(a []float32) float32

Float32Sum returns Σ aᵢ over a float32 slice, accumulated in float32. The sum of an empty slice is 0.

func Max

func Max(a []float64) float64

Max returns the largest element of a. It panics if a is empty. If a contains a NaN the result is NaN.

func Min

func Min(a []float64) float64

Min returns the smallest element of a. It panics if a is empty. If a contains a NaN the result is NaN.

func Scale

func Scale(alpha float64, x []float64)

Scale computes x[i] *= alpha for every element of x.

func Sum

func Sum(a []float64) float64

Sum returns Σ aᵢ. The sum of an empty slice is 0.

func Xpby

func Xpby(x []float64, beta float64, y []float64)

Xpby computes y[i] = x[i] + beta*y[i] for i < len(x), fused: y[i] = math.FMA(beta, y[i], x[i]). It is the search-direction update of conjugate gradient (p = r + βp). It panics if len(y) < len(x).