fastCD - v1.3.1

fastCD is a lean C kernel for Cayley–Dickson multiplication over real and complex numbers and finite fields, with a NumPy-compatible Python extension, written for researchers and engineers in need of faster Cayley–Dickson multiplication. The key advantage of fastCD is that it implements a new O(N log N) algorithm which, to the author's knowledge, has never been presented before. See the accompanying paper for more information on the mathematical foundations and rigorous derivation of correctness and time complexity. fastCD may not be completely suitable for every application, but it is licensed under the MIT License so that it may be improved upon and adapted to the user's needs completely, with no questions asked.

Besides the quasilinear algorithm, fastCD includes direct Cayley–Dickson multiplication and the algorithm of Cariow and Cariowa. All three are available for real coefficients, while complex and finite-field coefficients use direct and quasilinear. Quasilinear has the best asymptotic complexity of the three, but it is not the fastest at every low dimension, so the default algorithm for each dimension was chosen from single-core benchmarks of single and batched products. Over the real numbers, direct is the default for single products at dimension 4, and direct or Cariow–Cariowa is the default at dimensions 8 and 16; see Benchmark.

In an exact-reference accuracy study, the median floating-point error of the quasilinear product on pooled Gaussian and uniform inputs is within a factor of two of direct multiplication at every tested dimension; see Numerical stability.

The canonical source repository is git.tanuki-cd.com/algebraity/fastCD.

Convention

A_0 := F,   A_n := A_(n-1) x A_(n-1),   (a,b)(c,d) := (ac - conjugate(d) b, d a + b conjugate(c))

The coefficient field F is the real numbers (the default), the complex numbers, or an odd-characteristic finite field F_(p^d). An element is a contiguous array of N := 2^n coefficients in the recursive standard basis: coordinate zero is the scalar part, and the two halves of the array are the pair (a, b). Each coefficient is one double over R, two doubles with the real part first over C (complex128 in NumPy), and d uint32 residues over F_(p^d). Conventions and results gives the full conventions and explains the algorithms: the basis signs, the quadratic form, how the quasilinear recursion works, and the exact operation counts.

Using the C library

The C library can be used by including the fastcd.h header file and using fastCD functions within a normal C program; <fastcd/complex.h> and <fastcd/finite.h> declare the complex and finite-field arithmetic. Here is an example of supported usage over the real numbers.

#include <stdlib.h>
#include <fastcd/fastcd.h>

int main(void) {
    const size_t N = 16;
    double x[16] = {1.0, 2.0};
    double y[16] = {3.0, 4.0};
    double product[16];
    size_t length;

    if (fastcd_mul_workspace_size(N, &length) != FASTCD_SUCCESS) return 1;
    double* scratch = length == 0 ? NULL : malloc(length * sizeof(*scratch));
    if (length != 0 && scratch == NULL) return 1;

    FastCDWorkspaceF64 workspace = {scratch, length};
    FastCDStatus status = fastcd_mul_f64(product, x, y, N,
            length == 0 ? NULL : &workspace);
    free(scratch);
    return status == FASTCD_SUCCESS ? 0 : 1;
}

For more information on how to get started using fastCD in C and to familiarize yourself with the API, take a look at C quick start and C API reference.

Using Python

fastCD can be used in Python via the package pyfastcd. It requires CPython 3.10 or newer and NumPy 1.26 or newer. To install it from PyPI with pip, run:

pip install pyfastcd

The package can also be built from source from this repository. This also requires a C11 compiler, Python development headers, and CMake 3.18 or newer. To build and install it from the repository root, run:

python -m pip install .

After importing pyfastcd, you can use the Algebra, FiniteField, Element, and ElementArray classes. Here's an example of what you can do over the real numbers; Algebra(dimension=8, char=13, deg=5) or Algebra(dimension=8, char=0, deg=2) selects octonions over F_(13^5) or over the complex numbers.

import numpy as np
import pyfastcd as fcd

O = fcd.Algebra(dimension=8)
a = O([1, 2, 3, 4, 5, 6, 7, 8])
b = O.basis(2)

product = a * b
square = a ** 2
commutator = a.commutator(b)

x = O.asarray(np.ones((1000, 8)))
y = x * b                                  # batch product, leading axes broadcast
z = O.multiply(x, b, algorithm="quasilinear")

An Element is one immutable Cayley–Dickson element. An ElementArray holds a batch of elements in a NumPy coefficient array of shape (..., N) (or (..., N, d) over an extension field) and evaluates batch operations in compiled code. See the Python quick start and the Python API reference for more information on how each class works.

Building the C library

A C11 compiler and CMake 3.18 or newer are required to build. Unix-like systems link the system math library. From the repository root, run:

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --parallel
cmake -E chdir build ctest --output-on-failure
cmake --install build --prefix /desired/prefix

