iii. Convolution Variants

Convolutions beyond the 1D sum rule: 2D, digit-wise in any base (generalised XOR / OR / AND), arbitrary-length DFT and evaluation at a geometric progression (chirp-z / Bluestein), $(\min ,+)$ with convexity, and arithmetic in ${\mathbb{F}}_{2^{64}}$. Uses ntt, convolution, Poly, polyMul from NTT.

Variant	Method	Cost
2D: $c_{ij}=\sum a_{xy}{b}_{i-x,j-y}$	flatten rows with stride $M=m_1+m_2-1$, one 1D product	$O(NM\log NM)$
digit-wise sum mod $k$ (base-$k$ “XOR”)	$k$-point DFT per digit in ${\mathbb{Z}}_p[x]/(x^k-1)$, reduce by ${\Phi }_k$	$O(N{\log }_{k}N\cdot k^2)$
digit-wise max / min (base-$k$ OR / AND)	prefix / suffix sums per digit (zeta), multiply, differences	$O(N{\log }_{k}N)$
DFT of any length; evaluation at $1,z,z^2,\dots$	chirp-z: $jk=\left(\genfrac{}{}{0ex}{}{j+k}{2}\right)-\left(\genfrac{}{}{0ex}{}{j}{2}\right)-\left(\genfrac{}{}{0ex}{}{k}{2}\right)$, one convolution	$O((n+m)\log )$
$(\min ,+)$, both convex	merge slopes	$O(n+m)$
$(\min ,+)$, one convex	monotone minima (divide and conquer)	$O((n+m)\log )$
$(\min ,+)$, general	no truly subquadratic algorithm known	$O(nm)$
polynomials over ${\mathbb{F}}_2$ / ${\mathbb{F}}_{2^{64}}$	carry-less multiplication, GF64	$O(64)$ per product

Grid pattern matching with wildcards uses three 2D convolutions of $\sum (a-b{)}^{2}ab$ terms, exactly like 1D wildcard matching (Other Matching); an arbitrary modulus replaces convolution with convMod (FFT).

vector<vector<mint>> convolution2D(const vector<vector<mint>>& a, const vector<vector<mint>>& b) {
    if (a.empty() || b.empty() || a[0].empty() || b[0].empty()) return {};
    int n1 = sz(a), m1 = sz(a[0]), n2 = sz(b), m2 = sz(b[0]), N = n1 + n2 - 1, M = m1 + m2 - 1;
    Poly fa(n1 * M), fb(n2 * M);
    for (int i = 0; i < n1; i++) for (int j = 0; j < m1; j++) fa[i * M + j] = a[i][j];
    for (int i = 0; i < n2; i++) for (int j = 0; j < m2; j++) fb[i * M + j] = b[i][j];
    Poly fc = polyMul(fa, fb); fc.resize(N * M);
    vector<vector<mint>> c(N, vector<mint>(M));
    for (int i = 0; i < N; i++) for (int j = 0; j < M; j++) c[i][j] = fc[i * M + j];
    return c;
}
 
Poly chirpZ(const Poly& a, mint z, int m) {               // a(z^k) for k = 0..m-1, O((n+m) log)
    int n = sz(a); if (n == 0 || m == 0) return Poly(m, 0);
    auto tri = [&](ll t) { return z.pow(t * (t - 1) / 2 % 998244352); };   // z^{C(t,2)} (exponent mod p-1, z != 0)
    if (z == mint(0)) { Poly r(m, a[0]); if (m) { r[0] = 0; for (auto x : a) r[0] += x; } return r; }
    Poly A(n), B(n + m - 1);
    for (int j = 0; j < n; j++) A[j] = a[j] / tri(j);        // a_j z^{-C(j,2)}
    for (int t = 0; t < n + m - 1; t++) B[t] = tri(t);       // z^{C(t,2)}
    reverse(all(A));
    Poly C = polyMul(A, B), r(m);
    for (int k = 0; k < m; k++) r[k] = C[n - 1 + k] / tri(k);   // sum_j a_j z^{C(j+k,2) - C(j,2) - C(k,2)}
    return r;
}

Proof: 2D flattening and the chirp identity

