xii. Counting DP Techniques

Techniques for counting objects (mod a prime) where a direct DP would need too much state. The recurring ideas: insert elements in a clever order so the state stays small (relative ranks, connected blocks), compress huge numeric states to their $O(\sqrt{K})$ distinct floor values, factor the lattice of divisors or subsets into chains, and double the range when $n$ is astronomically large. Generating-function methods are in Generating Functions and Counting with Power Series.

Permutations by relative order

Counting permutations with constraints between neighbouring positions (an up/down signature, $|p_i-p_{i+1}|$ conditions): process positions left to right and remember only the rank of the last element among those placed so far.

// # permutations of 1..n with p_i < p_{i+1} where s[i] == '<' and > where '>'  (|s| = n - 1)
vector<mint> signature(const string& s) {
    int n = sz(s) + 1; vector<mint> dp{1};                 // dp[j]: last element has rank j among the first i
    for (int i = 1; i < n; i++) {
        vector<mint> pre(i + 1), nd(i + 1);
        for (int j = 0; j < i; j++) pre[j + 1] = pre[j] + dp[j];
        for (int j = 0; j <= i; j++) nd[j] = s[i - 1] == '<' ? pre[j] : pre[i] - pre[j];
        dp = nd;
    }
    return dp;                                             // sum of dp = answer
}

Proof: ranks are a complete state

Map the first $i$ values to their relative order (a permutation of $\mathrm{0..}i-1$). Appending a new value whose rank among the first $i+1$ is $j$ shifts all old ranks $\ge j$ up by one and leaves the relative order of the old elements unchanged; every pair (old pattern, $j$) arises from exactly one pattern of $i+1$ values. The new last element is larger than the previous last iff the previous last had old rank $<j$. So the count of patterns ending in rank $j$ depends only on the counts for the previous ranks, which is the recurrence; prefix sums make each layer $O(i)$.

Connected-component (insertion) DP

For constraints involving values of neighbours (zigzag permutations, costs $|p_i-p_{i+1}|$, “buildings visible”), insert values in increasing order into a set of $j$ blocks: contiguous pieces of the final permutation whose relative order is fixed and whose gaps will be filled by larger values later. The new value can

• start a new block (in one of the $j+1$ gaps),
• attach to an end of an existing block ($2j$ ends), or
• merge two neighbouring blocks ($j-1$ ways).

Each move fixes the value’s final neighbours relative to larger values: a merging value is larger than both neighbours, a new block’s value is smaller than both. Fixed endpoints restrict which gaps are available.

mint kangaroo(int n, int s, int t) {                       // zigzag permutations with p_1 = s, p_n = t (CEOI 2016)
    if (s == t) return n == 1;
    vector<vector<mint>> dp(n + 1, vector<mint>(n + 2)); dp[0][0] = 1;
    for (int i = 1; i <= n; i++) for (int j = 1; j <= i; j++) {
        if (i == s || i == t) dp[i][j] = dp[i - 1][j - 1] + dp[i - 1][j];   // new border block or extend the border block
        else {
            int gaps = j - (i > s) - (i > t);                 // gaps not blocked by an already placed end
            dp[i][j] = dp[i - 1][j - 1] * max(gaps, 0) + dp[i - 1][j + 1] * j;   // new block (local min) or merge (local max)
        }
    }
    return dp[n][1];
}

Costs like $\sum |p_i-p_{i+1}|$: when value $v$ is inserted, each of its open sides (sides that will later touch a larger value) contributes $-v$ and each closed side $+v$; track the running sum as an extra dimension (often bounded by “sum of open ends × remaining gap”, which keeps it small).

Proof: the block moves enumerate each permutation once

Fix a final permutation and read it while inserting values $1,2,\dots$. After inserting $\mathrm{1..}i$, the positions holding these values form maximal runs (blocks) separated by positions holding larger values; their left-to-right order is determined. Value $i+1$ either has no neighbour among $\mathrm{1..}i$ (new block, in a determined gap), exactly one (attaches to a determined end), or two (merges two determined blocks). So each permutation gives exactly one sequence of moves, and each move sequence ending in one block gives exactly one permutation. In the zigzag problem, an inner value is either a local max (both final neighbours smaller, hence already placed: merge) or a local min (both larger: new block); the ends $s$, $t$ are forced to the borders, which removes the outer gaps once they are placed.

Huge numeric states: floor values

A state like “remaining budget $\lfloor K/p\rfloor$” takes only $O(\sqrt{K})$ distinct values, because $\lfloor \lfloor K/a\rfloor /b\rfloor =\lfloor K/(ab)\rfloor$ and $\lfloor K/d\rfloor$ has at most $2\sqrt{K}$ distinct values. Index small values directly and large ones by $K/q$.

