machine_vectors.h – SIMD-accelerated operations on fixed-length vectors

Vector types and operations mapping onto the target’s SIMD instructions, together with the flint_sgemm and flint_dgemm matrix multiplication kernels built on top of them.

Backends are selected automatically: AVX2 (and AVX512 for the vec8dz family) on x86, NEON on ARM, and otherwise a generic backend, using GNU vector extensions where the compiler supports them and plain ISO C structs elsewhere. The AVX2 and NEON backends require a 64-bit word size, since their integer vectors have ulong lanes; a 32-bit build uses the generic backends, which provide only the floating-point types. The generic backends implement only the subset of the interface required by flint_sgemm/flint_dgemm; the full interface, as used by fft_small, still requires AVX2 or NEON.

For the vector operations to use the target’s instructions, FLINT must be built with appropriate compiler flags. configure chooses these from the detected CPU; --enable-avx2 and --enable-avx512 force them on.

Defining FLINT_MACHINE_VECTORS_FORCE_GENERIC before including this header selects the generic backend even on a target with AVX2 or NEON, and FLINT_MACHINE_VECTORS_STRICT_C additionally selects the ISO C tier over GNU vector extensions. These are intended for testing and profiling the portable code paths.

The strict ISO C tier has no way to express a vector operation, so whether its operations become SIMD instructions is entirely up to the compiler’s SLP vectorizer; its performance therefore varies between compilers and compiler versions, unlike the other backends. It exists so that the header works on compilers without GNU vector extensions, and is not the fallback used on GCC or clang.

The generic backends express a fused multiply-add as a * b + c, which a compiler fuses into an FMA instruction only when floating-point contraction is enabled. GCC in a strict ISO mode, which is how FLINT is built, does not contract by default; code using these operations in a performance-critical loop should request contraction, as machine_vectors/gemm.c does with #pragma GCC optimize ("fp-contract=fast"). The AVX2, AVX512 and NEON backends use fused intrinsics and are unaffected.

Some functions may require that vectors are aligned in memory.

Types

type vec1n

type vec2n

type vec4n

type vec8n

Vector with 1, 2, 4, or 8 ulong entries.

type vec1d

type vec2d

type vec4d

type vec8d

Vector with 1, 2, 4, or 8 double entries.

type vec1f

type vec4f

type vec8f

type vec16f

Vector with 1, 4, 8, or 16 float entries.

type vec8dz

type vec16dz

type vec16fz

type vec32fz

type vec8nz

Vectors backed by AVX512 registers, available only when building with AVX512F support: 8 or 16 double entries, 16 or 32 float entries, and 8 ulong entries respectively. The z suffix distinguishes these from the equally-named types built from pairs of narrower registers; vec8d, for instance, remains a pair of AVX2 registers, which gives more instruction level parallelism in existing code.

Printing

void vec4d_print(vec4d a)

void vec4n_print(vec4n a)

Access and conversions

vec1d vec1d_load(const double *a)

vec4d vec4d_load(const double *a)

vec8d vec8d_load(const double *a)

vec1d vec1d_load_aligned(const double *a)

vec4d vec4d_load_aligned(const double *a)

vec8d vec8d_load_aligned(const double *a)

vec1d vec1d_load_unaligned(const double *a)

vec4d vec4d_load_unaligned(const double *a)

vec8d vec8d_load_unaligned(const double *a)

vec4n vec4n_load_unaligned(const ulong *a)

vec8n vec8n_load_unaligned(const ulong *a)

void vec1d_store(double *z, vec1d a)

void vec4d_store(double *z, vec4d a)

void vec8d_store(double *z, vec8d a)

void vec1d_store_aligned(double *z, vec1d a)

void vec4d_store_aligned(double *z, vec4d a)

void vec8d_store_aligned(double *z, vec8d a)

void vec1d_store_unaligned(double *z, vec1d a)

void vec4d_store_unaligned(double *z, vec4d a)

void vec4n_store_unaligned(ulong *z, vec4n a)

void vec8d_store_unaligned(double *z, vec8d a)

double vec4d_get_index(vec4d a, const int i)

double vec8d_get_index(vec8d a, int i)

Extract the entry at index \(i\).

vec1d vec1d_set_d(double a)

vec4d vec4d_set_d(double a)

vec4n vec4n_set_n(ulong a)

