fft-dihedral

Fast Fourier transforms for the dihedral group over prime finite fields.

This package computes the Fourier transform of functions on

D_{2n} = <r, s | r^n = s^2 = e, srs^{-1} = r^{-1}>.

The fast path reduces the nonabelian dihedral transform to four cyclic number-theoretic transforms (NTTs). The inverse transform reconstructs two cyclic spectra and applies two inverse NTTs.

Install

From PyPI:

python -m pip install fft-dihedral

From a checkout:

python -m pip install -e ".[plot,dev]"

Quick Start

from fft_dihedral import (
    DEFAULT_MODULUS,
    dihedral_dft_fft,
    dihedral_inverse_fft,
    flatten_transform,
    root_of_unity,
)

n = 16
omega = root_of_unity(n, DEFAULT_MODULUS)
rotations = list(range(n))
reflections = [2 * k for k in range(n)]

transform = dihedral_dft_fft(rotations, reflections, DEFAULT_MODULUS, omega)
assert len(flatten_transform(transform)) == 2 * n
assert dihedral_inverse_fft(transform, omega) == (rotations, reflections)

Conventions

A function f: D_{2n} -> GF(p) is represented by two length-n arrays:

rotations[k]   = f(r^k)
reflections[k] = f(s r^k)

For a primitive nth root of unity omega, the two-dimensional irreducible representations satisfy

rho_j(r^k)   = [[omega^(j k), 0], [0, omega^(-j k)]]
rho_j(s r^k) = [[0, omega^(-j k)], [omega^(j k), 0]]

With the normalized finite-group DFT convention

fhat(rho) = (1 / |D_{2n}|) sum_g f(g) rho(g),

each two-dimensional Fourier coefficient is assembled from four cyclic transforms:

fhat(rho_j) = 1/(2n) * [
  [DFT_omega(rotations)[j],    DFT_omega^{-1}(reflections)[j]],
  [DFT_omega(reflections)[j],  DFT_omega^{-1}(rotations)[j]]
]

The one-dimensional representations are computed by direct character sums.

What The Four NTTs Do

The input vector is really two vectors:

rotations   = [f(1), f(r), ..., f(r^{n-1})]
reflections = [f(s), f(sr), ..., f(sr^{n-1})]

The FFT applies cyclic NTTs to those two vectors with both orientations of the root:

A^+ = NTT_omega(rotations)
A^- = NTT_omega^{-1}(rotations)
B^+ = NTT_omega(reflections)
B^- = NTT_omega^{-1}(reflections)

For each two-dimensional irrep rho_j, the matrix coefficient is assembled as:

1/(2n) * [[A^+[j], B^-[j]],
          [B^+[j], A^-[j]]]

The inverse FFT reverses this packing. The one-dimensional coefficients recover the j = 0 cyclic frequencies, and, for even n, the j = n/2 frequencies. The two-dimensional matrices recover the paired frequencies j and n-j. Once the full A^+ and B^+ spectra have been rebuilt, two inverse NTTs give back the rotation and reflection coefficient vectors.

Group Algebra Multiplication

Elements of the group algebra are represented by the same two arrays:

x = sum_k rotations[k] r^k + sum_k reflections[k] s r^k

The package provides both a quadratic reference product and an FFT product:

from fft_dihedral import (
    DEFAULT_MODULUS,
    dihedral_multiply_fft,
    dihedral_multiply_naive,
    root_of_unity,
)

n = 16
omega = root_of_unity(n, DEFAULT_MODULUS)
lhs_rotations = [1] * n
lhs_reflections = [2] * n
rhs_rotations = [3] * n
rhs_reflections = [4] * n

fast = dihedral_multiply_fft(
    lhs_rotations,
    lhs_reflections,
    rhs_rotations,
    rhs_reflections,
    DEFAULT_MODULUS,
    omega,
)
slow = dihedral_multiply_naive(
    lhs_rotations,
    lhs_reflections,
    rhs_rotations,
    rhs_reflections,
    DEFAULT_MODULUS,
)
assert fast == slow