The default configuration builds the static library, the test suite, and the benchmark runner together. -DBUILD_SHARED_LIBS=ON builds a shared library, and -DBUILD_TESTING=OFF -DFASTCD_BUILD_BENCHMARK=OFF builds the library alone. The installation provides the <fastcd/*.h> headers, libfastcd, the license, and a CMake package: use find_package(fastcd CONFIG REQUIRED) and link fastcd::fastcd.

Benchmark

All timings below were measured on 2026-11-09 on a Dell PowerEdge T440 with an Intel Xeon Gold 6148, using one logical CPU per run, GCC 16.2.1 and Linux 7.2.11. The host is shared, and no governor, turbo, or SMT settings were changed.

Only the C kernel is timed, with storage already allocated: Python overhead, allocation, and verification are outside the timed region. Every product of every operand pair was also checked against an independent implementation.

The measured source is commit 33addc9, built from a clean checkout; the archived build records list the SHA-256 of every measured source file. The raw observations, environment, commands, and build records are archived in benchmark-results/t440-20261109/, and the benchmark documentation describes the protocol and how to rerun it.

Real coefficients

Mean nanoseconds per product with two-sided 95% Student-t intervals over 24 timing blocks of 10,000 operand pairs, in batch mode (one call for all 10,000 products, two products per SIMD instruction where the kernels support it) and single mode (one validated call per product). All 600,000 verification checks passed with zero mismatches. fastest mean names the lowest mean, or a tie when the paired per-block ratio of the two lowest means does not separate them, and default is the algorithm that automatic selection uses. The automatic choices agree with every cell of this run; single products at dimension 2 are a tie between direct and quasilinear, and the default there is quasilinear.

Batch products (ns per product)

Single products (ns per product)

Finite fields and complex coefficients

The same protocol compares direct and quasilinear multiplication over C and over the nine finite fields in the table below, which are the field benchmark's default set: the prime fields F_3, F_5, F_13, F_65521 and F_2147483647 and the extension fields F_(3^5), F_(5^2), F_(13^5) and F_(65521^2). Each cell uses 1,000 operand pairs and 24 timing blocks. Every pair is verified before timing, against independent schoolbook polynomial arithmetic over the finite fields and against a long-double reference over C; all 100,000 pairs passed. Each cell below gives the quasilinear mean in nanoseconds per single product and, in parentheses, the paired direct/quasilinear ratio, so a ratio above 1 favors quasilinear. Direct is faster at dimension 2 everywhere and at dimension 4 over C and five of the nine fields. From dimension 8 upward quasilinear is faster in every cell except two ties: F_(3^5) at dimension 64 and single products over F_(65521^2) at dimension 16. Automatic selection uses direct for N <= 4 and quasilinear otherwise, over C and every finite field. From dimension 16 upward, batch and single means differ by at most 6%, except over C at dimension 16, where batching saves 10%. Each cell runs in one process, so its interval covers variation between timing blocks but not between processes: in repeated runs on this host, individual cells from dimension 64 upward differed by up to 17%, and direct/quasilinear ratios by up to 8%. Both modes, every dimension and the confidence intervals are in domains.md in the archive.

Single products over finite fields and C: quasilinear ns per product (direct/quasilinear)

Comparison with FLINT

FLINT (Fast Library for Number Theory) is an independent exact-arithmetic C library. The comparison times fastCD's finite-field arithmetic against FLINT 3.6.0 on two workloads with identical operands and defining polynomials, 1,000 operand pairs per cell and 24 paired blocks, every pair verified against the same independent implementation; all 132,000 pairs passed. The scalar workload is one field multiplication (fastcd_ff_field_mul against nmod_mul for prime fields or fq_nmod_mul for extensions). The CD workload is one complete product at each dimension, where both sides evaluate the same quasilinear schedule and only the coefficient arithmetic differs: fastCD's packed kernels against FLINT's field operations. Ratios are paired geometric means of FLINT time over fastCD time, so a ratio above 1 favors fastCD. As in the finite-field comparison, each cell runs in one process, so its interval does not include variation between processes or runs. Every dimension and the intervals for the CD cells are in flint.md in the archive.

Scalar field multiplication: ns per product

Complete CD multiplication: fastCD ns / FLINT ns (FLINT/fastCD)

Numerical stability

Every floating-point sum inside the quasilinear product whose length grows with the dimension is accumulated in blocks of eight with a balanced tree between blocks, and above dimension 256 the Cariow–Cariowa correction sums are combined by a balanced tree. Accuracy was measured against an exact reference on 13,584 operand pairs from 21 real and 19 complex input families. The error of a computed product fl(xy) is reported as eta = ||fl(xy) - xy||_2 / (u ||x||_2 ||y||_2), with u = 2^-53 the binary64 unit roundoff. On pooled Gaussian and uniform inputs, the median eta of the quasilinear product is within a factor of two of direct multiplication at every tested dimension:

The largest error observed anywhere in the corpus is 3.7 over R and 4.1 over C (direct: 2.9 and 2.2; sequential accumulation on the same corpus: 78 and 18). These figures describe the measured corpus, not a worst case. From dimension 64 up, balanced accumulation costs at most 6% more time on default real products, and 5–13% on complex products and on builds without the explicit SIMD kernels. The method, the input families and the archived summaries and charts are described in Numerical stability. The raw operands and observations are kept on the author's machine and are not included in Git or the source distributions; reruns on another NumPy version can draw different random inputs.

Documentation

• C quick start and C API reference
• Python quick start and Python API reference
• Conventions and results
• Benchmark
• Numerical stability

License and Attribution

All code is released under the MIT License and is presented with no warranty or guarantees. You may use it in any way for any purpose and are not restricted by the original author of the code. To cite the software, use the metadata in CITATION.cff. The changelog is in CHANGELOG.md.

The quasilinear algorithm is due to Harrison Lemley (me!), Quasilinear multiplication in the Cayley–Dickson tower. The Cariow–Cariowa algorithm follows Aleksandr Cariow and Galina Cariowa, An unified approach for developing rationalized algorithms for hypercomplex number multiplication, Przegląd Elektrotechniczny 91(2), 36–39 (2015).