2D. Entry $(i,j)$ is placed at $iM+j$ with $j<M$. A product term $a_{xy}{b}_{uv}$ lands at $(x+u)M+(y+v)$, and $y+v\le m_1+m_2-2<M$, so no term crosses into the next row: the 1D product read with stride $M$ is the 2D convolution. Chirp-z. $\left(\genfrac{}{}{0ex}{}{j+k}{2}\right)=\left(\genfrac{}{}{0ex}{}{j}{2}\right)+\left(\genfrac{}{}{0ex}{}{k}{2}\right)+jk$, so $z^{jk}=z^{\left(\genfrac{}{}{0ex}{}{j+k}{2}\right)}{z}^{-\left(\genfrac{}{}{0ex}{}{j}{2}\right)}{z}^{-\left(\genfrac{}{}{0ex}{}{k}{2}\right)}$ and $a(z^k)=z^{-\left(\genfrac{}{}{0ex}{}{k}{2}\right)}{\sum }_j{a}_j{z}^{-\left(\genfrac{}{}{0ex}{}{j}{2}\right)}{z}^{\left(\genfrac{}{}{0ex}{}{j+k}{2}\right)}$: with $A$ reversed, ${\sum }_j{A}_{n-1-j}^{\prime }{B}_{j+k}$ is coefficient $n-1+k$ of the product. This uses only integer exponents, so no square root of $z$ is needed (the classical Bluestein form $jk=\frac{j^2+k^2-(j-k{)}^2}{2}$ needs ${\omega }^{1/2}$, i.e. a $2n$-th root, over $\mathbb{C}$). With $z$ a primitive $n$-th root of unity it is a DFT of arbitrary length $n$. Exponents are reduced mod $p-1$ (Fermat). Checked against direct evaluation (and $z=0$).

DFT length	Use
power of two	Cooley–Tukey directly
smooth (small prime factors)	mixed-radix Cooley–Tukey
$n=n_1{n}_2$ coprime	Good–Thomas
prime	Rader (cyclic convolution over $(\mathbb{Z}/n{)}^{\ast }$) or Bluestein
arbitrary, or a partial spectrum / zoom FFT	Bluestein (chirp-z): about 3–6× a native transform

In contests you rarely need a DFT of a fixed non-power-of-two length (zero-padding a linear convolution is simpler); chirp-z earns its place for cyclic convolution of exact length $n$ when wraparound matters, and for evaluating at a geometric progression $z^0,z^1,\dots$ (“evaluate at powers of $q$” problems, $q$-analogues), which arbitrary multipoint evaluation (Evaluation and Interpolation) does only in $O(n{\log }^{2}n)$.

Digit-wise convolutions

For $N=k^d$ indices written in base $k$: $c_{x\oplus y}=\sum a_x{b}_y$ where $\oplus$ adds digits mod $k$ without carry ($k=2$: XOR), and the digit-wise max / min versions ($k=2$: OR / AND). The base-$k$ DFT needs a primitive $k$-th root of unity; instead of requiring $k\mid p-1$, compute in ${\mathbb{Z}}_p[x]/(x^k-1)$ with $w=x$ (multiplying by $w^e$ is a cyclic shift) and reduce the final values modulo the cyclotomic polynomial ${\Phi }_k$, where $x$ is a genuine primitive $k$-th root; each result is a constant there. Examples: CodeChef MDSWIN (ternary), XORTREEH (base $k$ with powers).