mint countProductAtMost(const vector<ll>& a, ll K) {       // # subsets (incl. empty) with product <= K
    if (K < 1) return 0;
    ll S = sqrtl(K); while ((S + 1) * (S + 1) <= K) S++; while (S * S > K) S--;
    int M = 2 * S + 2; auto id = [&](ll q) { return q <= S ? int(q) : int(S + K / q); };
    vector<ll> val(M); for (ll q = 0; q <= S; q++) val[q] = q; for (ll d = 1; d <= S; d++) val[S + d] = K / d;
    vector<mint> cnt(M); cnt[id(K)] = 1;                      // state = floor(K / product so far)
    for (ll x : a) { auto nx = cnt; for (int i = 0; i < M; i++) if (val[i] / x >= 1) nx[id(val[i] / x)] += cnt[i]; cnt = nx; }
    mint t = 0; for (auto& c : cnt) t += c; return t;
}

Proof: nested floors and the $2\sqrt{K}$ bound

Write $K=q\cdot ab+r$ with $0\le r<ab$; then $K/a=qb+r/a$ with $r/a<b$, so $\lfloor K/a\rfloor =qb+\lfloor r/a\rfloor$ with $\lfloor r/a\rfloor <b$, and dividing by $b$ gives $q$. Hence “product $\cdot$ rest $\le K$” depends on the product only through $\lfloor K/\text{product}\rfloor$. For $d\le \sqrt{K}$ there are at most $\sqrt{K}$ values $\lfloor K/d\rfloor$; for $d>\sqrt{K}$ the value is $<\sqrt{K}$. So at most $2\sqrt{K}$ distinct values.

Bounded XOR tuples

Count $(x_1,\dots ,x_n)$ with $0\le x_i\le a_i$ and $x_1\oplus \cdots \oplus x_n=X$ in $O(n\log \max a)$. At the highest bit $k$ set in some $a_i$: if some such variable puts 0 at bit $k$, its lower $k$ bits become completely free, so it can absorb any XOR of the lower bits (take the first such variable; it contributes factor 1, others contribute their counts). Otherwise all of them keep bit $k$; strip it and recurse.

ll countXor(vector<ll> a, ll X, ll mod) {
    ll res = 0;
    for (int k = 62; k >= 0; k--) {
        bool any = false; for (ll x : a) any |= x >> k & 1;
        if (!any) { if (X >> k & 1) return res; continue; }
        ll dp[2][2] = {{1, 0}, {0, 0}}, low = (1LL << k) % mod;   // [parity at bit k][free variable chosen]
        for (ll x : a) {
            ll nd[2][2] = {};
            for (int p = 0; p < 2; p++) for (int f = 0; f < 2; f++) if (ll c = dp[p][f]) {
                if (x >> k & 1) {
                    nd[p ^ 1][f] = (nd[p ^ 1][f] + c * ((x - (1LL << k) + 1) % mod)) % mod;   // keep bit k
                    nd[p][1] = (nd[p][1] + c * (f ? low : 1)) % mod;                             // drop it: lower bits free
                } else nd[p][f] = (nd[p][f] + c * ((x + 1) % mod)) % mod;
            }
            memcpy(dp, nd, sizeof dp);
        }
        res = (res + dp[X >> k & 1][1]) % mod;
        int ones = 0; for (ll x : a) ones += x >> k & 1;
        if ((ones & 1) != (X >> k & 1)) return res;
        for (ll& x : a) if (x >> k & 1) x ^= 1LL << k;
        X &= ~(1LL << k);
    }
    return (res + (X == 0)) % mod;
}

Proof: the free variable absorbs the lower bits

Split the tuples at bit $k$ by whether some variable with bit $k$ available chooses 0 there. In that case let $i$ be the first such variable: its value is $x_i=y<2^k\le a_i$ with $y$ arbitrary. Fix every other variable’s value (counted with their multiplicities: $a-2^k+1$ ways to keep the bit, $2^k$ ways to drop it, $a+1$ ways for variables without bit $k$). Then exactly one $y$ makes the XOR of the bits below $k$ equal to $X$‘s, and the bit-$k$ parity is determined by the others. So the count is the DP value with factor 1 for the first free variable. Otherwise every variable with bit $k$ available keeps it, the parity at bit $k$ is fixed, and the remaining freedom is exactly the problem on lower bits with $a_i-2^k$. Bits above the top bit must be 0 in $X$.

Divisor lattice DP

