i. Fast Fourier Transform

Evaluate at the $n$-th roots of unity, multiply pointwise, interpolate back: $O(n^2)$ convolution becomes $O(n\log n)$. This note has the complex FFT (real coefficients, exact integer products by rounding, arbitrary modulus by splitting), the convolution toolkit, Karatsuba and big-number multiplication. Exact modular convolution is the NTT, which most later notes use.

$\text{DFT}:X_k={\sum }_{j<n}{x}_j{\omega }^{jk},\omega =e^{2\pi i/n};\text{inverse}:x_j=\frac{1}{n}{\sum }_k{X}_k{\omega }^{-jk};c=a\ast b⟺C_k=A_k{B}_k(\text{cyclic, length }n).$

evaluate   A and B at the n-th roots of unity    O(n log n)
multiply   pointwise                             O(n)
interpolate (inverse transform)                  O(n log n)

	complex FFT	NTT (NTT)
arithmetic	double	integers mod $p=c{2}^k+1$
exactness	approximate, round at the end	exact
limits	about 53 bits of precision in the result	length $\le 2^k$; needs a modulus (or CRT)
speed	comparable (complex mult = 4 real)	slightly slower per butterfly with %
use	real data, exact integers with small products, arbitrary modulus by splitting	everything modulo 998244353

FFT variant	Idea
radix-2 DIT (below)	split by index parity; bit-reversal on input
radix-2 DIF	split by output parity; bit-reversal on output
radix-4 / split-radix	fewer multiplications, about 20% faster
Good–Thomas	$n=n_1{n}_2$ coprime: no twiddle factors (Theoretical Algorithms)
Rader	prime $n$ via a cyclic convolution over the multiplicative group
Bluestein (chirp-z)	any $n$, and evaluation at a geometric progression (Convolution Variants)

using cd = complex<double>;
void fft(vector<cd>& a, bool invert = false) {              // n a power of two; roots built in long double
    int n = sz(a); if (n <= 1) return;
    static vector<complex<long double>> R(2, 1); static vector<cd> rt(2, 1);   // rt[k + j] = e^{i pi j / k}
    for (static int k = 2; k < n; k *= 2) {
        R.resize(n); rt.resize(n); auto x = polar(1.0L, acosl(-1) / k);
        for (int i = k; i < 2 * k; i++) rt[i] = R[i] = i & 1 ? R[i / 2] * x : R[i / 2];
    }
    for (int i = 1, j = 0; i < n; i++) { int bit = n >> 1; for (; j & bit; bit >>= 1) j ^= bit; j ^= bit; if (i < j) swap(a[i], a[j]); }
    for (int k = 1; k < n; k *= 2)
        for (int i = 0; i < n; i += 2 * k) for (int j = 0; j < k; j++) {
            cd z = rt[j + k] * a[i + j + k];
            a[i + j + k] = a[i + j] - z; a[i + j] += z;
        }
    if (invert) { reverse(a.begin() + 1, a.end()); for (auto& x : a) x /= n; }
}
 
vector<double> convReal(const vector<double>& a, const vector<double>& b) {   // one forward + one inverse FFT (packing)
    if (a.empty() || b.empty()) return {};
    int need = sz(a) + sz(b) - 1, n = 1; while (n < need) n *= 2;
    vector<cd> f(n); for (int i = 0; i < sz(a); i++) f[i].real(a[i]); for (int i = 0; i < sz(b); i++) f[i].imag(b[i]);
    fft(f);
    vector<cd> p(n);
    for (int k = 0; k < n; k++) {                           // A_k = (F_k + conj F_{-k}) / 2, B_k = (F_k - conj F_{-k}) / 2i
        cd x = f[k], y = conj(f[(n - k) & (n - 1)]);
        p[k] = (x + y) * (x - y) / cd(0, 4);
    }
    fft(p, true);
    vector<double> res(need); for (int i = 0; i < need; i++) res[i] = p[i].real();
    return res;
}
 
vector<ll> convExact(const vector<ll>& a, const vector<ll>& b) {   // exact while every result (and sum |a||b|) stays below ~1e15
    vector<double> x(all(a)), y(all(b)); auto r = convReal(x, y);
    vector<ll> res(sz(r)); for (int i = 0; i < sz(r); i++) res[i] = llround(r[i]);
    return res;
}
 