template <class Z> struct DigitRing {                      // Z[x] / (x^k - 1), reduced modulo Phi_k at the end
    int k; vector<ll> phi;                                  // Phi_k, monic, integer coefficients
    explicit DigitRing(int base) : k(base) {
        auto cyclo = [&](auto&& self, int n) -> vector<ll> {
            vector<ll> p(n + 1, 0); p[0] = -1, p[n] = 1;    // x^n - 1, divided by Phi_d for d | n, d < n
            for (int d = 1; d < n; d++) if (n % d == 0) {
                auto q = self(self, d); int dq = sz(q) - 1; vector<ll> res(sz(p) - dq, 0);
                for (int i = sz(p) - 1; i >= dq; i--) { ll c = p[i]; res[i - dq] = c; for (int j = 0; j <= dq; j++) p[i - dq + j] -= c * q[j]; }
                p = res;
            }
            return p;
        };
        phi = cyclo(cyclo, k);
    }
    vector<Z> mul(const vector<Z>& a, const vector<Z>& b) const {
        vector<Z> c(k, Z(0)); for (int i = 0; i < k; i++) if (a[i] != Z(0)) for (int j = 0; j < k; j++) c[(i + j) % k] += a[i] * b[j]; return c;
    }
    Z reduce(vector<Z> a) const {
        int d = sz(phi) - 1;
        for (int i = k - 1; i >= d; i--) { Z c = a[i]; if (c == Z(0)) continue; for (int j = 0; j <= d; j++) a[i - d + j] -= c * Z(phi[j]); }
        return a[0];
    }
    void transform(vector<vector<Z>>& a, bool inverse) const {
        int n = sz(a); vector<vector<Z>> tmp(k, vector<Z>(k));
        for (int len = 1; len < n; len *= k) for (int i = 0; i < n; i += len * k) for (int j = 0; j < len; j++) {
            for (int p = 0; p < k; p++) {
                fill(all(tmp[p]), Z(0));
                for (int q = 0; q < k; q++) {
                    int e = p * q % k; if (inverse) e = (k - e) % k;
                    const auto& src = a[i + j + q * len];
                    for (int c = 0; c < k; c++) tmp[p][(c + e) % k] += src[c];
                }
            }
            for (int p = 0; p < k; p++) a[i + j + p * len] = tmp[p];
        }
    }
};
template <class Z> vector<Z> digitAddPower(vector<Z> a, ll e, int k, const vector<Z>& b = {}) {   // a^e (e >= 1), or a * b if b given
    int n = 1; while (n < max(sz(a), sz(b))) n *= k;
    DigitRing<Z> ring(k);
    auto lift = [&](const vector<Z>& v) { vector<vector<Z>> r(n, vector<Z>(k, Z(0))); for (int i = 0; i < sz(v); i++) r[i][0] = v[i]; ring.transform(r, false); return r; };
    auto A = lift(a);
    if (!b.empty()) { auto B = lift(b); for (int i = 0; i < n; i++) A[i] = ring.mul(A[i], B[i]); }
    else for (int i = 0; i < n; i++) { vector<Z> res(k, Z(0)), base = A[i]; res[0] = 1; for (ll t = e; t > 0; t >>= 1, base = ring.mul(base, base)) if (t & 1) res = ring.mul(res, base); A[i] = res; }
    ring.transform(A, true);
    Z invN = Z(1) / Z(n); vector<Z> c(n);
    for (int i = 0; i < n; i++) c[i] = ring.reduce(A[i]) * invN;
    return c;
}
template <class Z> vector<Z> digitAddConvolution(const vector<Z>& a, const vector<Z>& b, int k) { return digitAddPower(a, 1, k, b); }
template <class Z> vector<Z> digitMaxConvolution(vector<Z> a, vector<Z> b, int k, bool useMin = false) {   // digit-wise max (or min)
    int n = 1; while (n < max(sz(a), sz(b))) n *= k;
    a.resize(n, Z(0)), b.resize(n, Z(0));
    auto zeta = [&](vector<Z>& f, bool inverse) {
        for (int len = 1; len < n; len *= k) for (int i = 0; i < n; i += len * k) for (int j = 0; j < len; j++) {
            if (!useMin) { if (!inverse) for (int p = 1; p < k; p++) f[i + j + p * len] += f[i + j + (p - 1) * len];
                           else for (int p = k - 1; p > 0; p--) f[i + j + p * len] -= f[i + j + (p - 1) * len]; }
            else { if (!inverse) for (int p = k - 2; p >= 0; p--) f[i + j + p * len] += f[i + j + (p + 1) * len];
                   else for (int p = 0; p + 1 < k; p++) f[i + j + p * len] -= f[i + j + (p + 1) * len]; }
        }
    };
    zeta(a, false), zeta(b, false);
    for (int i = 0; i < n; i++) a[i] *= b[i];
    zeta(a, true); return a;
}
template <class Z> vector<Z> digitMinConvolution(const vector<Z>& a, const vector<Z>& b, int k) { return digitMaxConvolution(a, b, k, true); }

Proof: digit-wise transforms