Functions on divisors of $N$ (or on $\mathrm{1..}n$ with divisibility) are functions on a product of chains (one per prime), so zeta/Möbius transforms are prefix sums/differences along each prime, exactly like SOS DP along each bit (SOS).

template<class T> void divisorZeta(vector<T>& a, const vector<int>& primes) {     // a[i] = sum_{d | i} a[d], 1..n
    int n = sz(a) - 1; for (int p : primes) for (int i = 1; (ll)i * p <= n; i++) a[i * p] += a[i];
}
template<class T> void divisorMobius(vector<T>& a, const vector<int>& primes) {
    int n = sz(a) - 1; for (int p : primes) for (int i = n / p; i >= 1; i--) a[i * p] -= a[i];
}
template<class T> void multipleZeta(vector<T>& a, const vector<int>& primes) {    // a[i] = sum_{i | j} a[j]
    int n = sz(a) - 1; for (int p : primes) for (int i = n / p; i >= 1; i--) a[i] += a[i * p];
}
template<class T> void multipleMobius(vector<T>& a, const vector<int>& primes) {
    int n = sz(a) - 1; for (int p : primes) for (int i = 1; (ll)i * p <= n; i++) a[i] -= a[i * p];
}
// gcd convolution c[k] = sum_{gcd(i,j)=k} a[i] b[j]: multipleZeta both, multiply pointwise, multipleMobius
// lcm convolution (k <= n): divisorZeta both, multiply, divisorMobius

Typical uses: number of pairs with $\gcd =d$ (multiple-Möbius of $\left(\genfrac{}{}{0ex}{}{c_d}{2}\right)$, where $c_d$ = count of multiples of $d$), $\sum \gcd (a_i,a_j)$, counting arrays with a given gcd or lcm, necklaces (Burnside). For a single $N\le 10^{12}$, do the same over the sorted divisor list of $N$ with $O(d(N)\omega (N))$ steps.

Proof: prime-by-prime prefix sums give the divisor zeta transform

Write $i=\prod p^{e_p}$. After processing primes $p_1..p_r$, $a[i]$ holds $\sum a[d]$ over divisors $d$ of $i$ that agree with $i$ on the exponents of all primes except $p_1..p_r$ (where they may be smaller). Processing $p_{r+1}$ in increasing $i$ adds $a[i/p]$ (already updated for $p_{r+1}$) to $a[i]$, which is a prefix sum along the $p_{r+1}$-exponent chain; by induction this extends the set of allowed divisors to also vary the $p_{r+1}$ exponent. After all primes, $a[i]={\sum }_{d\mid i}a[d]$. The Möbius version undoes each prefix sum by a difference in reverse order. The cost is ${\sum }_{p}n/p=O(n\log \log n)$. Gcd convolution works because ${\sum }_{k\mid g}c[g]=({\sum }_{k\mid i}{a}_i)({\sum }_{k\mid j}{b}_j)$ (both $k\mid i$ and $k\mid j$ iff $k\mid \gcd (i,j)$).

Doubling over huge $n$

If the answer for $\{\mathrm{1..}n\}$ can be updated by “append element $n+1$” and by “combine $\{\mathrm{1..}p\}$ with $\{p+\mathrm{1..2}p\}$”, compute it from the binary expansion of $n$ in $O(\log n)$ combines. Example: elementary symmetric sums $e_k(1,\dots ,n)$ for all $k\le K$ (these are unsigned Stirling numbers of the first kind $\left[\genfrac{}{}{0ex}{}{n+1}{n+1-k}\right]$), mod a prime $>K$, in $O(K^2\log n)$:

$e_i(p+1,\dots ,2p)={\sum }_{t=0}^i{e}_t(1,\dots ,p)p^{i-t}\left(\genfrac{}{}{0ex}{}{p-t}{i-t}\right).$

vector<mint> elemSym(ll n, int K) {                        // e_0..e_K of {1, ..., n}
    vector<mint> inv(K + 2, 1); for (int i = 1; i <= K + 1; i++) inv[i] = mint(1) / i;
    auto rec = [&](auto&& self, ll m) -> vector<mint> {
        vector<mint> ans(K + 1); if (m == 0) { ans[0] = 1; return ans; }
        if (m & 1) { auto l = self(self, m - 1); for (int i = 0; i <= K; i++) ans[i] = l[i] + (i ? l[i - 1] * mint(m) : 0); return ans; }
        ll p = m / 2; auto l = self(self, p); vector<mint> r(K + 1), pw(K + 1, 1);
        for (int i = 1; i <= K; i++) pw[i] = pw[i - 1] * mint(p);
        for (int j = 0; j <= K && j <= p; j++) {
            mint c = 1;                                     // C(p - j, t)
            for (int t = 0; j + t <= K && j + t <= p; t++) { r[j + t] += l[j] * pw[t] * c; c *= mint(p - j - t) * inv[t + 1]; }
        }
        for (int i = 0; i <= K; i++) for (int j = 0; i + j <= K; j++) ans[i + j] += l[i] * r[j];
        return ans;
    };
    return rec(rec, n);
}

