Awesome Go

simd

CategoryScience and Data Analysis
SubcategoryScience and Data Analysis
Stars26

Native Go vector and SIMD operations on slices with multi-architecture assembly acceleration

About simd

This large README is shown as a plain-text preview. Read the full README at the source

# simd

[![Go Reference](https://pkg.go.dev/badge/github.com/tphakala/simd.svg)](https://pkg.go.dev/github.com/tphakala/simd)
[![codecov](https://codecov.io/gh/tphakala/simd/branch/main/graph/badge.svg)](https://codecov.io/gh/tphakala/simd)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)

A high-performance SIMD (Single Instruction, Multiple Data) library for Go providing vectorized operations on float64, float32, float16, int32, int16, int8, complex128, and complex64 slices.

## Features

- **Pure Go assembly** - Native Go assembler, simple cross-compilation
- **Runtime CPU detection** - Automatically selects optimal implementation (AVX-512, AVX+FMA, AVX without FMA, SSE2, NEON, NEON+FP16, or pure Go); the minimum amd64 SIMD tier is per-package (see [Architecture Support](#architecture-support))
- **Zero allocations** - All operations work on pre-allocated slices
- **150+ operations** - Arithmetic, reduction, statistical, vector, signal processing, activation functions, integer DSP, and complex number operations
- **Multi-architecture** - AMD64 (AVX-512/AVX+FMA/AVX/SSE2, c64 needs SSE4.1) and ARM64 (NEON/NEON+FP16) with pure Go fallback
- **Half-precision support** - Native FP16 SIMD on ARM64 with FP16 extension (Apple Silicon, Cortex-A55+); F16C-accelerated conversions on AMD64
- **Tunable dispatch** - `SIMD_DISABLE` env var masks feature tiers at startup (avoid AVX-512 downclocking, exercise lower tiers, benchmark tier-vs-tier)
- **Thread-safe** - All functions are safe for concurrent use

## Installation

```bash
go get github.com/tphakala/simd
```

Requires Go 1.25+

## Quick Start

```go
package main

import (
    "fmt"
    "github.com/tphakala/simd/cpu"
    "github.com/tphakala/simd/f64"
)

func main() {
    fmt.Println("SIMD:", cpu.Info())

    // Vector operations
    a := []float64{1, 2, 3, 4, 5, 6, 7, 8}
    b := []float64{8, 7, 6, 5, 4, 3, 2, 1}

    // Dot product
    dot := f64.DotProduct(a, b)
    fmt.Println("Dot product:", dot) // 120

    // Element-wise operations
    dst := make([]float64, len(a))
    f64.Add(dst, a, b)
    fmt.Println("Sum:", dst) // [9, 9, 9, 9, 9, 9, 9, 9]

    // Statistical operations
    mean := f64.Mean(a)
    stddev := f64.StdDev(a)
    fmt.Printf("Mean: %.2f, StdDev: %.2f\n", mean, stddev)

    // Vector operations
    f64.Normalize(dst, a)
    fmt.Println("Normalized:", dst)

    // Distance calculation
    dist := f64.EuclideanDistance(a, b)
    fmt.Println("Distance:", dist)
}
```

## Packages

### `cpu` - CPU Feature Detection

```go
import "github.com/tphakala/simd/cpu"

fmt.Println(cpu.Info())        // "AMD64 AVX-512", "AMD64 AVX+FMA", "AMD64 AVX", "AMD64 SSE2", "AMD64 (scalar)", "ARM64 NEON+FP16", or "ARM64 NEON"
                               // SVE-capable ARM64 hosts append " (SVE detected, unused)" - the library runs the NEON path
fmt.Println(cpu.HasAVX())      // true/false
fmt.Println(cpu.HasAVX2())     // true/false
fmt.Println(cpu.HasFMA())      // true/false
fmt.Println(cpu.HasAVX512VL()) // true/false (AVX-512 F+VL)
fmt.Println(cpu.HasNEON())     // true/false
fmt.Println(cpu.HasFP16())     // true/false (ARM64 half-precision SIMD)
fmt.Println(cpu.HasPCLMULQDQ()) // true/false (x86 carry-less multiply)
fmt.Println(cpu.HasF16C())     // true/false (x86 half<->single conversion)
fmt.Println(cpu.HasPMULL())    // true/false (ARM64 polynomial multiply)
```

#### Disabling feature tiers with `SIMD_DISABLE`

Set the `SIMD_DISABLE` environment variable before the process starts to mask
detected CPU features. This is useful for forcing a lower tier on parts where
heavy AVX-512 use causes frequency downclocking, exercising the SSE2/NEON/pure-Go
paths locally, and benchmarking tiers against each other on one machine.

The value is a comma-separated, case-insensitive list of tokens, read once at
program start. Each token clears its own flag plus everything that depends on it:

| Token       | Clears                                    |
| ----------- | ----------------------------------------- |
| `avx512`    | AVX512F, AVX512VL                         |
| `avxvnni`   | AVXVNNI only                              |
| `avx2`      | AVX2, AVXVNNI (and the `avx512` set)      |
| `avx`       | AVX, FMA, F16C (and the `avx2` set)       |
| `fma`       | FMA only                                  |
| `sse42`     | SSE42 (and the `avx` set)                 |
| `sse41`     | SSE41 (and the `sse42` set)               |
| `ssse3`     | SSSE3 (and the `sse41` set)               |
| `sse3`      | SSE3 (and the `ssse3` set)                |
| `pclmulqdq` | PCLMULQDQ only                            |
| `neon`      | NEON, FP16, SVE, SVE2, PMULL, DOTPROD     |
| `fp16`      | FP16 only                                 |
| `sve`       | SVE, SVE2                                 |
| `pmull`     | PMULL only                                |
| `dotprod`   | DOTPROD only                              |
| `all`       | every flag (forces the pure-Go path)      |

F16C is VEX-encoded and only detected alongside AVX, so it clears with the `avx`
cascade (and therefore with every `sse*` token and `all`); `avx2`, `fma`, and
`avx512` sit above AVX and leave F16C set. AVX-VNNI (the VEX-encoded `VPDPWSSD`
tier `i16.XCorr` uses) is detected only alongside AVX2 and its dispatch sits
above it, so `avx2` and every token that cascades through it also clear AVXVNNI,
while the `avxvnni` token clears only that tier (handy for A/B'ing the VNNI
`XCorr` kernel against the plain AVX2 one on one machine).

Unknown tokens are ignored (the library never panics or writes to stderr on env
input). `cpu.Info()` reflects the cleared flags.

```sh
SIMD_DISABLE=avx512 go test ./...   # run as if the CPU had no AVX-512
SIMD_DISABLE=all go test ./...      # force the pure-Go path everywhere
```

The variable must be set before the process starts; it cannot be toggled at
runtime, because the SIMD packages cache their selected kernels during package
init based on the features visible at that moment (function pointers on amd64,
capability flags on arm64).

### `crc` - Cyclic Redundancy Checks

```go
import "github.com/tphakala/simd/crc"

// CRC-16 over poly 0x8005 (init 0, MSB-first, no reflection), the unreflected
// 0x8005 parameterization FLAC uses; folded 16 bytes at a time with PCLMULQDQ
// (amd64) / PMULL (arm64), scalar slice-by-16 fallback.
sum := crc.Checksum16(p) // bit-identical to the scalar reference, zero-alloc
```

| Function          | Description                              | Acceleration                       |
| ----------------- | ---------------------------------------- | ---------------------------------- |
| `Checksum16(p)`   | CRC-16 (poly 0x8005, MSB-first; used by FLAC) | PCLMULQDQ / PMULL carry-less fold |

### `f64` - float64 Operations

**Scope:** `f64` carries the FLAC/LPC and scientific double-precision surface,
including `Autocorrelate` (lag-vectorized LPC autocorrelation) and the split-format
FFT butterfly building blocks (`ButterflyComplex`, `ButterflyComplexStage`,
`ButterflyComplexStage4`, `RealFFTUnpack`, `RealFFTPower`) that a
double-precision FFT/STFT path needs. The broader
audio/ML helpers (PCM conversions, the general split-format complex ops such as
`MulComplex` / `AbsSqComplex`, indexed/strided dot products) live in `f32`
instead, so the two float surfaces remain intentionally asymmetric.

Prefer `ButterflyComplexStage` over `ButterflyComplex` when driving a whole
transform. It takes the stage rather than one block, so the per-block call
overhead disappears and the short spans that cannot fill a vector along `j`
vectorize across blocks instead. Measured over one stage of a 1024-point
transform on an x86 core, that collapses the cost spread across spans from
roughly 30x to under 2x, so the small-span stages stop dominating the transform.

| Category        | Function                            | Description                   | SIMD Width                          |
| --------------- | ----------------------------------- | ----------------------------- | ----------------------------------- |
| **Arithmetic**  | `Add(dst, a, b)`                    | Element-wise addition         | 8x (AVX-512) / 4x (AVX) / 2x (NEON) |
|                 | `Sub(dst, a, b)`                    | Element-wise subtraction      | 8x / 4x / 2x                        |
|                 | `Mul(dst, a, b)`                    | Element-wise multiplication   | 8x / 4x / 2x                        |
|                 | `Div(dst, a, b)`                    | Element-wise division         | 8x / 4x / 2x                        |
|                 | `Scale(dst, a, s)`                  | Multiply by scalar            | 8x / 4x / 2x                        |
|                 | `AddScalar(dst, a, s)`              | Add scalar                    | 8x / 4x / 2x                        |
|                 | `SubFromScalar(dst, a, s)`          | Scalar minus vector           | 8x / 4x / 2x (composed SIMD)        |
|                 | `FMA(dst, a, b, c)`                 | Fused multiply-add: a\*b+c    | 8x / 4x / 2x                        |
|                 | `MulAdd(dst, a, b)`                 | Fused accumulate: dst += a\*b | 8x / 4x / 2x                        |
|                 | `AddScaled(dst, alpha, s)`          | dst += alpha\*s (axpy)        | 8x / 4x / 2x                        |
| **Unary**       | `Abs(dst, a)`                       | Absolute value                | 8x / 4x / 2x                        |
|                 | `Neg(dst, a)`                       | Negation                      | 8x / 4x / 2x                        |
|                 | `Sqrt(dst, a)`                      | Square root                   | 8x / 4x / 2x                        |
|                 | `Reciprocal(dst, a)`                | Reciprocal (1/x)              | 8x / 4x / 2x                        |
|                 | `Round(dst, src)`                   | Round half away from zero     | 4x (AVX) / 2x (NEON) / Go fallback  |
| **Reduction**   | `DotProduct(a, b)`                  | Dot product                   | 8x / 4x / 2x                        |
|                 | `WeightedSum(w, src)`               | Weighted sum Σ(wᵢ·srcᵢ)       | 8x / 4x / 2x                        |
|                 | `SumOfSquares(src)`                 | Sum of squares Σ(srcᵢ²)       | 8x / 4x / 2x                        |
|                 | `Sum(a)`                            | Sum of elements               | 8x / 4x / 2x                        |
|                 | `Min(a)`                            | Minimum value                 | 8x / 4x / 2x                        |
|                 | `Max(a)`                            | Maximum value                 | 8x / 4x / 2x                        |
|                 | `MaxAbs(a)`                         | Max absolute value (∞-norm)   | 8x / 4x / 2x                        |
|                 | `MinIdx(a)`                         | Index of minimum value        | Pure Go                             |
|                 | `MaxIdx(a)`                         | Index of maximum value        | Pure Go                             |
| **Statistical** | `Mean(a)`                           | Arithmetic mean               | 8x / 4x / 2x                        |
|                 | `Variance(a)`                       | Population variance           | 8x / 4x / 2x                        |
|                 | `StdDev(a)`                         | Standard deviation            | 8x / 4x / 2x                        |
| **Vector**      | `EuclideanDistance(a, b)`           | L2 distance                   | 8x / 4x / 2x                        |
|                 | `Normalize(dst, a)`                 | Unit vector normalization     | 8x / 4x / 2x                        |
|                 | `CumulativeSum(dst, a)`             | Running sum                   | Sequential                          |
| **Range**       | `Clamp(dst, a, min, max)`           | Clamp to range                | 8x / 4x / 2x                        |
| **Activation**  | `Sigmoid(dst, src)`                 | Sigmoid: 1/(1+e^-x)           | 4x (AVX2) / 2x (NEON)               |
|                 | `ReLU(dst, src)`                    | Rectified Linear Unit         | 4x (AVX) / 2x (NEON)                |
|                 | `Tanh(dst, src)`                    | Hyperbolic tangent            | 4x (AVX2) / 2x (NEON)               |
|                 | `Exp(dst, src)`                     | Exponential e^x               | 4x (AVX2) / 2x (NEON)               |
|                 | `ClampScale(dst, src, min, max, s)` | Fused clamp and scale         | 4x (AVX) / 2x (NEON)                |
| **Transcendental** | `Log(dst, src)`                  | Natural log ln(x)             | 4x (AVX2+FMA) / 2x (NEON)           |
|                 | `Log2(dst, src)` / `Log10(dst, src)`| Base-2 / base-10 log          | 4x (AVX2+FMA) / 2x (NEON)           |
|                 | `Pow(dst, src, exp)`                | Scalar power x^exp (PCEN, dB) | 4x (AVX2+FMA) / 2x (NEON)           |
|                 | `PowElem(dst, base, exp)`           | Elementwise base^exp          | 4x (AVX2+FMA) / 2x (NEON)           |
| **Batch**       | `DotProductBatch(r, rows, v)`       | Multiple dot products         | 8x / 4x / 2x                        |
| **Signal**      | `ConvolveValid(dst, sig, k)`        | FIR filter / convolution      | 8x / 4x / 2x                        |
|                 | `ConvolveValidMulti(dsts, sig, ks)` | Multi-kernel convolution      | 8x / 4x / 2x                        |
|                 | `ConvolveValidMaxAbs(sig, k)`       | Fused FIR abs-max peak (no scratch) | 8x / 4x / 2x                  |
|                 | `ConvolveValidMaxAbsMulti(sig, ks)` | Multi-kernel abs-max peak (true-peak) | 8x / 4x / 2x                |
|                 | `ConvolveDecimate(dst,sig,k,f,p)`   | Strided FIR downsample (decimate) | 8x / 4x / 2x                    |
|                 | `AccumulateAdd(dst, src, off)`      | Overlap-add: dst[off:] += src | 8x / 4x / 2x                        |
|                 | `Autocorrelate(autoc, x, maxLag)`   | LPC autocorrelation Σ x[i]·x[i-lag] (bit-exact) | 4x (AVX2) / 2x (NEON)     |
| **Complex/FFT** | `ButterflyComplex(uRe,uIm,lRe,lIm,twRe,twIm)` | FFT butterfly with twiddle multiply | 4x (AVX+FMA) / 2x (NEON)   |
|                 | `ButterflyComplexStage(re,im,span,twRe,twIm)` | One whole radix-2 DIT stage (any span) | 4x (AVX+FMA) / 2x (NEON)   |
|                 | `ButterflyComplexStage4(re,im,span,tw1..tw3)` | One whole radix-4 DIT stage (two radix-2 stages in one pass) | 4x (AVX+FMA) / 2x (NEON)   |
|                 | `RealFFTUnpack(outRe,outIm,zRe,zIm,twRe,twIm)` | Real-FFT even/odd unpack step | 4x (AVX2+FMA) / 2x (NEON)   |
|                 | `RealFFTPower(dst,zRe,zIm,twRe,twIm)`         | Fused real-FFT power spectrum \|X_k\|^2 (single pass) | 4x (AVX2+FMA) / 2x (NEON)   |
| **Audio**       | `Interleave2(dst, a, b)`            | Pack stereo: [L,R,L,R,...]    | 4x / 2x                             |
|                 | `Deinterleave2(a, b, src)`          | Unpack stereo to channels     | 4x / 2x                             |
|                 | `InterleaveN(dst, srcs)`            | Pack N planar streams (any N; N-stream Interleave2) | N=2,4,8 AVX, N=3,6 AVX2 / N=2,3,4 NEON; else Go |
|                 | `DeinterleaveN(dsts, src)`          | Unpack N interleaved streams (any N) | N=2,4,8 AVX, N=3,6 AVX2 / N=2,3,4 NEON; else Go |
|                 | `CubicInterpDot(hist,a,b,c,d,x)`    | Fused cubic interp dot product| 4x / 2x                             |

`DotProductBatch` scores its `[][]float64` rows in groups of four, keeping the
query vector resident in registers across each group via a fused 4-row kernel on
AMD64 (AVX-512 and AVX+FMA) and ARM64 NEON instead of re-loading it per row.
Short, ragged, or sub-SIMD-width rows fall back to the per-row dot product, with
identical results.

`Autocorrelate` computes the LPC autocorrelation `autoc[lag] = Σ x[i]·x[i-lag]`
used by FLAC-style encoders. It vectorizes across lags (one accumulator lane per
lag, never fusing the multiply-add), so each lag's sum keeps the exact left-to-right
order of the scalar loop and every build emits byte-identical results to the pure-Go
reference. The AVX2 path accumulates four lags per YMM, NEON two lags per V register;
non-AVX2/NEON CPUs and short blocks use the scalar reference.

#### STFT (fused real-input short-time Fourier transform)

`STFTPlan` is the spectral front-end's missing middle: the library already covers
the post-FFT power spectrum (`c128.AbsSq`), mel projection (`DotProductBatch`),
and PCEN / log-mel normalization (`Exp`, `Mul`, `Log`), but not the transform.

Both `f64` and `f32` provide it (with `complex64` output for `f32`).

```go
plan, _ := f64.NewSTFTPlan(1024)               // power-of-two nfft; reuse across calls
bins := plan.NumBins()                         // nfft/2 + 1 (Hermitian half-spectrum)
nFrames := plan.NumFrames(len(signal), hop, f64.PadZero)

spec := make([][]complex128, nFrames)          // caller-owned output, one row per frame
for i := range spec { spec[i] = make([]complex128, bins) }
plan.STFT(spec, signal, hann, hop, f64.PadZero) // fills spec; returns frames written

// Flat, frame-contiguous power (stride NumBins) feeds DotProductBatch directly
// as a mel-filterbank projection, with no per-frame allocation:
power := make([]float64, nFrames*bins)
plan.STFTPowerInto(power, signal, hann, hop, f64.PadZero)
for f := range nFrames {
    f64.DotProductBatch(mel[f], filterbank, power[f*bins:(f+1)*bins])
}
```

The transform uses a half-length complex FFT (rfft, ~2x cheaper than a full
complex FFT), keeps the twiddle/bit-reversal plan resident, and fuses the window
multiply into the frame pack (and the `|.|^2` power step in `STFTPower` /
`STFTPowerInto`). The `PadMode` argument selects the framing convention: `NoPad`
is the no-padding case (frame `f` is `signal[f*hop : f*hop+nfft]`, matching
librosa `stft(..., center=False)`), while `PadZero` and `PadReflect` center each
frame with `nfft/2` of zero or reflect padding per side, matching librosa
`center=True` (`pad_mode="constant"` / `"reflect"`). `NumFrames` reports the frame
count for a given pad mode so you can size buffers. The centered output is pinned
against a librosa golden vector in the tests. The plan is allocation-free across
calls; a plan holds transform scratch, so use one plan per goroutine. The
transform is a power-of-two `nfft` rfft whose every arithmetic pass over the
frame runs through the vector primitives: the frame is packed with
`Deinterleave2` and `Mul`, the FFT core is radix-4 (`ButterflyComplexStage4`,
with at most one trailing radix-2 `ButterflyComplexStage`), and the real-input
unravel is `RealFFTUnpack` or, for the power spectrum, `RealFFTPower`; the
bit-reversal reorder, padded edge frames and rows shorter than `NumBins` stay
scalar. The output is tolerance-stable, not bit-stable, across CPU tiers and
library versions.

The inverse direction is available too. `RFFT` and `IRFFT` are the single-frame
forward and inverse real FFT: `RFFT` transforms one windowed `nfft`-sample frame
to its Hermitian half-spectrum (bit-for-bit a `NoPad` STFT row, same pipeline),
and `IRFFT` inverts it back to `nfft` real samples scaled by `1/nfft`, following
the `numpy.fft.irfft` convention that the imaginary parts of the DC and Nyquist
bins are ignored, so `IRFFT(RFFT(x))` round-trips within tolerance. `ISTFT` is
the batched synthesis inverse of `STFT`: it overlap-adds the windowed inverse
frames and normalizes by the squared-window overlap (the librosa/scipy WOLA
convention), trimming per the same `PadMode`, so `ISTFT(STFT(x))` reconstructs
`x` for any window and hop whose squared overlap never vanishes. All three are
allocation-free and reuse the plan scratch (so one plan per goroutine), and

Frequently Asked Questions

What is simd?

simd is a Science and Data Analysis library for the Go programming language. Native Go vector and SIMD operations on slices with multi-architecture assembly acceleration

How do I install simd?

Install simd with the Go module system using `go get tphakala/simd`. Check the repository for the current installation instructions.

What category does simd belong to?

simd is listed under Science and Data Analysis, specifically Science and Data Analysis.

← Back to Science and Data Analysis