# Elementwise and Reductions

> NumKong for C

`#include <numkong/scalar.h>`

`#include <numkong/each.h>`

`#include <numkong/reduce.h>`

`#include <numkong/trigonometry.h>`

## nk_*_fma_best

Scalar fused multiply-add, a × b + c.

```c
nk_f32_t nk_f32_fma_best(
    nk_f32_t a,
    nk_f32_t b,
    nk_f32_t c,
    nk_capability_t capabilities);
```

- `nk_f32_fma_best`: a, b, c `nk_f32_t`
- `nk_f64_fma_best`: a, b, c `nk_f64_t`
- `nk_f16_fma_best`: a, b, c `nk_f16_t`

- `nk_f32_fma` backends: serial, neon, haswell, rvv, v128relaxed, powervsx
- `nk_f64_fma` backends: serial, neon, haswell, rvv, v128relaxed, powervsx
- `nk_f16_fma` backends: serial, neonhalf, haswell, sapphire

- `a` (in): Multiplicand.
- `b` (in): Multiplier.
- `c` (in): Addend.
- `capabilities` (in): The CPU's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports; serial always stays.

See also: C++ `std::fma`

See also: Rust `f32::mul_add`

## nk_each_scale_*_best

Element-wise scale with shift: result[i] = alpha * a[i] + beta.

```c
nk_status_t nk_each_scale_f64_best(
    nk_f64_t const *a,
    nk_size_t n,
    nk_f64_t const *alpha,
    nk_f64_t const *beta,
    nk_f64_t *result,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_scale_f64_best`: a, result `nk_f64_t`, alpha, beta `nk_f64_t`
- `nk_each_scale_f32_best`: a, result `nk_f32_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_f16_best`: a, result `nk_f16_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_bf16_best`: a, result `nk_bf16_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_i8_best`: a, result `nk_i8_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_u8_best`: a, result `nk_u8_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_i16_best`: a, result `nk_i16_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_u16_best`: a, result `nk_u16_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_i32_best`: a, result `nk_i32_t`, alpha, beta `nk_f64_t`
- `nk_each_scale_u32_best`: a, result `nk_u32_t`, alpha, beta `nk_f64_t`
- `nk_each_scale_i64_best`: a, result `nk_i64_t`, alpha, beta `nk_f64_t`
- `nk_each_scale_u64_best`: a, result `nk_u64_t`, alpha, beta `nk_f64_t`
- `nk_each_scale_e4m3_best`: a, result `nk_e4m3_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_e5m2_best`: a, result `nk_e5m2_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_e2m3_best`: a, result `nk_e2m3_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_e3m2_best`: a, result `nk_e3m2_t`, alpha, beta `nk_f32_t`
- `nk_each_scale_f32c_best`: a, result `nk_f32c_t`, alpha, beta `nk_f32c_t`
- `nk_each_scale_f64c_best`: a, result `nk_f64c_t`, alpha, beta `nk_f64c_t`