Proof: the shift formula

$e_i(p+1,\dots ,2p)={\sum }_{|S|=i}{\prod }_{j\in S}(p+j)$ over $S\subseteq \{\mathrm{1..}p\}$. Expanding each product chooses a subset $T\subseteq S$ whose elements contribute $j$ and the rest contribute $p$: the term is $p^{i-|T|}{\prod }_{j\in T}j$. Each $T$ of size $t$ lies in $\left(\genfrac{}{}{0ex}{}{p-t}{i-t}\right)$ sets $S$ of size $i$. Summing gives the formula. Then $e(\mathrm{1..2}p)=e(\mathrm{1..}p)\ast e(p+\mathrm{1..2}p)$ as polynomials (truncated at degree $K$).

Dynamic submask counting

A multiset of $B$-bit masks with insert/erase and “how many elements are submasks of $s$” in $O(2^{B/2})$ per operation: split bits into halves; table[hi][lo] = number of elements whose high half equals hi and whose low half is a submask of lo. Insertion enumerates supermasks of its low half; a query enumerates submasks of its high half.

template<int B> struct DynSubmask {
    static constexpr int H = B / 2, L = B - H, FL = (1 << L) - 1;
    vector<int> t = vector<int>(1 << B);
    void update(int x, int d) { int hi = x >> L, lo = x & FL; for (int m = lo;; m = (m + 1) | lo) { t[hi << L | m] += d; if (m == FL) break; } }
    int countSub(int s) const {                              // # elements a with (a & s) == a
        int hi = s >> L, lo = s & FL, r = 0;
        for (int m = hi;; m = (m - 1) & hi) { r += t[m << L | lo]; if (!m) break; }
        return r;
    }
};

Proof: the split table

$a$ is a submask of $s$ iff $\text{hi}(a)\subseteq \text{hi}(s)$ and $\text{lo}(a)\subseteq \text{lo}(s)$. The query enumerates every candidate $\text{hi}(a)=m\subseteq \text{hi}(s)$ and reads $t[m][\text{lo}(s)]$, which counts elements with high half exactly $m$ and low half a submask of $\text{lo}(s)$ (insertions added 1 to every supermask of their low half). Both loops enumerate at most $2^{B/2}$ masks.

More counting patterns

• Counting distinct objects, not constructions: fix a canonical construction (e.g. the lexicographically smallest way, or “the first occurrence”) and count only those. Distinct subsequences: $d{p}_i=2d{p}_{i-1}-d{p}_{\text{prev}(i)-1}$.
• First bad event: count sequences avoiding something by subtracting, for each possible first violation, (ways to reach it cleanly) × (anything after) (Classic DP, grid paths).
• Moments over all subsets: to sum $f(S{)}^2$ over all subsets, carry (count, $\sum f$, $\sum f^2$) per state; adding an element with value $x$ maps $(c,s_1,s_2)\to (c,s_1+cx,s_2+2x{s}_1+c{x}^2)$. Higher powers via binomial expansion.
• Union sizes of all families of sets in $O(n+m{2}^m)$: an element misses the union of mask iff its membership set is a submask of the complement, so $u_{\text{mask}}=n-{\sum }_{s\subseteq \overline{\text{mask}}}cn{t}_s$, a subset zeta transform (intersections use the superset transform instead). Typical use: choose a family of sets of minimum cost covering at least $k$ elements.
• Counting labelled structures by their connected pieces: fix the component containing vertex 1 (Graph Counting).
• Permutations by inserting values in increasing order into growing contiguous blocks (a different “connected-component DP” from the graph-counting one above): $dp[i][j]$ = ways to place values $\mathrm{1..}i$ as $j$ blocks (contiguous pieces of the eventual permutation, order between blocks fixed, gaps filled later). Value $i+1$ can start a new block, extend an existing block’s end, or merge two adjacent blocks; which moves are legal, and their multiplicity (number of available ends/gaps), comes from the problem’s adjacency constraint. $O(n^2)$. Same pattern underlies permutations with a cost on $|p_i-p_{i+1}|$, “zigzag permutations with fixed endpoints” (fixed endpoints can only create/extend the outermost blocks), and problems like CF 1515E (Phoenix and Computers) and JOI’s “Skyscraper”.