vec8d vec8d_set_d(double a)

vec8n vec8n_set_n(ulong a)

Set all entries to the same value.

vec4d vec4d_set_d4(double a0, double a1, double a2, double a3)

vec4n vec4n_set_n4(ulong a0, ulong a1, ulong a2, ulong a3)

vec8d vec8d_set_d8(double a0, double a1, double a2, double a3, double a4, double a5, double a6, double a7)

Create vector from distinct entries.

vec1n vec1d_convert_limited_vec1n(vec1d a)

vec2n vec2d_convert_limited_vec2n(vec2d a)

vec4n vec4d_convert_limited_vec4n(vec4d a)

Given that each entry in the input vector is an exact integer in [0, 2^{52}), convert it to ulong. Note that vec1d and vec2d functions are only available on NEON.

vec2d vec2n_convert_limited_vec2d(vec2n a)

vec4d vec4n_convert_limited_vec4d(vec4n a)

vec8d vec8n_convert_limited_vec8d(vec8n a)

The inverse conversion, from ulong entries to exact double entries. Requires the same assumption.

Permutations

vec4d vec4d_unpacklo(vec4d a, vec4d b)

vec4d vec4d_unpackhi(vec4d a, vec4d b)

vec4d vec4d_permute_0_2_1_3(vec4d a)

vec4d vec4d_permute_3_1_2_0(vec4d a)

vec4d vec4d_permute_3_2_1_0(vec4d a)

vec4d vec4d_permute2_0_2(vec4d a, vec4d b)

vec4d vec4d_permute2_1_3(vec4d a, vec4d b)

vec4d vec4d_unpack_lo_permute_0_2_1_3(vec4d u, vec4d v)

vec4d vec4d_unpack_hi_permute_0_2_1_3(vec4d u, vec4d v)

vec4d vec4d_unpackhi_permute_3_1_2_0(vec4d u, vec4d v)

vec4d vec4d_unpacklo_permute_3_1_2_0(vec4d u, vec4d v)

VEC4D_TRANSPOSE(z0, z1, z2, z3, a0, a1, a2, a3)

Sets the rows z to the transpose of the 4x4 matrix given by rows a.

Comparisons

int vec1d_same(double a, double b)

int vec4d_same(vec4d a, vec4d b)

int vec8d_same(vec8d a, vec8d b)

Check whether the vectors are equal.

vec4d vec4d_cmp_ge(vec4d a, vec4d b)

vec4d vec4d_cmp_gt(vec4d a, vec4d b)

Entrywise comparisons.

Arithmetic and basic operations

vec1d vec1d_round(vec1d a)

vec4d vec4d_round(vec4d a)

vec8d vec8d_round(vec8d a)

vec1d vec1d_zero()

vec4d vec4d_zero()

vec8d vec8d_zero()

vec1d vec1d_one()

vec4d vec4d_one()

vec8d vec8d_one()

vec1d vec1d_add(vec1d a, vec1d b)

vec1d vec1d_sub(vec1d a, vec1d b)

vec4d vec4d_add(vec4d a, vec4d b)

vec4d vec4d_sub(vec4d a, vec4d b)

vec4n vec4n_add(vec4n a, vec4n b)

vec4n vec4n_sub(vec4n a, vec4n b)

vec8d vec8d_add(vec8d a, vec8d b)

vec8d vec8d_sub(vec8d a, vec8d b)

vec1d vec1d_addsub(vec1d a, vec1d b)

vec4d vec4d_addsub(vec4d a, vec4d b)

vec1d vec1d_neg(vec1d a)

vec4d vec4d_neg(vec4d a)

vec8d vec8d_neg(vec8d a)

vec1d vec1d_abs(vec1d a)

vec4d vec4d_abs(vec4d a)

vec1d vec1d_max(vec1d a, vec1d b)

vec1d vec1d_min(vec1d a, vec1d b)

vec4d vec4d_max(vec4d a, vec4d b)

vec4d vec4d_min(vec4d a, vec4d b)

vec8d vec8d_max(vec8d a, vec8d b)

vec8d vec8d_min(vec8d a, vec8d b)

vec1d vec1d_mul(vec1d a, vec1d b)

vec4d vec4d_mul(vec4d a, vec4d b)

vec8d vec8d_mul(vec8d a, vec8d b)