- `nk_each_scale_f64` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_f32` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_scale_f16` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_scale_bf16` backends: serial, neonbfdot, v128relaxed, haswell, skylake, rvv, cuda, ampere
- `nk_each_scale_i8` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_scale_u8` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_scale_i16` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_u16` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_i32` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_u32` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_i64` backends: serial, neon, skylake, rvv, cuda
- `nk_each_scale_u64` backends: serial, neon, skylake, rvv, cuda
- `nk_each_scale_e4m3` backends: serial, neon, haswell, skylake, rvv, cuda, ada
- `nk_each_scale_e5m2` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_e2m3` backends: serial, neon, cuda
- `nk_each_scale_e3m2` backends: serial, neon, cuda
- `nk_each_scale_f32c` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_scale_f64c` backends: serial, neon, haswell, skylake, rvv, cuda

- `a` (in): The input vector.
- `n` (in): The number of elements in the vector.
- `alpha` (in): Pointer to the scaling factor (type depends on input precision).
- `beta` (in): Pointer to the shift (bias) value (type depends on input precision).
- `result` (out): The output vector.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_sum_*_best

Element-wise sum: result[i] = a[i] + b[i].

```c
nk_status_t nk_each_sum_f64_best(
    nk_f64_t const *a,
    nk_f64_t const *b,
    nk_size_t n,
    nk_f64_t *result,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_sum_f64_best`: a, b, result `nk_f64_t`
- `nk_each_sum_f32_best`: a, b, result `nk_f32_t`
- `nk_each_sum_f16_best`: a, b, result `nk_f16_t`
- `nk_each_sum_bf16_best`: a, b, result `nk_bf16_t`
- `nk_each_sum_i8_best`: a, b, result `nk_i8_t`
- `nk_each_sum_u8_best`: a, b, result `nk_u8_t`
- `nk_each_sum_i16_best`: a, b, result `nk_i16_t`
- `nk_each_sum_u16_best`: a, b, result `nk_u16_t`
- `nk_each_sum_i32_best`: a, b, result `nk_i32_t`
- `nk_each_sum_u32_best`: a, b, result `nk_u32_t`
- `nk_each_sum_i64_best`: a, b, result `nk_i64_t`
- `nk_each_sum_u64_best`: a, b, result `nk_u64_t`
- `nk_each_sum_e4m3_best`: a, b, result `nk_e4m3_t`
- `nk_each_sum_e5m2_best`: a, b, result `nk_e5m2_t`
- `nk_each_sum_e2m3_best`: a, b, result `nk_e2m3_t`
- `nk_each_sum_e3m2_best`: a, b, result `nk_e3m2_t`
- `nk_each_sum_f32c_best`: a, b, result `nk_f32c_t`
- `nk_each_sum_f64c_best`: a, b, result `nk_f64c_t`

- `nk_each_sum_f64` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_sum_f32` backends: serial, neon, v128, haswell, skylake, rvv, cuda, rocm
- `nk_each_sum_f16` backends: serial, neonhalf, v128relaxed, haswell, sapphire, rvv, cuda, rocm
- `nk_each_sum_bf16` backends: serial, neonbfdot, v128, haswell, skylake, rvv, cuda, ampere, rocm
- `nk_each_sum_i8` backends: serial, neon, v128, haswell, icelake, rvv, cuda
- `nk_each_sum_u8` backends: serial, neon, v128, haswell, icelake, rvv, cuda
- `nk_each_sum_i16` backends: serial, neon, haswell, icelake, rvv, cuda
- `nk_each_sum_u16` backends: serial, neon, haswell, icelake, rvv, cuda
- `nk_each_sum_i32` backends: serial, neon, haswell, icelake, rvv, cuda
- `nk_each_sum_u32` backends: serial, neon, haswell, icelake, rvv, cuda
- `nk_each_sum_i64` backends: serial, neon, icelake, rvv, cuda
- `nk_each_sum_u64` backends: serial, neon, icelake, rvv, cuda
- `nk_each_sum_e4m3` backends: serial, neon, neonfhm, haswell, skylake, sapphire, rvv, cuda, ada
- `nk_each_sum_e5m2` backends: serial, neon, neonfhm, haswell, skylake, rvv, cuda
- `nk_each_sum_e2m3` backends: serial, neon, cuda
- `nk_each_sum_e3m2` backends: serial, neon, cuda
- `nk_each_sum_f32c` backends: serial, cuda
- `nk_each_sum_f64c` backends: serial, cuda

- `a` (in): The first input vector.
- `b` (in): The second input vector.
- `n` (in): The number of elements in the vectors.
- `result` (out): The output vector.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_blend_*_best

Weighted sum: result[i] = alpha * a[i] + beta * b[i].

```c
nk_status_t nk_each_blend_f64_best(
    nk_f64_t const *a,
    nk_f64_t const *b,
    nk_size_t n,
    nk_f64_t const *alpha,
    nk_f64_t const *beta,
    nk_f64_t *result,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_blend_f64_best`: a, b, result `nk_f64_t`, alpha, beta `nk_f64_t`
- `nk_each_blend_f32_best`: a, b, result `nk_f32_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_f16_best`: a, b, result `nk_f16_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_bf16_best`: a, b, result `nk_bf16_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_i8_best`: a, b, result `nk_i8_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_u8_best`: a, b, result `nk_u8_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_i16_best`: a, b, result `nk_i16_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_u16_best`: a, b, result `nk_u16_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_i32_best`: a, b, result `nk_i32_t`, alpha, beta `nk_f64_t`
- `nk_each_blend_u32_best`: a, b, result `nk_u32_t`, alpha, beta `nk_f64_t`
- `nk_each_blend_i64_best`: a, b, result `nk_i64_t`, alpha, beta `nk_f64_t`
- `nk_each_blend_u64_best`: a, b, result `nk_u64_t`, alpha, beta `nk_f64_t`
- `nk_each_blend_e4m3_best`: a, b, result `nk_e4m3_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_e5m2_best`: a, b, result `nk_e5m2_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_e2m3_best`: a, b, result `nk_e2m3_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_e3m2_best`: a, b, result `nk_e3m2_t`, alpha, beta `nk_f32_t`
- `nk_each_blend_f32c_best`: a, b, result `nk_f32c_t`, alpha, beta `nk_f32c_t`
- `nk_each_blend_f64c_best`: a, b, result `nk_f64c_t`, alpha, beta `nk_f64c_t`

- `nk_each_blend_f64` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_blend_f32` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_blend_f16` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_blend_bf16` backends: serial, neonbfdot, v128relaxed, haswell, skylake, rvv, cuda, ampere
- `nk_each_blend_i8` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_blend_u8` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_blend_i16` backends: serial, cuda
- `nk_each_blend_u16` backends: serial, cuda
- `nk_each_blend_i32` backends: serial, cuda
- `nk_each_blend_u32` backends: serial, cuda
- `nk_each_blend_i64` backends: serial, cuda
- `nk_each_blend_u64` backends: serial, cuda
- `nk_each_blend_e4m3` backends: serial, neon, haswell, skylake, rvv, cuda, ada
- `nk_each_blend_e5m2` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_blend_e2m3` backends: serial, neon, cuda
- `nk_each_blend_e3m2` backends: serial, neon, cuda
- `nk_each_blend_f32c` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_blend_f64c` backends: serial, neon, haswell, skylake, rvv, cuda

- `a` (in): The first input vector.
- `b` (in): The second input vector.
- `n` (in): The number of elements in the vectors.
- `alpha` (in): Pointer to the first weight (type depends on input precision).
- `beta` (in): Pointer to the second weight (type depends on input precision).
- `result` (out): The output vector.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_fma_*_best

Fused multiply-add: result[i] = alpha * a[i] * b[i] + beta * c[i].

```c
nk_status_t nk_each_fma_f64_best(
    nk_f64_t const *a,
    nk_f64_t const *b,
    nk_f64_t const *c,
    nk_size_t n,
    nk_f64_t const *alpha,
    nk_f64_t const *beta,
    nk_f64_t *result,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_fma_f64_best`: a, b, c, result `nk_f64_t`, alpha, beta `nk_f64_t`
- `nk_each_fma_f32_best`: a, b, c, result `nk_f32_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_f16_best`: a, b, c, result `nk_f16_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_bf16_best`: a, b, c, result `nk_bf16_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_i8_best`: a, b, c, result `nk_i8_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_u8_best`: a, b, c, result `nk_u8_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_i16_best`: a, b, c, result `nk_i16_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_u16_best`: a, b, c, result `nk_u16_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_i32_best`: a, b, c, result `nk_i32_t`, alpha, beta `nk_f64_t`
- `nk_each_fma_u32_best`: a, b, c, result `nk_u32_t`, alpha, beta `nk_f64_t`
- `nk_each_fma_i64_best`: a, b, c, result `nk_i64_t`, alpha, beta `nk_f64_t`
- `nk_each_fma_u64_best`: a, b, c, result `nk_u64_t`, alpha, beta `nk_f64_t`
- `nk_each_fma_e4m3_best`: a, b, c, result `nk_e4m3_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_e5m2_best`: a, b, c, result `nk_e5m2_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_e2m3_best`: a, b, c, result `nk_e2m3_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_e3m2_best`: a, b, c, result `nk_e3m2_t`, alpha, beta `nk_f32_t`
- `nk_each_fma_f32c_best`: a, b, c, result `nk_f32c_t`, alpha, beta `nk_f32c_t`
- `nk_each_fma_f64c_best`: a, b, c, result `nk_f64c_t`, alpha, beta `nk_f64c_t`

- `nk_each_fma_f64` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_fma_f32` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_fma_f16` backends: serial, neon, neonfhm, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_fma_bf16` backends: serial, neonbfdot, v128relaxed, haswell, skylake, rvv, cuda, ampere
- `nk_each_fma_i8` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_fma_u8` backends: serial, neon, v128relaxed, haswell, skylake, rvv, cuda
- `nk_each_fma_i16` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_fma_u16` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_fma_i32` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_fma_u32` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_fma_i64` backends: serial, neon, skylake, rvv, cuda
- `nk_each_fma_u64` backends: serial, neon, skylake, rvv, cuda
- `nk_each_fma_e4m3` backends: serial, neon, neonfhm, haswell, skylake, rvv, cuda, ada
- `nk_each_fma_e5m2` backends: serial, neon, neonfhm, haswell, skylake, rvv, cuda
- `nk_each_fma_e2m3` backends: serial, neon, cuda
- `nk_each_fma_e3m2` backends: serial, neon, cuda
- `nk_each_fma_f32c` backends: serial, neon, haswell, skylake, rvv, cuda
- `nk_each_fma_f64c` backends: serial, neon, haswell, skylake, rvv, cuda

- `a` (in): The first input vector.
- `b` (in): The second input vector.
- `c` (in): The third input vector.
- `n` (in): The number of elements in the vectors.
- `alpha` (in): Pointer to the scaling factor for a[i] * b[i] (type depends on input precision).
- `beta` (in): Pointer to the scaling factor for c[i] (type depends on input precision).
- `result` (out): The output vector.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_swiglu_*_best

Fused SwiGLU: y = silu(gate) × up, with a scale on the gate and on the output.

```c
nk_status_t nk_each_swiglu_f32_best(
    nk_f32_t const *gate,
    nk_f32_t const *up,
    nk_f32_t *y,
    nk_size_t rows,
    nk_size_t columns,
    nk_size_t gate_stride,
    nk_size_t up_stride,
    nk_size_t y_stride,
    nk_f32_t gate_scale,
    nk_f32_t output_scale,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_swiglu_f32_best`: gate, up, y `nk_f32_t`
- `nk_each_swiglu_f16_best`: gate, up, y `nk_f16_t`
- `nk_each_swiglu_bf16_best`: gate, up, y `nk_bf16_t`
- `nk_each_swiglu_e4m3_best`: gate, up, y `nk_e4m3_t`

- `nk_each_swiglu_f32` backends: neon, serial, haswell, skylake, cuda, rocm, metal
- `nk_each_swiglu_f16` backends: neon, serial, haswell, skylake, metal
- `nk_each_swiglu_bf16` backends: neon, serial, haswell, skylake, cuda, ampere, rocm, metal
- `nk_each_swiglu_e4m3` backends: neon, serial, haswell, skylake, cuda, ada, rocm, metal

SiLU is not linear, so FP8 codes need `gate_scale` to reach true values first, while the up and output scales of FP8 tensors fold into `output_scale` as `up_scale / y_scale`.

- `gate` (in): The gate input matrix of shape `rows` by `columns`.
- `up` (in): The up input matrix, same shape as `gate`; NULL collapses to plain SiLU.
- `y` (out): The output matrix, same shape and dtype as the inputs; may alias `gate`.
- `rows` (in): The number of rows in each matrix.
- `columns` (in): The number of columns in each matrix.
- `gate_stride` (in): The row stride of `gate` in bytes.
- `up_stride` (in): The row stride of `up` in bytes.
- `y_stride` (in): The row stride of `y` in bytes.
- `gate_scale` (in): Multiplies every gate element before the SiLU; 1 for unscaled inputs.
- `output_scale` (in): Multiplies every result before the store; 1 for unscaled inputs.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_rmsnorm_*_best

Grouped RMSNorm: y = x × rsqrt(mean(x²) + epsilon) × gamma.

```c
nk_status_t nk_each_rmsnorm_f32_best(
    nk_f32_t const *x,
    nk_f32_t const *gamma,
    nk_f32_t *y,
    nk_size_t rows,
    nk_size_t groups,
    nk_size_t columns,
    nk_size_t x_stride,
    nk_size_t y_stride,
    nk_f32_t epsilon,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_rmsnorm_f32_best`: x, y `nk_f32_t`
- `nk_each_rmsnorm_f16_best`: x, y `nk_f16_t`
- `nk_each_rmsnorm_bf16_best`: x, y `nk_bf16_t`
- `nk_each_rmsnorm_e4m3_best`: x, y `nk_e4m3_t`

- `nk_each_rmsnorm_f32` backends: neon, serial, haswell, skylake, cuda, rocm
- `nk_each_rmsnorm_f16` backends: neon, serial, haswell, skylake, metal
- `nk_each_rmsnorm_bf16` backends: neon, serial, haswell, skylake, genoa, cuda, ampere, rocm, metal
- `nk_each_rmsnorm_e4m3` backends: neon, serial, haswell, skylake, genoa, cuda, ada, rocm, metal

Scaling x leaves y unchanged but for `epsilon`, so FP8 codes need no input scale: pass epsilon in units of the codes squared, and fold the output scale of FP8 y into `gamma`.

Types narrower than F32 stay in F32 throughout: the sum of squares, its mean, epsilon, the inverse root and the products with x and gamma, so rows whose squares overflow F32 normalize to zeros. F32 rows sum their squares and take the inverse root in F64, then multiply in F32 through a power of two that keeps tiny and huge rows in range.

- `x` (in): Input matrix; each row holds `groups` separately normalized `columns-vectors`.
- `gamma` (in): Per-column gain, length `columns`, shared by groups and rows; NULL is unit.
- `y` (out): Output matrix, same shape and dtype as `x`; may alias `x` for in-place use.
- `rows` (in): Number of rows in the input and output matrices.
- `groups` (in): Number of independent normalization groups per row.
- `columns` (in): Number of columns per group.
- `x_stride` (in): Row stride of `x` in bytes; groups pack `columns` apart.
- `y_stride` (in): Row stride of `y` in bytes.
- `epsilon` (in): Variance epsilon added before the reciprocal square root.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_rmscast_*_best

Grouped RMSNorm of dot products back into the type they multiplied: y = x × rsqrt(mean(x²) + epsilon) × gamma.

```c
nk_status_t nk_each_rmscast_bf16_best(
    nk_f32_t const *x,
    nk_f32_t const *gamma,
    nk_bf16_t *y,
    nk_size_t rows,
    nk_size_t groups,
    nk_size_t columns,
    nk_size_t x_stride,
    nk_size_t y_stride,
    nk_f32_t epsilon,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_each_rmscast_bf16_best`: x `nk_f32_t`, y `nk_bf16_t`
- `nk_each_rmscast_f16_best`: x `nk_f32_t`, y `nk_f16_t`
- `nk_each_rmscast_e4m3_best`: x `nk_f32_t`, y `nk_e4m3_t`
- `nk_each_rmscast_e5m2_best`: x `nk_f32_t`, y `nk_e5m2_t`
- `nk_each_rmscast_e2m3_best`: x `nk_f32_t`, y `nk_e2m3_t`
- `nk_each_rmscast_e3m2_best`: x `nk_f32_t`, y `nk_e3m2_t`
- `nk_each_rmscast_f32_best`: x `nk_f64_t`, y `nk_f32_t`
- `nk_each_rmscast_i8_best`: x `nk_i32_t`, y `nk_i8_t`
- `nk_each_rmscast_u8_best`: x `nk_u32_t`, y `nk_u8_t`

- `nk_each_rmscast_bf16` backends: neon, serial, haswell, skylake, cuda, ampere, rocm, metal
- `nk_each_rmscast_f16` backends: neon, serial, cuda, rocm, metal
- `nk_each_rmscast_e4m3` backends: neon, serial, haswell, skylake, cuda, ada, rocm, metal
- `nk_each_rmscast_e5m2` backends: serial, cuda, ada, rocm, metal
- `nk_each_rmscast_e2m3` backends: serial, cuda, rocm, metal
- `nk_each_rmscast_e3m2` backends: serial, cuda, rocm, metal
- `nk_each_rmscast_f32` backends: neon, serial, haswell, skylake, cuda, rocm
- `nk_each_rmscast_i8` backends: serial, cuda, rocm, metal
- `nk_each_rmscast_u8` backends: serial, cuda, rocm, metal

Reads what [`nk_dot_output_dtype`](https://www.unum.cloud/docs/numkong/c/dot-products.md#nk_dot_output_dtype) names for the output type, F32 for BF16, F16, E4M3, E5M2, E2M3 and E3M2, F64 for F32, I32 for I8 and U32 for U8, normalizes it like [`nk_each_rmsnorm_f32_best`](#nk_each_rmsnorm_f32_best), and rounds every output like `nk_cast` of that pair, NaN payloads included, so a GEMM output becomes the next GEMM's input in one pass. Epsilon adds to mean(x²) of the raw input. Integer outputs are saturate(round_to_nearest_even(x / rms(x) × gamma)), where gamma carries the output quantization scale, and U8 suits only non-negative activations.

F32 inputs stay in F32 throughout, like the RMSNorm of types narrower than F32. F64 inputs stay in F64, from a sum of squares rescaled where squares leave its range through the products with gamma, rounding once into F32. I32 and U32 inputs sum in F64 and round their mean square to F32, where epsilon adds and everything after runs.

- `x` (in): Input matrix of the wider type; each row holds `groups` separate vectors.
- `gamma` (in): Per-column gain, length `columns`, shared by groups and rows; NULL is unit.
- `y` (out): Output matrix, the shape of `x` in the type the kernel names.
- `rows` (in): Number of rows in the input and output matrices.
- `groups` (in): Number of independent normalization groups per row.
- `columns` (in): Number of columns per group.
- `x_stride` (in): Row stride of `x` in bytes; groups pack `columns` apart.
- `y_stride` (in): Row stride of `y` in bytes.
- `epsilon` (in): Variance epsilon added before the reciprocal square root.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_each_scale_input_dtype

Returns the scalar parameter dtype for elementwise scale/blend/fma operations.

```c
nk_dtype_t nk_each_scale_input_dtype(nk_dtype_t dtype);
```

## nk_each_error_bound

Returns the error bound of elementwise scale, blend and FMA before their results round into `dtype:` zero for integers, which match exactly, and else per [`nk_accumulation_error_bound`](https://www.unum.cloud/docs/numkong/c/scalar-types.md#nk_accumulation_error_bound) of their coefficients.

```c
nk_f64_t nk_each_error_bound(nk_dtype_t dtype);
```

## nk_each_rmsnorm_reduction_dtype

Returns the type RMSNorm and RMSCast reduce in, never compensated: one tier above the narrower of `input_dtype` and `output_dtype`, F64 when both are at least 32 bits and F32 otherwise.

```c
nk_dtype_t nk_each_rmsnorm_reduction_dtype(
    nk_dtype_t input_dtype,
    nk_dtype_t output_dtype);
```

## nk_each_rmsnorm_error_bound

Returns the RMSNorm and RMSCast error bound: `columns` positive squares sum within `columns` unit roundoffs, halved by the inverse root into eighths of the tier bound, plus 3 per tier.

```c
nk_f64_t nk_each_rmsnorm_error_bound(
    nk_dtype_t input_dtype,
    nk_dtype_t output_dtype,
    nk_size_t columns);
```

## nk_each_swiglu_error_bound

Returns the error bound of SwiGLU outputs before they round into `dtype:` 5 roundings, through both descales, the exponent, the division and the product with the up value, each per [`nk_each_error_bound`](#nk_each_error_bound).

```c
nk_f64_t nk_each_swiglu_error_bound(nk_dtype_t dtype);
```

## nk_*_find_kernel

Finds the each kernel of `kind` for `dtype`, from the best of `capabilities`.

```c
nk_status_t nk_each_find_kernel(
    nk_kernel_kind_t kind,
    nk_dtype_t dtype,
    nk_capability_t capabilities,
    nk_kernel_punned_t *kernel,
    nk_capability_t *capability);
```

- `nk_each_find_kernel`
- `nk_reduce_find_kernel`
- `nk_trigonometry_find_kernel`

- `kind` (in): The verb whose kernel to find.
- `dtype` (in): The input type the kernel reads.
- `capabilities` (in): The mask to pick the best capability from.
- `kernel` (out): The kernel, or null when no capability in `capabilities` has it.
- `capability` (out): The capability the kernel belongs to, or zero.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k)
- [`nk_missing_library_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_library_k) Header-only.

## nk_reduce_moments_*_best

Horizontal moments reduction (sum + sum-of-squares) over a strided array.

```c
nk_status_t nk_reduce_moments_f64_best(
    nk_f64_t const *data,
    nk_size_t count,
    nk_size_t stride,
    nk_f64_t *sum_ptr,
    nk_f64_t *sumsq_ptr,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_reduce_moments_f64_best`: data `nk_f64_t`, sum_ptr `nk_f64_t`, sumsq_ptr `nk_f64_t`
- `nk_reduce_moments_f32_best`: data `nk_f32_t`, sum_ptr `nk_f64_t`, sumsq_ptr `nk_f64_t`
- `nk_reduce_moments_i8_best`: data `nk_i8_t`, sum_ptr `nk_i64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_u8_best`: data `nk_u8_t`, sum_ptr `nk_u64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_i16_best`: data `nk_i16_t`, sum_ptr `nk_i64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_u16_best`: data `nk_u16_t`, sum_ptr `nk_u64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_i32_best`: data `nk_i32_t`, sum_ptr `nk_i64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_u32_best`: data `nk_u32_t`, sum_ptr `nk_u64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_i64_best`: data `nk_i64_t`, sum_ptr `nk_i64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_u64_best`: data `nk_u64_t`, sum_ptr `nk_u64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_f16_best`: data `nk_f16_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_bf16_best`: data `nk_bf16_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_e4m3_best`: data `nk_e4m3_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_e5m2_best`: data `nk_e5m2_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_e2m3_best`: data `nk_e2m3_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_e3m2_best`: data `nk_e3m2_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_e2m1_best`: data `nk_e2m1x2_t`, sum_ptr `nk_f32_t`, sumsq_ptr `nk_f32_t`
- `nk_reduce_moments_i4_best`: data `nk_i4x2_t`, sum_ptr `nk_i64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_u4_best`: data `nk_u4x2_t`, sum_ptr `nk_u64_t`, sumsq_ptr `nk_u64_t`
- `nk_reduce_moments_u1_best`: data `nk_u1x8_t`, sum_ptr `nk_u64_t`, sumsq_ptr `nk_u64_t`

- `nk_reduce_moments_f64` backends: serial, neon, haswell, skylake, rvv, v128, cuda, rocm
- `nk_reduce_moments_f32` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_i8` backends: serial, neon, neonsdot, haswell, skylake, icelake, sierra, rvv, v128, cuda, rocm
- `nk_reduce_moments_u8` backends: serial, neon, neonsdot, haswell, skylake, icelake, alder, sierra, rvv, v128, cuda, rocm
- `nk_reduce_moments_i16` backends: serial, neon, haswell, skylake, icelake, alder, rvv, v128, cuda, rocm
- `nk_reduce_moments_u16` backends: serial, neon, haswell, skylake, alder, rvv, v128, cuda, rocm
- `nk_reduce_moments_i32` backends: serial, neon, haswell, skylake, rvv, v128, cuda, rocm
- `nk_reduce_moments_u32` backends: serial, neon, haswell, skylake, rvv, v128, cuda, rocm
- `nk_reduce_moments_i64` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_u64` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_f16` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_bf16` backends: serial, neonbfdot, haswell, skylake, genoa, rvv, v128, cuda, rocm
- `nk_reduce_moments_e4m3` backends: serial, neon, neonfhm, haswell, skylake, genoa, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_e5m2` backends: serial, neon, neonfhm, haswell, skylake, genoa, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_e2m3` backends: serial, neon, neonsdot, haswell, skylake, icelake, alder, sierra, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_e3m2` backends: serial, neon, haswell, skylake, icelake, alder, rvv, v128relaxed, cuda, rocm
- `nk_reduce_moments_e2m1` backends: serial, cuda, rocm
- `nk_reduce_moments_i4` backends: serial, neonsdot, haswell, skylake, cuda, rocm
- `nk_reduce_moments_u4` backends: serial, neonsdot, haswell, skylake, cuda, rocm
- `nk_reduce_moments_u1` backends: serial, neon, haswell, skylake, cuda, rocm

- `data` (in): Pointer to the input data.
- `count` (in): Counts dimensions, a multiple of the values per byte.
- `stride` (in): Byte stride between elements, `sizeof(*data)` for contiguous arrays.
- `sum_ptr` (out): Output sum.
- `sumsq_ptr` (out): Output sum of squares.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_reduce_minmax_*_best

Horizontal min+max reduction with argmin/argmax over a strided array.

```c
nk_status_t nk_reduce_minmax_f64_best(
    nk_f64_t const *data,
    nk_size_t count,
    nk_size_t stride,
    nk_f64_t *min_value_ptr,
    nk_size_t *min_index_ptr,
    nk_f64_t *max_value_ptr,
    nk_size_t *max_index_ptr,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_reduce_minmax_f64_best`: data `nk_f64_t`, min_value_ptr, max_value_ptr `nk_f64_t`
- `nk_reduce_minmax_f32_best`: data `nk_f32_t`, min_value_ptr, max_value_ptr `nk_f32_t`
- `nk_reduce_minmax_i8_best`: data `nk_i8_t`, min_value_ptr, max_value_ptr `nk_i8_t`
- `nk_reduce_minmax_u8_best`: data `nk_u8_t`, min_value_ptr, max_value_ptr `nk_u8_t`
- `nk_reduce_minmax_i16_best`: data `nk_i16_t`, min_value_ptr, max_value_ptr `nk_i16_t`
- `nk_reduce_minmax_u16_best`: data `nk_u16_t`, min_value_ptr, max_value_ptr `nk_u16_t`
- `nk_reduce_minmax_i32_best`: data `nk_i32_t`, min_value_ptr, max_value_ptr `nk_i32_t`
- `nk_reduce_minmax_u32_best`: data `nk_u32_t`, min_value_ptr, max_value_ptr `nk_u32_t`
- `nk_reduce_minmax_i64_best`: data `nk_i64_t`, min_value_ptr, max_value_ptr `nk_i64_t`
- `nk_reduce_minmax_u64_best`: data `nk_u64_t`, min_value_ptr, max_value_ptr `nk_u64_t`
- `nk_reduce_minmax_f16_best`: data `nk_f16_t`, min_value_ptr, max_value_ptr `nk_f16_t`
- `nk_reduce_minmax_bf16_best`: data `nk_bf16_t`, min_value_ptr, max_value_ptr `nk_bf16_t`
- `nk_reduce_minmax_e4m3_best`: data `nk_e4m3_t`, min_value_ptr, max_value_ptr `nk_e4m3_t`
- `nk_reduce_minmax_e5m2_best`: data `nk_e5m2_t`, min_value_ptr, max_value_ptr `nk_e5m2_t`
- `nk_reduce_minmax_e2m3_best`: data `nk_e2m3_t`, min_value_ptr, max_value_ptr `nk_e2m3_t`
- `nk_reduce_minmax_e3m2_best`: data `nk_e3m2_t`, min_value_ptr, max_value_ptr `nk_e3m2_t`
- `nk_reduce_minmax_i4_best`: data `nk_i4x2_t`, min_value_ptr, max_value_ptr `nk_i8_t`
- `nk_reduce_minmax_u4_best`: data `nk_u4x2_t`, min_value_ptr, max_value_ptr `nk_u8_t`
- `nk_reduce_minmax_u1_best`: data `nk_u1x8_t`, min_value_ptr, max_value_ptr `nk_u8_t`

- `nk_reduce_minmax_f64` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_f32` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_i8` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_u8` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_i16` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_u16` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_i32` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_u32` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_i64` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_u64` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_f16` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_bf16` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_e4m3` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_e5m2` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_e2m3` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_e3m2` backends: serial, neon, haswell, skylake, rvv, v128relaxed, cuda, rocm
- `nk_reduce_minmax_i4` backends: serial, cuda, rocm
- `nk_reduce_minmax_u4` backends: serial, cuda, rocm
- `nk_reduce_minmax_u1` backends: serial, cuda, rocm

- `data` (in): Pointer to the input data.
- `count` (in): Counts dimensions, a multiple of the values per byte.
- `stride` (in): Byte stride between elements, `sizeof(*data)` for contiguous arrays.
- `min_value_ptr` (out): Output minimum value.
- `min_index_ptr` (out): Output index of the minimum value.
- `max_value_ptr` (out): Output maximum value.
- `max_index_ptr` (out): Output index of the maximum value.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_reduce_moments_*_dtype

Returns the accumulator dtype for the `sum` output of reduce_moments.

```c
nk_dtype_t nk_reduce_moments_sum_dtype(nk_dtype_t dtype);
```

- `nk_reduce_moments_sum_dtype`
- `nk_reduce_moments_sumsq_dtype`

Float types accumulate into wider floats; signed ints into i64; unsigned ints into u64.

Same as sum except all integers (signed and unsigned) accumulate into u64.

## nk_reduce_moments_error_bound

Returns the error bound of both reduce_moments outputs, per [`nk_accumulation_error_bound`](https://www.unum.cloud/docs/numkong/c/scalar-types.md#nk_accumulation_error_bound) of their accumulators.

```c
nk_f64_t nk_reduce_moments_error_bound(nk_dtype_t dtype);
```

## nk_reduce_minmax_value_dtype

Returns the value dtype for reduce_minmax outputs.

```c
nk_dtype_t nk_reduce_minmax_value_dtype(nk_dtype_t dtype);
```

Standard types return themselves. Sub-byte types widen: i4 → i8, u4 → u8, u1 → u8.

## nk_trig_sin_*_best

Element-wise sine over f64 inputs in radians.

```c
nk_status_t nk_trig_sin_f64_best(
    nk_f64_t const *ins,
    nk_size_t n,
    nk_f64_t *outs,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_trig_sin_f64_best`: ins, outs `nk_f64_t`
- `nk_trig_sin_f32_best`: ins, outs `nk_f32_t`
- `nk_trig_sin_f16_best`: ins, outs `nk_f16_t`

- `nk_trig_sin_f64` backends: serial, neon, haswell, skylake, v128relaxed, rvv
- `nk_trig_sin_f32` backends: serial, neon, haswell, skylake, v128relaxed, rvv
- `nk_trig_sin_f16` backends: serial, neonhalf, svehalf, skylake, sapphire, rvv

- `ins` (in): Input array of angles in radians.
- `n` (in): Number of elements in the input/output arrays.
- `outs` (out): Output array of sine values.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_trig_cos_*_best

Element-wise cosine over f64 inputs in radians.

```c
nk_status_t nk_trig_cos_f64_best(
    nk_f64_t const *ins,
    nk_size_t n,
    nk_f64_t *outs,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_trig_cos_f64_best`: ins, outs `nk_f64_t`
- `nk_trig_cos_f32_best`: ins, outs `nk_f32_t`
- `nk_trig_cos_f16_best`: ins, outs `nk_f16_t`

- `nk_trig_cos_f64` backends: serial, neon, haswell, skylake, v128relaxed, rvv
- `nk_trig_cos_f32` backends: serial, neon, haswell, skylake, v128relaxed, rvv
- `nk_trig_cos_f16` backends: serial, neonhalf, svehalf, skylake, sapphire, rvv

- `ins` (in): Input array of angles in radians.
- `n` (in): Number of elements in the input/output arrays.
- `outs` (out): Output array of cosine values.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## nk_trig_atan_*_best

Element-wise arc-tangent over f64 inputs.

```c
nk_status_t nk_trig_atan_f64_best(
    nk_f64_t const *ins,
    nk_size_t n,
    nk_f64_t *outs,
    nk_capability_t capabilities,
    nk_stream_t stream);
```

- `nk_trig_atan_f64_best`: ins, outs `nk_f64_t`
- `nk_trig_atan_f32_best`: ins, outs `nk_f32_t`
- `nk_trig_atan_f16_best`: ins, outs `nk_f16_t`

- `nk_trig_atan_f64` backends: serial, neon, haswell, skylake, v128relaxed, rvv
- `nk_trig_atan_f32` backends: serial, neon, haswell, skylake, v128relaxed, rvv
- `nk_trig_atan_f16` backends: serial, neonhalf, svehalf, skylake, sapphire, rvv

- `ins` (in): Input array of input values.
- `n` (in): Number of elements in the input/output arrays.
- `outs` (out): Output array of arc-tangent values.
- `capabilities` (in): One device's capabilities, like [`nk_cpu_capabilities_enabled`](https://www.unum.cloud/docs/numkong/c/runtime-dispatch.md#nk_cpu_capabilities_enabled) reports.
- `stream` (in): Null on the CPU, or the GPU stream of that device to queue on.

Statuses: - [`nk_missing_kernel_k`](https://www.unum.cloud/docs/numkong/c/quickstart.md#nk_missing_kernel_k) No capability in `capabilities` has it.

## Numerical Stability

Integer sum is elementwise a[i]+b[i], clamped to the type's range. Serial widens to i64 and clamps on store, while NEON uses hardware saturating adds, SQADD/UQADD. An f16/bf16/FP8 sum is promoted to f32, added, and truncated back, so double rounding is possible. Scale/blend/fma use float alpha/beta arithmetic, round to nearest with ties to even, then clamp. f32/f64 operations are native precision with no widening.

## Relevant x86 Instructions

FP16 conversions, VCVTPH2PS/VCVTPS2PH, serve f16 scale/sum/blend/fma operations, converting to f32 for arithmetic and back. The 6-7 cycle latency is amortized over vector-width elements. Saturating integer adds, VPADDSW/VPADDUSW, protect i16/u16 sums from overflow without branching. FMA, VFMADD231PS, is the workhorse for scale (alpha*x+beta) and blend (alpha*a+beta*b).

| Intrinsic | Instruction | Icelake | Genoa |
|:-|:-|:-|:-|
| _mm512_cvtph_ps | VCVTPH2PS (ZMM, YMM) | 7cy @ p0+p5 | 6cy @ p12+p23 |
| _mm512_cvtps_ph | VCVTPS2PH (YMM, ZMM, I8) | 7cy @ p0+p5 | 7cy @ p12+p23 |
| _mm256_adds_epi16 | VPADDSW (YMM, YMM, YMM) | 1cy @ p01 | n/a |
| _mm256_adds_epu16 | VPADDUSW (YMM, YMM, YMM) | 1cy @ p01 | n/a |
| _mm512_fpclass_ps_mask | VFPCLASSPS (K, ZMM, I8) | 3cy @ p5 | 5cy @ p01 |
| _mm256_fmadd_ps | VFMADD231PS (YMM, YMM, YMM) | 4cy @ p01 | 4cy @ p01 |

## Relevant ARM NEON/SVE Instructions

On ARM, i8/u8 and f16 scale/blend/fma widen to f32 intermediates, as f16 ones would round where the serial kernels do not. Saturating adds (SQADD/UQADD) handle integer overflow. FMLA provides fused multiply-add for floating-point scale/blend/fma.

| Intrinsic | Instruction | M1 Firestorm | Graviton 3 | Graviton 4 |
|:-|:-|:-|:-|:-|
| vfmaq_f32 | FMLA.S (vec) | 4cy @ V0123 | 4cy @ V0123 | 4cy @ V0123 |
| vqaddq_s16 | SQADD (vec) | 3cy @ V0123 | 2cy @ V0123 | 2cy @ V0123 |
| vqaddq_u16 | UQADD (vec) | 3cy @ V0123 | 2cy @ V0123 | 2cy @ V0123 |
| vcvtq_f32_s32 | SCVTF (vec) | 3cy @ V0123 | 3cy @ V01 | 3cy @ V01 |
| vcvtnq_s32_f32 | FCVTNS (vec) | 3cy @ V0123 | 3cy @ V01 | 3cy @ V01 |

## References

__See also.__

x86 intrinsics: [https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html](https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html)

__See also.__

Arm intrinsics: [https://developer.arm.com/architectures/instruction-sets/intrinsics/](https://developer.arm.com/architectures/instruction-sets/intrinsics/)

## Numerical stability

All accumulations are performed with stable techniques and __saturation__ in mind. Single-precision inputs are aggregated in double-precision. Double-precision inputs are handled with __Neumaier-like__ compensated summation schemes. Mini-floats are propagated to more hardware-friendly types. And integers are handled with proper saturation logic, as opposed to simple pairwise saturation, meaning that if several extremely large values are followed by equal negative values, the sum will be zero.

All MinMax scans are performed with respect to NaN values beyond simple total ordering. All positive and negative NaN values are masked out on the fly and can never be included in the output. For empty or NaN-only inputs, the returned argmin/argmax positions will be set to sentinel value `NUMKONG_SIZE_MAX`.

## Reduction Strategy

The key insight is that `_mm512_reduce_add_ps()` and similar intrinsics are actually serial operations, they don't parallelize the reduction across lanes. The correct approach is:

- Accumulate vertically in SIMD registers throughout the entire loop
- Perform a single horizontal reduction at the very end, reconstructing the lane positions

```
__m512 sum_f32x16 = _mm512_setzero_ps();
for (...) {
    __m512 data_f32x16 = _mm512_loadu_ps(data_pointer);
    sum_f32x16 = _mm512_add_ps(sum_f32x16, data_f32x16);
}
// Single horizontal reduce at the end only
nk_f32_t result = nk_reduce_add_f32x16_skylake_(sum_f32x16);
```

## Stride Handling Strategies

- stride == sizeof(scalar): Contiguous SIMD loads with masked tail
- Large stride with gather support: Use gather instructions (32/64-bit types)
- Otherwise: Serial fallback

## Argmin/Argmax Strategy

Single-pass algorithm tracking both value and index in SIMD registers:

```
__m512 min_f32x16 = _mm512_set1_ps(NUMKONG_F32_INF);
__m512i min_idx_i32x16 = _mm512_setzero_si512();
__m512i current_idx_i32x16 = _mm512_setr_epi32(0,1,2,3,...,15);
__m512i step_i32x16 = _mm512_set1_epi32(16);
for (...) {
    __m512 data_f32x16 = _mm512_loadu_ps(data_pointer);
    __mmask16 lt_mask = _mm512_cmp_ps_mask(data_f32x16, min_f32x16, _CMP_LT_OQ);
    min_f32x16 = _mm512_mask_mov_ps(min_f32x16, lt_mask, data_f32x16);
    min_idx_i32x16 = _mm512_mask_mov_epi32(min_idx_i32x16, lt_mask, current_idx_i32x16);
    current_idx_i32x16 = _mm512_add_epi32(current_idx_i32x16, step_i32x16);
}
```

## Accuracy

- `f64` sin, cos and atan: within 2 ULP of the correctly rounded result, so not faithful.
- `f32` sin and cos: within 2 ULP for |x| ≤ 10⁴, next to their zeros too; atan within 3 ULP.
- `f16:` shorter polynomials evaluated in `f32`, within 1 ULP of the correctly rounded result.
- `f16` on NEON FP16, SVE FP16 and Sapphire Rapids: evaluated in `f16`, also within 1 ULP.
- Their sine and cosine reduce in `f16` up to |x| ≤ 256, and in `f32` past it.
The measurements behind these bounds are in `include/numkong/trigonometry/README.md`.

## GLibC IEEE-754-compliant Math Functions

The GNU C Library, GLibC, provides IEEE-754-compliant math functions, like single-precision `sinf` and `cosf` or double-precision `sin` and `cos`. Those are accurate to ~0.55 ULP, units in the last place, but can be slow to evaluate. They combine techniques like:

- Taylor series expansions for small values.
- Table lookups combined with corrections for moderate values.
- Accurate modulo reduction for large values.
The precomputed tables may be the hardest part to accelerate with SIMD, as they contain 440x values, each 64-bit wide.

__See also.__

glibc argument reduction: [https://github.com/lattera/glibc/blob/895ef79e04a953cac1493863bcae29ad85657ee1/sysdeps/ieee754/dbl-64/branred.c#L54](https://github.com/lattera/glibc/blob/895ef79e04a953cac1493863bcae29ad85657ee1/sysdeps/ieee754/dbl-64/branred.c#L54)

__See also.__

glibc sine: [https://github.com/lattera/glibc/blob/895ef79e04a953cac1493863bcae29ad85657ee1/sysdeps/ieee754/dbl-64/s_sin.c#L84](https://github.com/lattera/glibc/blob/895ef79e04a953cac1493863bcae29ad85657ee1/sysdeps/ieee754/dbl-64/s_sin.c#L84)

## Approximation Algorithms

There are several ways to approximate trigonometric functions, and the choice depends on the target hardware and the desired precision.

Taylor Series approximation is a series expansion of a sum of its derivatives at a target point. It's easy to derive for differentiable functions and works well for functions smooth around the expansion point, but can perform poorly for functions with singularities or rapid, high-frequency oscillations.

Pade approximations are rational functions that approximate a function by a ratio of polynomials. They often converge faster than Taylor for functions with singularities or steep changes, and approximate both smooth and rational functions well, but can be more computationally intensive to evaluate, and can have holes, undefined points.

Moreover, most approximations can be combined with Horner's methods of evaluating polynomials to reduce the number of multiplications and additions, and to improve the numerical stability. In trigonometry, the Payne-Hanek Range Reduction is another technique used to reduce the argument to a smaller range, where the approximation is more accurate.

## Optimization Notes

The following optimizations were evaluated but did not yield performance improvements.

Estrin's scheme for polynomial evaluation is a tree-based approach that reduces the dependency depth from N sequential FMAs to log2(N), by computing powers of x in parallel with partial sums. For an 8-term polynomial, it reduces depth from 7 to 3, but benchmarks showed a ~20% regression: the extra MUL operations for computing x², x⁴, x⁸ hurt throughput more than the reduced dependency depth helps latency. For large arrays, out-of-order execution across loop iterations already hides FMA latency, making throughput the bottleneck.

RCPPS with Newton-Raphson refinement, a fast reciprocal approximation at ~4 cycles with one refinement iteration for ~22-bit precision, was tested as an alternative to VDIVPS at ~11 cycles. It did not improve performance when combined with Estrin's scheme, likely because the division is not on the critical path when processing large arrays.

## Relevant x86 Instructions

Polynomial evaluation (Horner's method) for sin/cos/tan uses chained FMAs - the 4-cycle latency is hidden by out-of-order execution across iterations. Range reduction uses VRNDSCALE for fast rounding (notably 3x faster on Genoa than Ice Lake). VFPCLASS detects NaN/Inf inputs for special case handling. Division appears in tangent's final step but isn't on the critical path.

| Intrinsic | Instruction | Icelake | Genoa |
|:-|:-|:-|:-|
| _mm512_roundscale_ps | VRNDSCALEPS (ZMM, ZMM, I8) | 8cy @ p0+p0 | 3cy @ p23 |
| _mm512_roundscale_pd | VRNDSCALEPD (ZMM, ZMM, I8) | 8cy @ p0+p0 | 3cy @ p23 |
| _mm512_fpclass_ps_mask | VFPCLASSPS (K, ZMM, I8) | 3cy @ p5 | 5cy @ p01 |
| _mm512_fmadd_ps | VFMADD231PS (ZMM, ZMM, ZMM) | 4cy @ p0 | 4cy @ p01 |
| _mm256_fmadd_ps | VFMADD231PS (YMM, YMM, YMM) | 4cy @ p01 | 4cy @ p01 |
| _mm256_div_ps | VDIVPS (YMM, YMM, YMM) | ~11cy @ p0 | ~11cy @ p01 |
| _mm256_div_pd | VDIVPD (YMM, YMM, YMM) | ~13cy @ p0 | ~13cy @ p01 |

## Relevant ARM NEON/SVE Instructions

ARM implementations use the same Horner polynomial approach with FMLA chains. FRINTA provides fast rounding for range reduction. The 4-cycle FMA latency with 4 inst/cycle throughput allows excellent pipelining when processing multiple elements.

| Intrinsic | Instruction | M1 Firestorm | Graviton 3 | Graviton 4 |
|:-|:-|:-|:-|:-|
| vfmaq_f32 | FMLA.S (vec) | 4cy @ V0123 | 4cy @ V0123 | 4cy @ V0123 |
| vfmaq_f64 | FMLA.D (vec) | 4cy @ V0123 | 4cy @ V0123 | 4cy @ V0123 |
| vrndaq_f32 | FRINTA.S | 2cy @ V0123 | 2cy @ V01 | 2cy @ V01 |

## References

__See also.__

x86 intrinsics: [https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html](https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html)

__See also.__

Arm intrinsics: [https://developer.arm.com/architectures/instruction-sets/intrinsics/](https://developer.arm.com/architectures/instruction-sets/intrinsics/)

## Overview

Contains following element-wise operations:

- Scale (Multiply) with shift: result[i] = α · a[i] + β
- Sum (Add): result[i] = a[i] + b[i]
- Blend: result[i] = α · a[i] + β · b[i]
- FMA (Fused Multiply-Add): result[i] = α · a[i] · b[i] + β · c[i]Beyond their obvious usecases, those can be reused for vector-scalar math and other operations:

- Scale with β = 0 for a pure multiply.
- Sum is equivalent to WSum with α = β = 1.
- Average is WSum with α = β = 0.5.
- Elementwise multiply is FMA with β = 0.For dtypes:

- f64c: 64-bit complex × 64-bit complex scales
- f32c: 32-bit complex × 32-bit complex scales
- f64: 64-bit IEEE floating point × 64-bit scales
- f32: 32-bit IEEE floating point × 32-bit scales
- f16: 16-bit IEEE floating point × 32-bit scales
- bf16: 16-bit brain floating point × 32-bit scales
- e4m3: 8-bit e4m3 floating point × 32-bit scales
- e5m2: 8-bit e5m2 floating point × 32-bit scales
- e2m3: 8-bit e2m3 floating point (MX) × 32-bit scales
- e3m2: 8-bit e3m2 floating point (MX) × 32-bit scales
- i8/u8: 8-bit integers × 32-bit scales
- i16/u16: 16-bit integers × 32-bit scales
- i32/u32: 32-bit integers × 64-bit scales
- i64/u64: 64-bit integers × 64-bit scalesFor hardware architectures:

- Arm: NEON, NEON+F16, NEON+BF16
- x86: Haswell, Skylake, Ice Lake, Sapphire Rapids
- RISC-V: RVVProvides horizontal reduction operations over vectors with:

- `nk_reduce_moments_*` —sum + sum-of-squares in one pass
- `nk_reduce_minmax_*` —min + max with argmin/argmax in one pass
- Runtime dispatch for runtime ISA selectionFor dtypes:

- f64: 64-bit IEEE floating point numbers
- f32: 32-bit IEEE floating point numbers
- f16: 16-bit IEEE floating point numbers
- bf16: 16-bit brain floating point numbers
- e4m3: 8-bit e4m3 floating point numbers
- e5m2: 8-bit e5m2 floating point numbers
- e2m3: 8-bit e2m3 floating point numbers (MX)
- e3m2: 8-bit e3m2 floating point numbers (MX)
- i8: 8-bit signed integers
- u8: 8-bit unsigned integers
- i16: 16-bit signed integers
- u16: 16-bit unsigned integers
- i32: 32-bit signed integers
- u32: 32-bit unsigned integers
- i64: 64-bit signed integers
- u64: 64-bit unsigned integers
- i4: 4-bit signed integers (packed pairs)
- u4: 4-bit unsigned integers (packed pairs)
- u1: 1-bit binary (packed octets)For hardware architectures:

- Arm: NEON, NEON+F16, NEON+FHM, NEON+BF16, NEON+SDOT
- x86: Haswell, Skylake, Ice Lake, Genoa, Sierra Forest
- RISC-V: RVV
- WASM: V128, V128RelaxedContains:

- Sine and Cosine approximations: fast for `f32` vs accurate for `f64`
- Tangent and the 2-argument arctangent: fast for `f32` vs accurate for `f64`For dtypes:

- 64-bit IEEE-754 floating point
- 32-bit IEEE-754 floating point
- 16-bit IEEE-754 floating pointFor hardware architectures:

- Arm: NEON, NEON FP16, SVE FP16
- x86: Haswell, Skylake, Sapphire RapidsThose functions partially complement the `each.h` module, and are necessary for the `geospatial.h` module, among others. Both Haversine and Vincenty's formulas require trigonometric functions, and those are the most expensive part of the computation.

__See also.__

SLEEF: [https://sleef.org/](https://sleef.org/)