vector<ll> convMod(const vector<ll>& a, const vector<ll>& b, ll M) {   // any M < 2^30, entries in [0, M); 4 FFTs
    if (a.empty() || b.empty()) return {};
    int need = sz(a) + sz(b) - 1, n = 1; while (n < need) n *= 2;
    const int cut = 1 << 15;
    vector<cd> P(n), Q(n);
    for (int i = 0; i < sz(a); i++) P[i] = cd(a[i] / cut, a[i] % cut);   // hi + i lo
    for (int i = 0; i < sz(b); i++) Q[i] = cd(b[i] / cut, b[i] % cut);
    fft(P); fft(Q);
    vector<cd> X(n), Y(n);
    for (int k = 0; k < n; k++) {
        int r = (n - k) & (n - 1);
        cd a1 = (P[k] + conj(P[r])) / 2.0, a0 = (P[k] - conj(P[r])) / cd(0, 2);
        cd b1 = (Q[k] + conj(Q[r])) / 2.0, b0 = (Q[k] - conj(Q[r])) / cd(0, 2);
        X[k] = a1 * b1 + cd(0, 1) * (a1 * b0); Y[k] = a0 * b1 + cd(0, 1) * (a0 * b0);
    }
    fft(X, true); fft(Y, true);
    vector<ll> res(need);
    for (int i = 0; i < need; i++) {
        ll hh = llround(X[i].real()) % M, hl = llround(X[i].imag()) % M, lh = llround(Y[i].real()) % M, ll_ = llround(Y[i].imag()) % M;
        res[i] = ((hh * cut % M * cut + (hl + lh) % M * cut + ll_) % M + M) % M;
    }
    return res;
}

Proof: Cooley–Tukey, the inverse, the convolution theorem and the packing tricks

Butterfly. Split by parity: $X_k={\sum }_j{x}_{2j}{\omega }^{2jk}+{\omega }^k{\sum }_j{x}_{2j+1}{\omega }^{2jk}=E_k+{\omega }^k{O}_k$, where $E,O$ are DFTs of length $n/2$ because ${\omega }^2$ is a primitive $(n/2)$-th root of unity; since ${\omega }^{n/2}=-1$, $X_{k+n/2}=E_k-{\omega }^k{O}_k$. Hence $T(n)=2T(n/2)+O(n)=O(n\log n)$. The recursion reads the input in bit-reversed order, so permuting first and running the levels bottom-up (rt[k + j] $=e^{i\pi j/k}$, a primitive $2k$-th root to the power $j$) is the same computation. Inverse: ${\sum }_k{\omega }^{(j-l)k}=n[j\equiv l]$ (geometric series of a root of unity $\ne 1$), so applying the DFT with ${\omega }^{-1}$ and dividing by $n$ inverts it; the DFT with ${\omega }^{-1}$ is the DFT with $\omega$ followed by the index reversal $j\to -j\mathrm{mod}n$, which is what the code does. Convolution theorem: ${\sum }_k{A}_k{B}_k{\omega }^{-lk}={\sum }_{i,j}{a}_i{b}_j{\sum }_k{\omega }^{(i+j-l)k}=n{\sum }_{i+j\equiv l}{a}_i{b}_j$, a cyclic convolution; padding to $n\ge |a|+|b|-1$ makes it the ordinary one (no index wraps).

Packing two real inputs. For real $u$, $\overline{U_{-k}}=U_k$. With $f=a+ib$, $F_k=A_k+i{B}_k$ and $\overline{F_{-k}}=A_k-i{B}_k$, so $A_k=\frac{F_k+\overline{F_{-k}}}{2}$, $B_k=\frac{F_k-\overline{F_{-k}}}{2i}$ and $A_k{B}_k=\frac{(x+y)(x-y)}{4i}$ with $x=F_k$, $y=\overline{F_{-k}}$. Splitting (convMod): write $a=a_1\cdot 2^{15}+a_0$, $b=b_1\cdot 2^{15}+b_0$; then $ab=a_1{b}_1{2}^{30}+(a_1{b}_0+a_0{b}_1)2^{15}+a_0{b}_0$. Each of the four convolutions has terms below $2^{30}$, so for $n\le 2^{20}$ the exact values stay below $2^{50}$, inside double precision with long-double roots. Linearity of the inverse transform lets one complex inverse carry two real results (real and imaginary parts), giving 2 forward + 2 inverse FFTs. Precision: the rounding error grows roughly like $\epsilon \log n\cdot \max |a|\max |b|\cdot n$; with $n=10^6$ and coefficients $10^9$ a plain FFT loses the low bits, which is why convExact has a bound and convMod splits. All three were checked against the $O(n^2)$ product (random lengths up to $2\cdot 10^5$, modulus $10^9+7$).

Convolution: recognising and choosing

$c_k={\sum }_{i+j=k}{a}_i{b}_j$ is the single most important operation of this chapter: recognising one turns a double loop into $O(n\log n)$. The answer at index $k$ sums over all pairs of indices combining to $k$.