vec1d vec1d_half(vec1d a)

vec4d vec4d_half(vec4d a)

vec1d vec1d_div(vec1d a, vec1d b)

vec4d vec4d_div(vec4d a, vec4d b)

vec8d vec8d_div(vec8d a, vec8d b)

vec1d vec1d_fmadd(vec1d a, vec1d b, vec1d c)

vec4d vec4d_fmadd(vec4d a, vec4d b, vec4d c)

vec8d vec8d_fmadd(vec8d a, vec8d b, vec8d c)

vec1d vec1d_fmsub(vec1d a, vec1d b, vec1d c)

vec4d vec4d_fmsub(vec4d a, vec4d b, vec4d c)

vec8d vec8d_fmsub(vec8d a, vec8d b, vec8d c)

vec1d vec1d_fnmadd(vec1d a, vec1d b, vec1d c)

vec4d vec4d_fnmadd(vec4d a, vec4d b, vec4d c)

vec8d vec8d_fnmadd(vec8d a, vec8d b, vec8d c)

vec1d vec1d_fnmsub(vec1d a, vec1d b, vec1d c)

vec4d vec4d_fnmsub(vec4d a, vec4d b, vec4d c)

vec8d vec8d_fnmsub(vec8d a, vec8d b, vec8d c)

vec1d vec1d_blendv(vec1d a, vec1d b, vec1d c)

vec4d vec4d_blendv(vec4d a, vec4d b, vec4d c)

vec8d vec8d_blendv(vec8d a, vec8d b, vec8d c)

vec4n vec4n_bit_shift_right(vec4n a, ulong b)

vec8n vec8n_bit_shift_right(vec8n a, ulong b)

vec4n vec4n_bit_and(vec4n a, vec4n b)

vec8n vec8n_bit_and(vec8n a, vec8n b)

Modular arithmetic

These functions are used internally by the small-prime FFT. Some double variants assume an odd modulus \(n < 2^{50}\). Other assumptions are not yet documented.

int vec1d_same_mod(vec1d a, vec1d b, vec1d n, vec1d ninv)

int vec4d_same_mod(vec4d a, vec4d b, vec4d n, vec4d ninv)

Return whether \(a\) and \(b\) are the same mod \(n\).

vec1d vec1d_reduce_pm1no_to_0n(vec1d a, vec1d n)

vec1d vec4d_reduce_pm1no_to_0n(vec4d a, vec4d n)

vec8d vec8d_reduce_pm1no_to_0n(vec8d a, vec8d n)

Return \(a \bmod n\) reduced to \([0,n)\) assuming \(a \in (-n,n)\).

vec1d vec1d_reduce_to_pm1n(vec1d a, vec1d n, vec1d ninv)

vec4d vec4d_reduce_to_pm1n(vec4d a, vec4d n, vec4d ninv)

vec8d vec8d_reduce_to_pm1n(vec8d a, vec8d n, vec8d ninv)

Return \(a \bmod n\) reduced to \([-n,n]\).

vec1d vec1d_reduce_to_pm1no(vec1d a, vec1d n, vec1d ninv)

vec4d vec4d_reduce_to_pm1no(vec4d a, vec4d n, vec4d ninv)

vec8d vec8d_reduce_to_pm1no(vec8d a, vec8d n, vec8d ninv)

Return \(a \bmod n\) reduced to \((-n,n)\).

vec1d vec1d_reduce_0n_to_pmhn(vec1d a, vec1d n)

vec4d vec4d_reduce_0n_to_pmhn(vec4d a, vec4d n)

Return \(a \bmod n\) reduced to \([-n/2, n/2]\) given \(a \in [0,n]\).

vec1d vec1d_reduce_pm1n_to_pmhn(vec1d a, vec1d n)

vec4d vec4d_reduce_pm1n_to_pmhn(vec4d a, vec4d n)

vec8d vec8d_reduce_pm1n_to_pmhn(vec8d a, vec8d n)

Return \(a \bmod n\) reduced to \([-n/2, n/2]\) given \(a \in [-n,n]\).

vec1d vec1d_reduce_2n_to_n(vec1d a, vec1d n)

vec4d vec4d_reduce_2n_to_n(vec4d a, vec4d n)

vec8d vec8d_reduce_2n_to_n(vec8d a, vec8d n)