Digit-wise sum. The group $({\mathbb{Z}}_k{)}^d$ has characters ${\chi }_y(x)=w^{⟨x,y⟩}$ for a primitive $k$-th root of unity $w$; the DFT ${\hat{a}}_y={\sum }_x{a}_x{w}^{⟨x,y⟩}$ turns the group convolution into a pointwise product, and the inverse uses $w^{-1}$ and division by $N$ (orthogonality ${\sum }_y{w}^{⟨x,y⟩}=N[x=0]$). The transform factorises digit by digit, which is the loop over len. In $R={\mathbb{Z}}_p[x]/(x^k-1)$ we use $w=x$; all identities used except orthogonality hold in $R$, and the final values, mapped to ${\mathbb{Z}}_p[x]/{\Phi }_k$ where $x$ is a primitive $k$-th root of unity, satisfy orthogonality there (${\sum }_{e<k}{x}^{ej}=\frac{x^{kj}-1}{x^j-1}=0$ for $k\nmid j$, since ${\Phi }_k\mid x^k-1$ but ${\Phi }_k$ and $x^j-1$ are coprime modulo $p$ when $p\nmid k$). So the reduction of each result modulo ${\Phi }_k$ is the true integer convolution value, a constant; dividing by $N$ needs $p\nmid k$. Powers $a^{\oplus e}$ exponentiate each transformed coordinate. Digit-wise max: prefix sums along every digit give ${\hat{a}}_y={\sum }_{x\le y}{a}_x$ (digit-wise $\le$); $\max (x,x^{\prime })\le y$ iff $x\le y$ and $x^{\prime }\le y$, so $\hat{c}=\hat{a}\hat{b}$, and taking differences along each digit inverts the prefix sums (min symmetrically with suffix sums). Checked against brute force for $k=2,3,5,6$.

(min, +) convolution

$c_k={\min }_{i+j=k}{a}_i+b_j$. Both convex: $c$ is convex and its slope sequence is the merged sorted slopes of $a$ and $b$ (a Minkowski sum, Minkowski Sum). $a$ convex, $b$ arbitrary: the optimal $j$ for row $k$ is monotone in $k$, so divide and conquer on rows with shrinking column ranges (monotone minima / SMAWK) works. Max-plus with concave inputs: negate. Typical use: merging knapsack-like DP layers where one side is convex (items of one kind sorted by value), and $(\min ,+)$ products inside Aliens’ trick (Aliens Trick).

template <class T> vector<T> minPlusConvexConvex(const vector<T>& a, const vector<T>& b) {
    int n = sz(a), m = sz(b); if (n == 0 || m == 0) return {};
    vector<T> c(n + m - 1); c[0] = a[0] + b[0];
    for (int k = 1, i = 0, j = 0; k < n + m - 1; k++) {
        if (j == m - 1 || (i < n - 1 && a[i + 1] - a[i] < b[j + 1] - b[j])) i++; else j++;
        c[k] = a[i] + b[j];
    }
    return c;
}
template <class T> vector<T> minPlusConvexArbitrary(const vector<T>& a, const vector<T>& b) {   // a convex
    int n = sz(a), m = sz(b); if (n == 0 || m == 0) return {};
    vector<T> c(n + m - 1);
    auto rec = [&](auto&& self, int kl, int kr, int jl, int jr) -> void {   // rows [kl, kr) have an optimal j in [jl, jr]
        if (kl >= kr) return;
        int k = (kl + kr) / 2, lo = max(jl, k - n + 1), hi = min(jr, k), best = lo;
        for (int j = lo + 1; j <= hi; j++) if (a[k - j] + b[j] < a[k - best] + b[best]) best = j;
        c[k] = a[k - best] + b[best];
        self(self, kl, k, jl, best); self(self, k + 1, kr, best, jr);
    };
    rec(rec, 0, n + m - 1, 0, m - 1);
    return c;
}

Proof: the two (min, +) algorithms

Convex + convex. A pair $(i,j)$ is reached by taking $k$ unit steps from $(0,0)$, each increasing $i$ or $j$; the cost increases by the slope taken. Since both slope sequences are non-decreasing, taking the smaller next slope first is optimal for every prefix length (exchange argument: any other order of the same multiset of steps has a prefix sum that is at least as large), so $c_k$ is $a_0+b_0$ plus the $k$ smallest available slopes in merge order, which the loop computes. Convex + arbitrary. For $k<k^{\prime }$ and $j<j^{\prime }$ (both feasible), convexity of $a$ gives $a_{k-j}+a_{k^{\prime }-j^{\prime }}\le a_{k-j^{\prime }}+a_{k^{\prime }-j}$ when the indices are in range (the four indices satisfy $k-j^{\prime }<k-j,k^{\prime }-j^{\prime }<k^{\prime }-j$ with equal sums), i.e. $M_{kj}=a_{k-j}+b_j$ is Monge on its feasible band. Hence a leftmost-tie optimal column of row $k$ is non-decreasing in $k$ (if row $k^{\prime }$ preferred $j<j_k^{\ast }$ strictly, adding the Monge inequality would contradict the optimality of $j_k^{\ast }$ in row $k$), and restricting rows left of $k$ to $[jl,best]$ and rows right of it to $[best,jr]$ loses nothing; the band constraint $j\in [k-n+1,k]$ is itself monotone in $k$. Each recursion level scans $O(n+m)$ columns: $O((n+m)\log (n+m))$. Checked against the $O(nm)$ definition on random convex / arbitrary inputs.