Because the transform is normalized by 1 / |D_{2n}|, multiplication in Fourier space is

(x y)^hat(rho) = |D_{2n}| xhat(rho) yhat(rho).

So the FFT product computes two transforms, multiplies each scalar or matrix block with the extra factor 2n, and then applies the inverse dihedral FFT.

Group Algebra Inversion

The same Fourier decomposition also gives fast inversion over compatible prime fields. In coefficient form, inverting an element means solving a dense 2n x 2n linear system. In Fourier form, the group algebra decomposes into scalar blocks and 2 x 2 matrix blocks, so inversion is blockwise:

from fft_dihedral import (
    DEFAULT_MODULUS,
    dihedral_invert_fft,
    dihedral_multiply_fft,
    root_of_unity,
)

n = 16
omega = root_of_unity(n, DEFAULT_MODULUS)
rotations = [0] * n
reflections = [0] * n
rotations[1] = 1

inv_rotations, inv_reflections = dihedral_invert_fft(
    rotations,
    reflections,
    DEFAULT_MODULUS,
    omega,
)

identity_rotations = [0] * n
identity_rotations[0] = 1
identity_reflections = [0] * n

assert dihedral_multiply_fft(
    rotations,
    reflections,
    inv_rotations,
    inv_reflections,
    DEFAULT_MODULUS,
    omega,
) == (identity_rotations, identity_reflections)

Because the transform is normalized by 1 / |D_{2n}|, Fourier-space inversion uses

(f^-1)^hat(rho) = |D_{2n}|^-2 * fhat(rho)^-1.

The package exposes this lower-level operation as invert_fourier_transform. Singular scalar or matrix blocks raise ValueError.

Supported Fields

The implementation supports prime fields GF(p) when:

• n is a power of two for the fast transform.
• p is prime.
• n | p - 1, so GF(p) contains a primitive nth root of unity.
• gcd(2n, p) = 1, so the normalization by 1/(2n) is defined.

The default field is GF(2013265921), where

2013265921 = 15 * 2^27 + 1

so it supports power-of-two rotation orders up to 2^27.

True extension fields GF(p^k) are not implemented yet. This is intentionally stricter than integer quotient rings like Z/p^kZ; those are not fields when k > 1.

Arbitrary Compatible Primes

Roots are selected automatically using Sympy:

from fft_dihedral import root_of_unity

omega = root_of_unity(16, 97)

CLI verification over a non-default prime:

fft-dihedral verify --n 16 --modulus 97

You can also pass an explicit root:

fft-dihedral verify --n 16 --modulus 97 --omega 8

Verify

python -m pytest
fft-dihedral verify --n 128
fft-dihedral verify --n 16 --modulus 97

Benchmark

fft-dihedral bench --min-exp 4 --max-exp 15 --repetitions 5 --naive-limit 512
fft-dihedral bench --min-exp 4 --max-exp 6 --modulus 257 --repetitions 3
fft-dihedral bench-mul --min-exp 4 --max-exp 12 --repetitions 5 --naive-limit 512
fft-dihedral bench-inv --min-exp 4 --max-exp 12 --repetitions 5

The benchmark prints FFT time / (2n log2(2n)). For an O(2n log(2n)) algorithm this ratio should stay roughly flat as n grows, modulo Python timing noise and cache effects.

The bench-mul command compares group-algebra multiplication by FFT against the direct quadratic product. The fast path computes two forward dihedral FFTs, multiplies the Fourier blocks, and applies one inverse dihedral FFT.

Plot

Install the plotting extra:

python -m pip install "fft-dihedral[plot]"

Then generate PNG and SVG timing plots:

fft-dihedral-plot --min-exp 4 --max-exp 14 --repetitions 3 --naive-limit 1024

Limitations

• No true extension-field GF(p^k) implementation yet.
• The radix-2 NTT is implemented in pure Python, so the Rust crate is much faster for production workloads.