Return \(a \bmod n\) reduced to \([0,n)\) given \(a \in [0,2n)\).

vec1d vec1d_reduce_to_0n(vec1d a, vec1d n, vec1d ninv)

vec4d vec4d_reduce_to_0n(vec4d a, vec4d n, vec4d ninv)

vec8d vec8d_reduce_to_0n(vec8d a, vec8d n, vec8d ninv)

Return \(a \bmod n\) reduced to \([0,n)\).

vec1d vec1d_mulmod(vec1d a, vec1d b, vec1d n, vec1d ninv)

vec4d vec4d_mulmod(vec4d a, vec4d b, vec4d n, vec4d ninv)

vec8d vec8d_mulmod(vec8d a, vec8d b, vec8d n, vec8d ninv)

Return \(ab \bmod n\) in \([-n,n]\) with assumptions.

vec1d vec1d_nmulmod(vec1d a, vec1d b, vec1d n, vec1d ninv)

vec4d vec4d_nmulmod(vec4d a, vec4d b, vec4d n, vec4d ninv)

vec8d vec8d_nmulmod(vec8d a, vec8d b, vec8d n, vec8d ninv)

Return \(ab \bmod n\) in \([-n,n]\) with assumptions.

vec4n vec4n_addmod(vec4n a, vec4n b, vec4n n)

vec8n vec8n_addmod(vec8n a, vec8n b, vec8n n)

Return \(a + b \bmod n\) in \([0,n)\)

vec4n vec4n_addmod_limited(vec4n a, vec4n b, vec4n n)

vec8n vec8n_addmod_limited(vec8n a, vec8n b, vec8n n)

Return \(a + b \bmod n\) in \([0,n)\), assuming that \(n < 2^{63}\).

Matrix multiplication

These functions compute a matrix product in single or double precision. They are always available: when FLINT is built with BLAS, the default is to call it, and otherwise FLINT’s own kernels are used. The intended use is as a building block for exact linear algebra over \(\mathbb{Z}\) and \(\mathbb{Z}/n\mathbb{Z}\), for example in nmod_mat_mul_blas() and fmpz_mat_mul_blas().

All of these functions compute \(C = AB\) for row-major matrices, where C is m by n, A is m by k and B is k by n, with ldc, lda and ldb the respective leading dimensions (the number of entries between the start of consecutive rows, which must be at least the number of columns). No transposition or accumulation is performed: the previous contents of C are overwritten, and \(k = 0\) sets C to zero. This is equivalent to cblas_sgemm or cblas_dgemm called with CblasRowMajor, CblasNoTrans, CblasNoTrans, alpha equal to 1 and beta equal to 0. Aliasing of C with A or B is not allowed.

The FLINT kernels are multithreaded internally according to flint_get_num_threads(), using FLINT’s thread pool. They handle arbitrary dimensions, including thin and unbalanced shapes, without requiring any padding or alignment of the input.

void flint_sgemm(slong m, slong k, slong n, const float *A, slong lda, const float *B, slong ldb, float *C, slong ldc)

void flint_dgemm(slong m, slong k, slong n, const double *A, slong lda, const double *B, slong ldb, double *C, slong ldc)

Sets \(C = AB\), calling either the BLAS or the FLINT implementation according to flint_gemm_use_blas.

void flint_sgemm_blas(slong m, slong k, slong n, const float *A, slong lda, const float *B, slong ldb, float *C, slong ldc)

void flint_dgemm_blas(slong m, slong k, slong n, const double *A, slong lda, const double *B, slong ldb, double *C, slong ldc)

Sets \(C = AB\) using cblas_sgemm or cblas_dgemm. These raise an exception if FLINT was built without BLAS support.

void flint_sgemm_fallback(slong m, slong k, slong n, const float *A, slong lda, const float *B, slong ldb, float *C, slong ldc)

void flint_dgemm_fallback(slong m, slong k, slong n, const double *A, slong lda, const double *B, slong ldb, double *C, slong ldc)

Sets \(C = AB\) using FLINT’s own kernels, which are always available.

int flint_gemm_use_blas

Selects the implementation used by flint_sgemm() and flint_dgemm(). It is initialized to 1 if FLINT was built with BLAS support and 0 otherwise, and may be set to either value at runtime, for example to compare the two implementations. Setting it to 1 in a build without BLAS support will result in an exception when a gemm is attempted.