Problem	Combination
polynomial multiplication; ways to get sum $k$ from two sets; subsequence-pattern counts	$i+j=k$
pairs at distance $k$	$i-j=k$ (reverse one array)
string matching with wildcards	correlation (Other Matching)
big integers	$i+j=k$ with carries
XOR / AND / OR pairs	$i\oplus j$, $i\&j$, $i\mid j$ (Bitwise Convolutions)
gcd / product pairs	$\gcd (i,j)=k$, $ij=k$: Dirichlet (Arithmetic Functions)

Combining rule	Transform	Cost
$i+j=k$	FFT / NTT	$O(n\log n)$
$i+j\equiv k(\mathrm{mod}n)$ (cyclic)	DFT of length exactly $n$ (Bluestein) or linear product then fold	$O(n\log n)$
$i\oplus j=k$	Walsh–Hadamard	$O(n\log n)$
$i\mid j=k$, $i\&j=k$	subset / superset zeta and Möbius	$O(n\log n)$
$i\cup j=k$ with $i\cap j=\varnothing$	subset convolution (ranked zeta)	$O(2^n{n}^2)$
$ij=k$	Dirichlet convolution	$O(n\log n)$
$\max (i,j)=k$	prefix sums	$O(n)$
digitwise sum mod $b$, digitwise max	base-$b$ DFT / zeta (Convolution Variants)	$O(n{\log }_{b}n\cdot b^2)$
$(\min ,+)$ / $(\max ,+)$, both convex	merge slopes (Minkowski)	$O(n+m)$
$(\min ,+)$, one convex	SMAWK / monotone minima	$O((n+m)\log )$
$(\min ,+)$, general	no truly subquadratic algorithm known	$O(nm)$

All transform-based ones share the shape transform(a); transform(b); c[i] = a[i] * b[i]; inverse(c), only the transform differs. Without a transform: sparse inputs with $s$ non-zeros cost $O(sn)$ directly; below $n\approx 32$–64 the quadratic loop wins; convex max-plus merges difference sequences (Minkowski Sum); prefix-shaped kernels are prefix sums. Multi-dimensional: transform along each axis, or flatten with a stride so rows do not wrap. Precision and modulus: small integer coefficients $\to$ NTT; arbitrary prime modulus $\to$ three NTT primes + CRT, or convMod; real coefficients $\to$ complex FFT; huge integers without modulus $\to$ split into 15-bit halves. The general lesson: when a problem says “count the ways to combine one thing from each of two sets”, write down the generating function; the answer is a coefficient of a product (Generating Functions).

Optimisation	Effect
precompute the roots once	removes a large constant (the code above caches them)
lazy reduction (Barrett / Montgomery, Fast Modular Multiplication)	fewer % in NTT
skip the inverse transform when one coefficient is needed	${\sum }_k{A}_k{B}_k{\omega }^{-lk}/n$ directly
reuse transforms; squaring needs one forward transform
pack two real sequences into one complex FFT	halves the forward work (convReal)

Karatsuba and big multiplication

Multiplying $n$-digit numbers is polynomial multiplication plus carry propagation. All of the ladder below is evaluate, multiply pointwise, interpolate: Karatsuba evaluates at $0,1,\infty$; Toom-3 at $0,1,-1,2,\infty$; Toom-$k$ at $2k-1$ small points; FFT at the $n$ roots of unity, where both evaluation and interpolation cost $O(n\log n)$ instead of $O(n)$ with a large constant, which is why FFT eventually wins. Complex multiplication in 3 real multiplications is the smallest instance of the trick; Strassen applies it to $2\times 2$ blocks.

Method	Time	Typical crossover (limbs)
schoolbook	$O(n^2)$	$<32$
Karatsuba	$O(n^{{\log }_{2}3})=O(n^{1.585})$	32 – 500
Toom-3 / Toom-4	$O(n^{1.465})$ / $O(n^{1.404})$ (Toom-3)	500 – 10000
FFT / NTT	$O(n\log n)$	$>10^4$
Schönhage–Strassen	$O(n\log n\log \log n)$	very large
Fürer; Harvey–van der Hoeven	$O(n\log n{2}^{O({\log }^{\ast }n)})$; $O(n\log n)$	galactic

