Math

Modular arithmetic, number theory, combinatorics, and convolution.

Formula reference

Binomial and multinomial identities

Unless stated otherwise, all parameters are nonnegative integers and \(\binom nk=0\) when \(k\notin[0,n]\).

\[ \binom nk=\binom{n-1}{k-1}+\binom{n-1}k,\qquad \binom nk=\binom n{n-k},\qquad \sum_{k=0}^{n}\binom nk=2^n. \]

Vandermonde and its central special case:

\[ \sum_k\binom rk\binom s{n-k}=\binom{r+s}n,\qquad \sum_{k=0}^{n}\binom nk^2=\binom{2n}n. \]

Hockey-stick, weighted sums, and binomial inversion:

\[ \sum_{i=r}^{n}\binom ir=\binom{n+1}{r+1},\qquad \sum_{k=0}^{n}k\binom nk=n2^{n-1},\qquad \sum_{k=0}^{n}k(k-1)\binom nk=n(n-1)2^{n-2}. \]

\[ g_n=\sum_{k=0}^{n}\binom nk f_k \iff f_n=\sum_{k=0}^{n}(-1)^{n-k}\binom nk g_k. \]

\[ \binom{n}{n_1,\ldots,n_m}=\frac{n!}{\prod_i n_i!} \quad\left(\sum_i n_i=n\right),\qquad (x_1+\cdots+x_m)^n =\sum_{n_1+\cdots+n_m=n}\binom{n}{n_1,\ldots,n_m}\prod_i x_i^{n_i}. \]

Stirling, Eulerian, and Bell numbers

The unsigned first-kind numbers count permutations with \(k\) cycles; the second-kind numbers count partitions into \(k\) nonempty unlabeled blocks:

\[ c(n,k)=c(n-1,k-1)+(n-1)c(n-1,k),\qquad \sum_{k=0}^{n}c(n,k)x^k=x(x+1)\cdots(x+n-1). \]

\[ S(n,k)=S(n-1,k-1)+kS(n-1,k) =\frac1{k!}\sum_{j=0}^{k}(-1)^{k-j}\binom{k}{j}j^n. \]

\[ x^n=\sum_{k=0}^{n}S(n,k)x(x-1)\cdots(x-k+1),\qquad B_n=\sum_{k=0}^{n}S(n,k),\qquad \sum_{n\ge0}B_n\frac{x^n}{n!}=e^{e^x-1}. \]

\[ \text{For }n\ge1,\quad \left\langle{n\atop k}\right\rangle =(n-k)\left\langle{n-1\atop k-1}\right\rangle +(k+1)\left\langle{n-1\atop k}\right\rangle =\sum_{j=0}^{k}(-1)^j\binom{n+1}{j}(k+1-j)^n. \]

Here \(\left\langle{n\atop k}\right\rangle\) counts permutations of \([n]\) with exactly \(k\) descents; it is zero outside \(0\le k<n\).

Common counting formulas

\[ C_n=\frac1{n+1}\binom{2n}{n},\qquad C_{n+1}=\sum_{i=0}^{n}C_iC_{n-i},\qquad D_n=(n-1)(D_{n-1}+D_{n-2}) =n!\sum_{k=0}^{n}\frac{(-1)^k}{k!}. \]

\[ \text{nonnegative stars and bars}=\binom{n+k-1}{k-1},\qquad \text{positive stars and bars}=\binom{n-1}{k-1}. \]

\[ p(0)=1,\quad p(n)=0\ (n<0),\qquad \sum_{n\ge0}p(n)x^n=\prod_{k\ge1}\frac1{1-x^k},\qquad p(n)=\sum_{k\ne0}(-1)^{k+1}p\!\left(n-\frac{k(3k-1)}2\right). \]

Burnside, necklaces modulo rotation with \(q\) colors, and labeled trees:

\[ |X/G|=\frac1{|G|}\sum_{g\in G}|\operatorname{Fix}(g)|,\qquad N_q(n)=\frac1n\sum_{d\mid n}\varphi(d)q^{n/d}. \]