GF(2^64)

Elements are 64-bit masks = polynomials over ${\mathbb{F}}_2$ modulo the irreducible $x^{64}+x^4+x^3+x+1$: addition is XOR, multiplication is carry-less multiplication then reduction. A field of characteristic 2 ($x+x=0$), so random evaluation of polynomials whose coefficients only matter mod 2 works (Schwartz–Zippel with error $\le \mathrm{deg}/2^{64}$), e.g. “is the number of perfect matchings odd” (CodeChef EVENPOLY): substitute random field elements into the Tutte / Edmonds matrix and test $\det \ne 0$. Nim multiplication (Sprague–Grundy) is another field of size $2^{64}$.

using u64 = unsigned long long; using u128 = unsigned __int128;
u128 clmul(u64 a, u64 b) { u128 r = 0; for (int i = 0; i < 64; i++) if (b >> i & 1) r ^= (u128)a << i; return r; }   // carry-less product
struct GF64 {
    u64 v = 0;
    GF64() {} explicit GF64(u64 x) : v(x) {}
    static u64 reduce(u128 x) {                              // mod x^64 + x^4 + x^3 + x + 1
        u64 hi = x >> 64, lo = (u64)x;
        u128 t = (u128)hi ^ ((u128)hi << 1) ^ ((u128)hi << 3) ^ ((u128)hi << 4);   // hi x^64 = hi (x^4 + x^3 + x + 1)
        u64 hi2 = t >> 64; lo ^= (u64)t;
        lo ^= hi2 ^ (hi2 << 1) ^ (hi2 << 3) ^ (hi2 << 4);   // hi2 has <= 4 bits, no further overflow
        return lo;
    }
    friend GF64 operator+(GF64 a, GF64 b) { return GF64(a.v ^ b.v); }
    friend GF64 operator-(GF64 a, GF64 b) { return GF64(a.v ^ b.v); }
    friend GF64 operator*(GF64 a, GF64 b) { return GF64(reduce(clmul(a.v, b.v))); }
    GF64& operator+=(GF64 o) { return *this = *this + o; }
    GF64& operator*=(GF64 o) { return *this = *this * o; }
    friend bool operator==(GF64 a, GF64 b) { return a.v == b.v; }
    GF64 pow(u64 e) const { GF64 r(1), b = *this; for (; e; e >>= 1, b *= b) if (e & 1) r *= b; return r; }
    GF64 inv() const { return pow(~0ULL - 1); }              // a^(2^64 - 2)
    GF64 operator/(GF64 o) const { return *this * o.inv(); }
};

Proof: GF(2^64) arithmetic

${\mathbb{F}}_2[x]/(f)$ is a field iff $f$ is irreducible; $x^{64}+x^4+x^3+x+1$ is the standard irreducible (indeed primitive) pentanomial-type modulus of degree 64, so the 64-bit masks form ${\mathbb{F}}_{2^{64}}$. Carry-less multiplication is polynomial multiplication over ${\mathbb{F}}_2$ (XOR adds coefficients). Reduction: $x^{64}\equiv x^4+x^3+x+1$, so the high word $h$ contributes $h\cdot (x^4+x^3+x+1)$, a polynomial of degree $\le 67$; its bits above 63 ($\le 4$ of them) are reduced once more, and $h_2(x^4+x^3+x+1)$ has degree $\le 7$, so it fits. Inverse: the multiplicative group has order $2^{64}-1$, so $a^{2^{64}-2}=a^{-1}$ (Lagrange). Verified: $a\cdot a^{-1}=1$, distributivity and associativity on random elements, and reduce against bitwise long division.

Related notes

• Bitwise Convolutions (XOR / OR / AND, subset convolution)
• Evaluation and Interpolation
• Aliens Trick