vector<ll> karatsuba(vector<ll> a, vector<ll> b) {          // exact integer polynomial product
    if (a.empty() || b.empty()) return {};
    int n = max(sz(a), sz(b));
    if (min(sz(a), sz(b)) <= 32) { vector<ll> r(sz(a) + sz(b) - 1, 0); for (int i = 0; i < sz(a); i++) for (int j = 0; j < sz(b); j++) r[i + j] += a[i] * b[j]; return r; }
    int m = (n + 1) / 2, need = sz(a) + sz(b) - 1; a.resize(2 * m, 0); b.resize(2 * m, 0);
    vector<ll> a0(a.begin(), a.begin() + m), a1(a.begin() + m, a.end()), b0(b.begin(), b.begin() + m), b1(b.begin() + m, b.end());
    auto z0 = karatsuba(a0, b0), z2 = karatsuba(a1, b1);
    for (int i = 0; i < m; i++) a0[i] += a1[i], b0[i] += b1[i];
    auto z1 = karatsuba(a0, b0);
    for (int i = 0; i < sz(z0); i++) z1[i] -= z0[i];
    for (int i = 0; i < sz(z2); i++) z1[i] -= z2[i];
    vector<ll> r(4 * m, 0);
    for (int i = 0; i < sz(z0); i++) r[i] += z0[i];
    for (int i = 0; i < sz(z1); i++) r[i + m] += z1[i];
    for (int i = 0; i < sz(z2); i++) r[i + 2 * m] += z2[i];
    r.resize(need); return r;
}
 
string multiplyDecimal(const string& s, const string& t) {   // non-negative decimal strings, base-100 limbs + FFT
    auto limbs = [](const string& x) { vector<ll> v; for (int i = sz(x); i > 0; i -= 2) v.push_back(stoi(x.substr(max(0, i - 2), i - max(0, i - 2)))); return v; };
    vector<ll> c = convExact(limbs(s), limbs(t));
    ll carry = 0;
    for (auto& x : c) { x += carry; carry = x / 100; x %= 100; }
    while (carry) c.push_back(carry % 100), carry /= 100;
    while (sz(c) > 1 && c.back() == 0) c.pop_back();
    string r = to_string(c.back());
    for (int i = sz(c) - 2; i >= 0; i--) { string d = to_string(c[i]); r += string(2 - sz(d), '0') + d; }
    return r;
}

Base choice. With base $B$ and $n$ limbs the convolution values reach $n{B}^2$; doubles hold about 53 bits, so $n{B}^2\ll 2^{53}$: base 100 is safe far beyond $10^6$ digits, base $10^4$ is marginal at $n=10^6$; NTT (or three primes) has no such limit. Karatsuba only wins above a threshold (32–64 limbs), below which recursion overhead loses to a tight schoolbook loop; every serious bignum library (GMP, Java BigInteger, Python int) switches between schoolbook, Karatsuba, Toom and FFT at tuned cutoffs.

Big-integer operation	Method	Cost
add, subtract	schoolbook	$O(n)$
multiply	as above	$O(n\log n)$
divide, modulo	Newton iteration on the reciprocal $x\leftarrow x(2-bx)$ (correct bits double)	$O(M(n))$
square root	Newton	$O(M(n))$
base conversion	divide and conquer on powers of the base	$O(M(n)\log n)$
gcd	binary gcd / half-gcd	$O(n^2)$ / $O(M(n)\log n)$
modular exponentiation	square and multiply with Montgomery	$O(M(n)\log e)$

In contests you rarely need big integers: answers modulo a prime, __int128 up to $1.7\cdot 10^{38}$ for comparisons like $ab$ vs $cd$, fast exponentiation, Legendre’s formula for factorial valuations; when exact huge outputs are required (exact Fibonacci or factorials, $10^4$-digit answers), a base-$10^9$ schoolbook implementation (Big Integers and Fractions) suffices up to a few thousand digits, and FFT multiplication beyond.

Proof: Karatsuba and carry propagation

With $x=x_1{B}^m+x_0$, $y=y_1{B}^m+y_0$: $xy=x_1{y}_1{B}^{2m}+(x_1{y}_0+x_0{y}_1)B^m+x_0{y}_0$ and $x_1{y}_0+x_0{y}_1=(x_1+x_0)(y_1+y_0)-x_1{y}_1-x_0{y}_0$, so three half-size products suffice; $T(n)=3T(n/2)+O(n)=\Theta (n^{{\log }_{2}3})$ by the master theorem. For polynomials $B$ is the variable $x$ and the identity is the same. Decimal product: the convolution of limb sequences gives the coefficients of the product in base 100 before carrying; propagating carries from the lowest limb keeps the value $\sum c_{i}100^i$ unchanged and makes every limb lie in $[0,100)$, the unique base-100 representation. Checked against __int128 products and the quadratic method.

Related notes

• NTT (exact modular convolution, three-prime CRT)
• Power Series (everything built on multiplication)
• Theoretical Algorithms (Toom-3, Good–Thomas, Rader)