\[ \#\text{ trees}=n^{n-2},\qquad \#\text{ trees with degrees }d_i =\frac{(n-2)!}{\prod_i(d_i-1)!} \quad\left(d_i\ge1,\ \sum_i d_i=2n-2\right). \]

For the Laplacian \(L\) of an undirected graph, let \(L^{(r)}\) denote the matrix obtained by deleting row and column \(r\). Then the Matrix-Tree theorem gives

\[ \tau(G)=\det L^{(r)}. \]

For a strongly connected Eulerian directed multigraph, fix an edge \(e\) whose tail is \(w\). Let \(t_w(G)\) be the number of spanning in-arborescences rooted at \(w\), and let \(\operatorname{EC}_e(G)\) count edge-distinct Euler circuits whose first edge is \(e\). The BEST theorem states

\[ \operatorname{EC}_e(G)=t_w(G) \prod_{v\in V}(\operatorname{outdeg}(v)-1)!. \]

Finite sums and inclusion-exclusion

\[ \sum_{i=0}^{n}x^i=\frac{x^{n+1}-1}{x-1},\qquad \sum_{i=0}^{n}ix^i =\frac{x-(n+1)x^{n+1}+nx^{n+2}}{(1-x)^2}. \]

\[ \sum_{i=1}^{n}i=\frac{n(n+1)}2,\qquad \sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}6,\qquad \sum_{i=1}^{n}i^3=\left(\frac{n(n+1)}2\right)^2. \]

\[ \left|\bigcup_{i=1}^{n}A_i\right| =\sum_{\varnothing\ne I\subseteq[n]}(-1)^{|I|+1} \left|\bigcap_{i\in I}A_i\right|. \]

Number theory

\[ \gcd(a,b)\operatorname{lcm}(a,b)=|ab|,\qquad \varphi(n)=n\prod_{p\mid n}\left(1-\frac1p\right),\qquad \sum_{d\mid n}\varphi(d)=n. \]

\[ \gcd(a,m)=1\Longrightarrow a^{\varphi(m)}\equiv1\pmod m,\qquad a^{p-1}\equiv1\pmod p\quad(p\text{ prime},\ p\nmid a). \]

\[ g(n)=\sum_{d\mid n}f(d) \iff f(n)=\sum_{d\mid n}\mu(d)g(n/d),\qquad \sum_{d\mid n}\mu(d)=[n=1]. \]

For coprime \(m,n\), if \(x\equiv a\pmod m\) and \(x\equiv b\pmod n\):

\[ x\equiv a+m\bigl((b-a)m^{-1}\bmod n\bigr)\pmod{mn}. \]

Generating functions and interpolation

\[ [x^n](A(x)B(x))=\sum_{k=0}^{n}a_kb_{n-k},\qquad a_n=\sum_{i=1}^{d}c_i a_{n-i} \Longrightarrow A(x)=\frac{P(x)}{1-c_1x-\cdots-c_dx^d}. \]

\[ P(x)=\sum_{i=0}^{n}y_i \prod_{\substack{0\le j\le n\\j\ne i}} \frac{x-x_j}{x_i-x_j},\qquad P(x_i)=y_i. \]

Newton identities, with power sums \(p_i=\sum_jx_j^i\) and elementary symmetric polynomials \(e_i\):

\[ ke_k=\sum_{i=1}^{k}(-1)^{i-1}e_{k-i}p_i. \]

Probability

\[ \mathbb E\!\left[\sum_iX_i\right]=\sum_i\mathbb E[X_i],\qquad \operatorname{Var}(X)=\mathbb E[X^2]-\mathbb E[X]^2,\qquad \mathbb E[\mathbf1_A]=\Pr(A). \]

\[ X\sim\operatorname{Bin}(n,p):\quad \Pr(X=k)=\binom nkp^k(1-p)^{n-k},\quad \mathbb E[X]=np,\quad\operatorname{Var}(X)=np(1-p). \]

\[ \Pr(X=k)_{\rm hypergeom} =\frac{\binom Kk\binom{N-K}{n-k}}{\binom Nn},\qquad \Pr(A\mid B)=\frac{\Pr(B\mid A)\Pr(A)}{\Pr(B)}. \]

Headers

Header	Summary
adjugate_matrix.hpp	Compute the adjugate of a square matrix over a field. Gauss--Jordan elimination tracks the row-operation matrix E with E A = R. For full rank, E = A^{-1}. For rank n-1, adj(R) is an outer product of the unique right-null vector and the final zero-row selector; multiplying by E recovers adj(A). Rank at most n-2 makes every (n-1)-minor zero.
batch_inverse.hpp	Invert a list of nonzero field elements with one division and O(n) multiplications.
big_integer_addition.hpp	Add two arbitrarily long signed integers represented in base 2..36. The result is canonical (no leading zeroes and no negative zero).
big_integer_division.hpp	Divide nonnegative decimal integers with base-1e9 limbs. Small divisors use normalized long division; large divisors obtain a fixed-point reciprocal by Newton doubling, using three-prime NTT for every large product, and correct the quotient by at most a constant number of steps.
big_integer_multiplication.hpp	Multiply two arbitrarily long signed integers in base 2..36. Digits are packed into limbs whose square fits exact three-prime CRT convolution; the product coefficients are then normalized by a linear carry pass.
binomial.hpp	Precomputed binomial coefficient table over a modular type T.
binomial_mod.hpp	Preprocess one fixed composite modulus for many binomial queries. Each prime-power factor stores prefix products with multiples of its prime removed. Splitting n! into complete residue blocks and recursively stripping one prime from every multiple yields its unit part and p-adic valuation in logarithmic time. The prime-power answers are joined by CRT.
bivariate_formal_power_series.hpp	Multiply two dense bivariate polynomials with rectangular truncation. Kronecker substitution stores x^i y^j at the univariate index iW+j, choosing W as the complete product width. Then no y-degree sum carries into the next x digit, so one ordinary convolution contains the exact bivariate product; unpacking only the requested rectangle performs the truncation.
characteristic_polynomial.hpp	Return coefficients c_i of det(xI-A)=sum c_i x^i over a field in O(n^3), using a similarity reduction to upper Hessenberg form.
combinatorial_sequences.hpp	Return B_0 through B_n. Their exponential generating function is exp(exp(x)-1), so one FPS exponential followed by factorial scaling yields all Bell numbers simultaneously.
convolution.hpp	Bitwise AND, OR, XOR, and subset convolutions for power-of-two indexed arrays.
convolution_f2_64.hpp	Element of GF(2^64) represented modulo x^64+x^4+x^3+x+1. Addition is XOR; multiplication uses a carry-less product and folds its high half twice through x^64=x^4+x^3+x+1.
convolution_large.hpp	Convolve modint sequences whose result exceeds one supported NTT. Split each input into half-transform-size blocks. Every block pair then fits in one transform; pairs with the same block-index sum share a frequency-domain accumulator, and their inverse transform is overlap-added at the corresponding coefficient offset. The optional transform size is useful for testing the block boundary and must be a supported power of two.
convolution_mod.hpp	Convolve nonnegative residues modulo an arbitrary modulus using three NTT primes. Requires result length <= 2^24 and each integer coefficient before reduction to be smaller than the product of the primes.
convolution_mod_2_64.hpp	Convolve unsigned 64-bit integers modulo 2^64. The convolution is computed modulo five NTT primes whose product exceeds every exact coefficient for input lengths at most 2^19. Garner's mixed-radix digits therefore recover the exact coefficient; accumulating those digits in a uint64_t keeps precisely its residue modulo 2^64.
decimal_power_mod.hpp	Compute base^exponent modulo an arbitrary positive 64-bit modulus from a decimal exponent of any length, without coprimality assumptions.
determinant_mod.hpp	Determinant of a square integer matrix modulo any positive modulus in O(n^3 log modulus), without requiring modular inverses.
dirichlet_prefix.hpp	Return Q_n={floor(n/i)} in increasing order. Consecutive equal quotients are represented once.
discrete_log.hpp	Find the smallest x >= 0 satisfying base^x = target (mod modulus), or nullopt if no solution exists.
distinct_subsequences.hpp	Count distinct nonempty subsequences over an arbitrary modular field.
divisor_summatory.hpp	Return sum_{i=1}^n tau(i) in O(sqrt(n)) quotient blocks.
divisor_table.hpp	Linear-sieve table of divisor count tau(n) and divisor sum sigma(n) for every integer through a limit.
dujiao_sieve.hpp	Summatory Mobius and Euler-phi functions with a linear-sieve prefix and quotient-block memoization for large arguments.
exponential_polynomial_sum.hpp	Return sum_{i>=0} r^i i^d over a finite field for r != 1, with 0^0=1. A degree-d polynomial is determined by its values at 0..d. Writing it in the consecutive-point Lagrange basis and summing each basis polynomial against the geometric series reduces the answer to factorial weights and prefix sums of r^i i^d.
extended_gcd.hpp	Extended Euclid, linear Diophantine equations, and Chinese remainder routines.
f2_matrix.hpp	Return the rank of a packed binary matrix. Gaussian elimination uses XOR as row addition; when the matrix is tall, transposition first minimizes the number of packed rows participating in elimination.
factorize.hpp	Deterministic 64-bit primality testing and integer factorization with Pollard's rho.
finite_field_polynomial_factorization.hpp	Cantor-Zassenhaus factorization of polynomials over arbitrary prime fields.
fixed_discrete_log.hpp	Preprocess discrete logarithms to one primitive root modulo a prime. Farey approximation writes every nonzero x as xq = +/-t (mod p) with q <= p^(1/3) and |t| <= p^(2/3). Logs of the small integers are built from a batched baby-step/giant-step pass for primes, multiplicativity for small composites, and p = floor(p/i)i + (p mod i) for the remaining interval. A query then needs only the precomputed logs of q and |t|.
floor_monoid_product.hpp	Universal Euclidean algorithm for noncommutative monoid products along a floor path.
floor_sum.hpp	Signed floor sums over an arbitrary half-open integer interval.
formal_power_series.hpp	Formal power series logarithm, exponential, power, and square root operations.
formal_power_series_composition.hpp	Compose two formal power series modulo x^n over 998244353. After shifting away g(0), transposition turns composition into simultaneous coefficient extraction from powers of g. A bivariate Bostan-Mori recursion halves the x degree at every level; reversing that recursion applies the transposed two-dimensional convolutions and reconstructs all coefficients.
fraction.hpp	Normalized signed 64-bit rational number; arithmetic results must fit int64, while comparison and intermediate products use signed 128-bit.
gaussian_integer.hpp	Euclidean division in Z[i], returning quotient and remainder with norm(remainder) < norm(divisor).
hafnian.hpp	Compute the hafnian of an even-order symmetric matrix over a commutative ring. A paired-vertex elimination keeps every remaining edge weight as a polynomial in a marker z. Eliminating vertices u,v adds the two possible paths i-u-v-j to each surviving edge; an include/exclude recurrence cancels structures that reuse a vertex. After n/2 eliminations, the coefficient of z^(n/2) contains exactly the weighted perfect matchings.
hex_big_integer_division.hpp	Divide nonnegative hexadecimal integers stored in 32-bit limbs. Karatsuba supplies subquadratic multiplication, while Burnikel-Ziegler splits a normalized 2n-by-n division into balanced 3n/2-by-n and n-by-n/2 subproblems; short blocks fall back to Knuth-style long division.
increasing_sequence_bounds.hpp	Convolution-accelerated counting of bounded nondecreasing sequences.
integer_kth_root.hpp	Return floor(value^(1/exponent)) for an unsigned 64-bit integer. Binary search uses an exponent-dependent upper bound. The comparison multiplies only while the next factor is at most value/current, so it is exact and never relies on floating-point rounding or overflowing products.
lagrange_interpolation.hpp	Evaluate the degree < values.size() polynomial known at consecutive points 0,1,... in O(n) over a field.
large_factorial.hpp	Compute several factorials modulo a fixed prime without a linear table. Split 1..N into blocks of length B about sqrt(N). The product inside one block is the degree-B polynomial P(x)=prod_{i=1}^B(x+i); a product tree constructs P and multipoint evaluation obtains P(0),P(B),P(2B),.... Prefix products answer every full block, followed by at most B direct factors.
linear_algebra.hpp	Matrix rank, determinant, linear-system solving, and nullspace construction over a field.
linear_recurrence.hpp	Recover the shortest linear recurrence of a sequence over a field.
matrix.hpp	Dense matrix multiplication in O(nmk).
min_of_mod_linear.hpp	Return an index attaining the minimum of (multiplier * x + addition) mod modulus for 0 <= x < count, together with that minimum. Continued-fraction neighbors describe every index at which the prefix minimum decreases; those indices form arithmetic runs, so only one run per Euclidean-algorithm step is inspected.
minmax_inclusion_exclusion.hpp	Coefficients for extended min-max inclusion-exclusion: (-1)^(s-k) C(s-1,k-1) for subset size s.
mod_kth_root.hpp	Return one x satisfying x^k = value modulo a prime, or nullopt; avoids a full discrete logarithm by extracting the prime-power factors of gcd(k, p - 1) independently.
mod_power_of_two_multiplicative_convolution.hpp	Convolve functions on the multiplicative monoid modulo 2^n. Every nonzero residue is 2^v times an odd unit. For modulus at least 8, the odd units form C_2 x C_(2^(n-v-2)), generated by -1 and 5. We group coefficients by the valuation v, Fourier-transform both unit-group coordinates, and add products from layers v_1,v_2 to layer min(n,v_1+v_2). The last layer represents zero.
mod_sqrt.hpp	Compute the smaller square root modulo a prime, or nullopt if no square root exists.
modular_multiplicative_convolution.hpp	Multiplicative convolution of arrays indexed modulo a prime p. The nonzero residues form a cyclic group generated by g, so mapping g^e to exponent e turns multiplication into addition modulo p-1. One ordinary NTT convolution followed by folding at p-1 computes every nonzero coefficient; pairs containing a zero index contribute only to coefficient zero.
montmort.hpp	Return derangement counts D_0 through D_n modulo modulus. Separating the position occupied by the image of the first element gives D_n = (n-1)(D_{n-1}+D_{n-2}), with D_0=1 and D_1=0.
multiplicative_function_sum.hpp	Sum a multiplicative function satisfying f(p^e)=ae+bp for every prime power. Quotient blocks first store both pi(x) and the sum of primes up to x. A recursion then fixes the least prime factor and its exponent; the remaining one-prime tail is aggregated from those two tables instead of visited individually. Each integer is represented by its unique prime factorization, so the accumulated contributions are disjoint and complete.
multivariate_convolution.hpp	Multiply dense multivariate polynomials truncated by x_i^n_i. Mixed-radix indices turn a monomial product into ordinary index addition, except that an overflowing coordinate creates a carry into the next digit. Color index j by sum floor(j/(n_1...n_i)) modulo k. The color difference between two factors and their sum equals the number of carries modulo k; since at most k-1 carries occur, retaining only the expected color removes every wrapped term. NTT performs all ordinary convolutions simultaneously.
multivariate_cyclic_convolution.hpp	Multiply dense multivariate polynomials modulo (1-x_1^n_1, ..., 1-x_k^n_k) over a prime field. A primitive root supplies an n_i-th root of unity for each coordinate. Applying a finite-field DFT along every coordinate changes circular convolution into pointwise multiplication; inverse coordinate transforms recover the coefficients. Each possibly non-power-of-two DFT is evaluated by a chirp convolution.
newton_basis.hpp	Convert monomial coefficients of f into coefficients b_i of f(x)=sum_i b_i product_{j<i}(x-p_j). Pad to a power of two and build the product tree of (x-p_i). Reversed formal division extracts the coefficient vector at the root. When descending, multiplying by the sibling product and taking the middle half transfers the contribution to the left child; the upper half transfers it to the right child. Leaves are the Newton coefficients in order.
nim_product.hpp	Compute the 64-bit nim product. Nimbers with twice as many bits form a quadratic extension of the lower field: for alpha=2^m, alpha^2=alpha xor 2^(m-1). Recursing on this identity generates all products of binary basis elements. Bilinearity then groups eight basis elements at a time into a 32 MiB byte-pair table, so one query needs 8*8 lookups and XORs.
pfaffian.hpp	Compute the Pfaffian of an even-order alternating matrix over a field. Each step moves a nonzero entry into the next 2-by-2 pivot block; the skew-symmetric Schur complement removes that pair, and the product of pivot entries is the Pfaffian. Swapping a paired index changes its sign.
polynomial.hpp	Dense polynomial arithmetic, inversion, division, and remainder operations.
polynomial_gcd.hpp	Fast half-gcd, modular inverse, and prime-field root finding for polynomials.
polynomial_multipoint.hpp	Product-tree multipoint evaluation and interpolation over a field.
polynomial_special_points.hpp	Recover f(c),...,f(c+m-1) from f(0),...,f(n-1). Lagrange weights turn every non-sampled value into one convolution with 1/(c+k-i); a sliding product supplies prod_j(c+k-j). Positions that coincide with an original sample are copied directly.
prime.hpp	A linear sieve providing primes, smallest prime factors, Euler phi, and Mobius values.
prime_binomial_table.hpp	Binomial coefficients modulo a runtime prime p for arguments below p. Store factorials and inverse factorials through N; Fermat inversion of N! followed by a backward sweep obtains every inverse using one power.
prime_summatory.hpp	Compute pi(n) and the sum of primes at most n with a quotient-block combinatorial sieve using O(sqrt(n)) memory.
primitive_root.hpp	Return the smallest primitive root modulo a prime.
q_binomial_prime.hpp	Answer Gaussian binomial coefficients at a fixed q modulo a prime. If d is the multiplicative order of q, q-Lucas gives [n choose k]_q = C(floor(n/d),floor(k/d)) [n mod d choose k mod d]_q. Ordinary factorials handle the first factor, while products of 1-q^i for i<d handle the second. q=0 and q=1 are treated by their direct limits.
range_product_changes.hpp	Enumerate the distinct products of ranges sharing a left endpoint.
rational_approximation.hpp	Find the closest reduced fractions on both sides of x/y whose numerators and denominators are at most limit. Alternating jumps between the two Stern--Brocot boundaries consumes one continued-fraction run at a time; each run length is maximized by integer binary search. Missing lower and upper answers are represented by 0/1 and 1/0.
rational_reconstruction.hpp	Reconstruct a reduced fraction numerator/denominator congruent to a residue modulo modulus within supplied bounds. Uniqueness is guaranteed when 2numerator_bounddenominator_bound < modulus.
segmented_sieve.hpp	Call callback(index, prime) for every prime not exceeding n and return pi(n). The callback is invoked in increasing prime order.
set_power_series.hpp	Compose an exponential generating function with a set power series. A coefficient indexed by a mask represents the square-free monomial whose variables are that mask. Splitting masks by their largest variable and by popcount turns every multiplication into ranked subset convolution. The subset zeta transform performs all disjoint decompositions at once, while the rank selects exactly the square-free terms. egf[k] is the multiplier of series^k / k!, and the constant term of series must be zero.
sparse_formal_power_series.hpp	Invert a sparse formal power series. From f g = 1, each new coefficient of g is a dot product against the nonconstant terms of f.
sparse_matrix_determinant.hpp	Compute a sparse matrix determinant over a finite field by scalar Wiedemann projection. Multiplying the columns by independent nonzero random values makes the characteristic polynomial equal the minimal polynomial with high probability. The sequence u^T(AD)^k v exposes that polynomial to Berlekamp--Massey; its constant coefficient is det(AD) up to the usual characteristic-polynomial sign. Dividing by the diagonal product recovers det(A). A projection of too small a degree is retried.
squarefree_count.hpp	Count square-free positive integers at most n. Mobius inversion handles divisors up to D = sqrt(n / n^(1/5)) directly. The remaining terms share only O(n^(1/5)) distinct square-root quotients; their Mertens values are recovered in descending order with the quotient-block recurrence, reducing both time and memory from O(sqrt(n)) to O(n^(2/5)).
stern_brocot.hpp	Run-length encoded L/R path from 1/1 to a positive reduced fraction.
stirling_small_prime.hpp	Query signed first-kind and second-kind Stirling numbers modulo a small fixed prime for indices as large as 10^18. Over F_p, the identities (x)_p=x^p-x and x^p=x split a large falling factorial or power into p-sized blocks. Repeating the split leaves one small Stirling number with indices below p and one binomial coefficient. Lucas' theorem evaluates the latter digit by digit. The constructor can omit either Stirling table when only one kind is needed.
subarray_gcd.hpp	For every right endpoint, group all nonempty suffixes by gcd in O(n log V) total groups.
subset_sum_count.hpp	Count subsets by their sum through the generating function F(x)=prod_i(1+x^s_i). Its logarithm is sum_s count(s) sum_{k>=1} (-1)^(k+1)x^(ks)/k, which is accumulated by a divisor loop. One formal-series exponential then recovers all coefficients of F up to the requested limit.
sum_two_squares.hpp	Decide and construct n = a^2 + b^2 for an unsigned 64-bit integer; also count ordered signed representations using the two-square theorem.
tetration_mod.hpp	Compute a power tower of height copies of base modulo modulus, with 0^0 = 1. Recursing through Euler-totient moduli keeps the exponent small; an extra "large" bit distinguishes a reduced exponent from a genuinely small one when the base is not coprime to the modulus.
xor_basis.hpp	Reduced linear basis over GF(2) for an unsigned integer type.