ix. Counting with Power Series

Whole rows and columns of combinatorial sequences in $O(n\log n)$ (Stirling, Bell, Bernoulli, partitions, power sums), Euler / set transforms (unlabelled multisets and sets: subset-sum counting, bounded knapsack counts), $O(nk)$ arithmetic on sparse series, and sums ${\sum }_{i<n}{a}^{i}f(i)$ of exponential times polynomial for $n\le 10^{18}$. Built on Power Series, taylorShift from Composition and Taylor Shift, lagrangeConsecutive from Evaluation and Interpolation and Comb (factorials up to the size used). The $O(n^2)$ tables and the combinatorial meaning are in Stirling, Bell and Partitions.

Sequence	Generating function	Computed as
unsigned $c(n,k)$, row $n$	$x(x+1)\cdots (x+n-1)$	doubling $F_{2m}(x)=F_m(x)F_m(x+m)$
$S(n,k)$, row $n$	$S(n,k)={\sum }_j\frac{j^n}{j!}\frac{(-1{)}^{k-j}}{(k-j)!}$	one convolution
column $k$: $S(\cdot ,k)$, $c(\cdot ,k)$	EGF $\frac{(e^x-1{)}^k}{k!}$, $\frac{(-\ln (1-x){)}^k}{k!}$	one power
Bell $B_0..B_N$	EGF $e^{e^x-1}$	one exp
Bernoulli $B_0..B_N$ ($B_1=-\frac{1}{2}$)	EGF $\frac{x}{e^x-1}$	one inverse
partitions $p(\mathrm{0..}N)$	${\prod }_k\frac{1}{1-x^k}$ = inverse of the pentagonal series	one inverse
power sums ${\sum }_{i<n}{i}^k$, all $k\le K$	EGF $\frac{e^{nx}-1}{e^x-1}$	one inverse
multisets of objects with counts $a_k$ (Euler transform)	${\prod }_k(1-x^k{)}^{-a_k}=\exp {\sum }_k{a}_k{\sum }_j\frac{x^{jk}}{j}$	harmonic sum + exp
sets (distinct objects)	${\prod }_k(1+x^k{)}^{a_k}=\exp {\sum }_k{a}_k{\sum }_j\frac{(-1{)}^{j+1}{x}^{jk}}{j}$	harmonic sum + exp

Poly stirling1Row(int n) {                                  // c(n, k), k = 0..n
    if (n == 0) return {1};
    Poly half = stirling1Row(n / 2), f = polyMul(half, taylorShift(half, mint(n / 2)));
    f.resize(2 * (n / 2) + 1);
    if (n % 2) { Poly g(n + 1, 0); for (int i = 0; i < sz(f); i++) g[i] += f[i] * mint(n - 1), g[i + 1] += f[i]; f = g; }   // times (x + n - 1)
    return f;
}
Poly stirling2Row(int n, const Comb& cb) {                  // S(n, k), k = 0..n
    Poly a(n + 1), b(n + 1);
    for (int i = 0; i <= n; i++) a[i] = mint(i).pow(n) * cb.fi[i], b[i] = i % 2 ? -cb.fi[i] : cb.fi[i];
    Poly c = polyMul(a, b); c.resize(n + 1); return c;
}
Poly stirlingColumn(int k, int N, bool first, const Comb& cb) {   // S(i, k) or c(i, k), i = 0..N
    Poly res(N + 1, 0); if (k > N) return res;
    int m = N - k + 1; Poly g(m);
    for (int j = 0; j < m; j++) g[j] = first ? cb.iv[j + 1] : cb.fi[j + 1];   // (-ln(1-x))/x or (e^x - 1)/x
    Poly gk = polyPow(g, k, m);
    for (int i = k; i <= N; i++) res[i] = gk[i - k] * cb.f[i] * cb.fi[k];
    return res;
}
Poly bellSeries(int N, const Comb& cb) {
    Poly e(N + 1, 0); for (int i = 1; i <= N; i++) e[i] = cb.fi[i];
    Poly b = polyExp(e, N + 1); for (int i = 0; i <= N; i++) b[i] *= cb.f[i];
    return b;
}
Poly bernoulli(int N, const Comb& cb) {                     // cb covers N + 1
    Poly g(N + 1); for (int i = 0; i <= N; i++) g[i] = cb.fi[i + 1];
    Poly b = polyInv(g, N + 1); for (int i = 0; i <= N; i++) b[i] *= cb.f[i];
    return b;
}
Poly partitionSeries(int N) {
    Poly q(N + 1, 0);
    for (int k = 0;; k++) {                                 // generalised pentagonal numbers for k, -k
        bool any = false;
        for (int s : {k, -k}) { if (k == 0 && s < 0) continue; ll e = (ll)s * (3 * s - 1) / 2; if (e > N) continue; any = true; q[e] = k % 2 ? mint(-1) : mint(1); }
        if (!any) break;
    }
    return polyInv(q, N + 1);
}
Poly powerSums(ll n, int K, const Comb& cb) {               // ps[k] = sum_{i<n} i^k, k = 0..K (0^0 = 1); cb covers K + 1
    Poly num(K + 1), den(K + 1); mint nn = mint(n), p = nn;
    for (int i = 0; i <= K; i++) num[i] = p * cb.fi[i + 1], den[i] = cb.fi[i + 1], p *= nn;
    Poly s = mulTrunc(num, polyInv(den, K + 1), K + 1);
    for (int i = 0; i <= K; i++) s[i] *= cb.f[i];
    return s;
}
Poly eulerTransform(const vector<mint>& a, int N, bool distinct) {   // a[k] objects of size k (k >= 1); multisets or sets, sizes 0..N
    vector<mint> iv(N + 2, 1); for (int i = 2; i <= N + 1; i++) iv[i] = -mint(998244353 / i) * iv[998244353 % i];
    Poly lg(N + 1, 0);
    for (int k = 1; k < sz(a) && k <= N; k++) if (a[k] != mint(0))
        for (int j = 1; (ll)k * j <= N; j++) { mint t = a[k] * iv[j]; lg[k * j] += distinct && j % 2 == 0 ? -t : t; }
    return polyExp(lg, N + 1);
}
Poly subsetSumCounts(const vector<int>& a, int t) {          // #subsets with each sum 0..t (a_i >= 1)
    vector<mint> cnt(t + 1, 0); for (int x : a) if (x <= t) cnt[x] += 1;
    return eulerTransform(cnt, t, true);
}
Poly boundedKnapsackCounts(const vector<int>& cnt, const vector<int>& lim, int t) {   // cnt[v] kinds of value v, each used 0..lim[v] times
    vector<mint> iv(t + 2, 1); for (int i = 2; i <= t + 1; i++) iv[i] = -mint(998244353 / i) * iv[998244353 % i];
    Poly lg(t + 1, 0);
    for (int v = 1; v < sz(cnt) && v <= t; v++) {
        if (!cnt[v]) continue;
        for (int j = 1; (ll)v * j <= t; j++) lg[v * j] += mint(cnt[v]) * iv[j];   // -log(1 - x^v)
        ll w = (ll)v * (lim[v] + 1);
        for (int j = 1; w * j <= t; j++) lg[w * j] -= mint(cnt[v]) * iv[j];       // + log(1 - x^w)
    }
    return polyExp(lg, t + 1);
}

Proof: the counting series

Stirling first kind: $x^{\overline{2m}}=x^{\overline{m}}(x+m{)}^{\overline{m}}$ and $(x+m{)}^{\overline{m}}=F_m(x+m)$, a Taylor shift; an odd $n$ multiplies by $x+n-1$. Each level costs $O(n\log n)$ and the sizes halve: $O(n\log n)$ total. Second kind: inclusion–exclusion $S(n,k)=\frac{1}{k!}{\sum }_j(-1{)}^{k-j}(\genfrac{}{}{0ex}{}{k}{j})j^n$ is the stated convolution. Columns: a set partition into $k$ blocks is a set of $k$ non-empty sets: EGF $\frac{(e^x-1{)}^k}{k!}$; a permutation with $k$ cycles is a set of $k$ cycles, EGF of one cycle ${\sum }_n(n-1)!\frac{x^n}{n!}=-\ln (1-x)$. Factoring $x^k$ out keeps the power well defined. Bell: a partition is a set of non-empty blocks, $\exp (e^x-1)$. Bernoulli: defined by $\frac{x}{e^x-1}=\sum B_n\frac{x^n}{n!}$, and $\frac{e^x-1}{x}=\sum \frac{x^i}{(i+1)!}$. Partitions: Euler’s pentagonal theorem ${\prod }_k(1-x^k)={\sum }_j(-1{)}^j{x}^{j(3j-1)/2}$ (proof in Stirling, Bell and Partitions). Power sums: ${\sum }_{i<n}{e}^{ix}=\frac{e^{nx}-1}{e^x-1}$, and $k![x^k]$ of it is ${\sum }_{i<n}{i}^k$; both numerator and denominator are divided by $x$. Euler / set transform: $\ln \frac{1}{1-x^k}={\sum }_j\frac{x^{jk}}{j}$ and $\ln (1+x^k)={\sum }_j\frac{(-1{)}^{j+1}{x}^{jk}}{j}$; the harmonic double loop has ${\sum }_{k}N/k=O(N\log N)$ terms. Bounded knapsack: $\frac{1-x^{v(L+1)}}{1-x^v}=1+x^v+\cdots +x^{vL}$. Each function was checked against the $O(n^2)$ tables, brute-force subset enumeration and direct power sums.

Sparse series

When $f$ has only $k$ non-zero terms, $1/f$, $\exp f$, $\log f$, $f^M$ and $\sqrt{f}$ cost $O(nk)$ without NTT (any prime modulus): each solves a first-order linear differential equation coefficient by coefficient ($fg=1$; $g^{\prime }=f^{\prime }g$; $f{g}^{\prime }=f^{\prime }$; $f{g}^{\prime }=M{f}^{\prime }g$). For $f^M$, factor out the lowest term $c{x}^s$: $f^M=c^M{x}^{sM}(f/(c{x}^s){)}^M$, zero if $sM\ge n$; $\sqrt{f}$ is $M=\frac{1}{2}$ with a Tonelli–Shanks root of $c$ and needs even $s$ (Library Checker ”… (Sparse)” problems).

using Sparse = vector<pair<int, mint>>;                     // (exponent, coefficient)
Sparse normalizeSparse(Sparse f) {                          // sort, merge equal exponents, drop zeros
    sort(all(f), [](auto& a, auto& b) { return a.first < b.first; }); Sparse r;
    for (auto [e, c] : f) { if (!r.empty() && r.back().first == e) r.back().second += c; else r.push_back({e, c}); }
    Sparse out; for (auto& pr : r) if (pr.second != mint(0)) out.push_back(pr);
    return out;
}
vector<mint> invTable(int n) { vector<mint> iv(n + 1, 1); for (int i = 2; i <= n; i++) iv[i] = -mint(998244353 / i) * iv[998244353 % i]; return iv; }
vector<mint> sparseInv(const Sparse& f, int n) {            // f_0 != 0
    mint f0 = 0; for (auto [i, c] : f) if (i == 0) f0 += c;
    mint i0 = f0.inv(); vector<mint> g(n, 0);
    for (int k = 0; k < n; k++) { mint s = k == 0 ? mint(1) : mint(0); for (auto [i, c] : f) if (i >= 1 && i <= k) s -= c * g[k - i]; g[k] = s * i0; }
    return g;
}
vector<mint> sparseExp(const Sparse& f, int n) {            // f_0 = 0: k g_k = sum i f_i g_{k-i}
    vector<mint> g(n, 0), iv = invTable(n); if (n) g[0] = 1;
    for (int k = 1; k < n; k++) { mint s = 0; for (auto [i, c] : f) if (i >= 1 && i <= k) s += mint(i) * c * g[k - i]; g[k] = s * iv[k]; }
    return g;
}
vector<mint> sparseLog(const Sparse& f, int n) {            // f_0 = 1: k g_k = k f_k - sum_{i>=1} f_i (k - i) g_{k-i}
    vector<mint> g(n, 0), iv = invTable(n), fk(n, 0);
    for (auto [i, c] : f) if (i < n) fk[i] += c;
    for (int k = 1; k < n; k++) { mint s = mint(k) * fk[k]; for (auto [i, c] : f) if (i >= 1 && i < k) s -= c * mint(k - i) * g[k - i]; g[k] = s * iv[k]; }
    return g;
}
vector<mint> sparsePowUnit(const Sparse& f, mint M, mint f0pow, int n) {   // f_0 != 0, M a field element (also 1/2); g_0 = f0pow
    mint f0 = 0; for (auto [i, c] : f) if (i == 0) f0 += c;
    vector<mint> g(n, 0), iv = invTable(n); if (n) g[0] = f0pow; mint i0 = f0.inv();
    for (int k = 1; k < n; k++) {                           // k f_0 g_k = sum_{i>=1} f_i (M i - (k - i)) g_{k-i}
        mint s = 0; for (auto [i, c] : f) if (i >= 1 && i <= k) s += c * (M * mint(i) - mint(k - i)) * g[k - i];
        g[k] = s * iv[k] * i0;
    }
    return g;
}
vector<mint> sparsePow(Sparse f, ll M, int n) {             // f^M, any f, M >= 0
    vector<mint> zero(n, 0);
    if (M == 0) { if (n) zero[0] = 1; return zero; }
    f = normalizeSparse(f);
    if (f.empty()) return zero;
    int s = f[0].first; if ((i128)s * M >= n) return zero;
    mint c = f[0].second; for (auto& [i, x] : f) i -= s;
    vector<mint> g = sparsePowUnit(f, mint(M % 998244353), c.pow(M), n - s * M);
    g.insert(g.begin(), s * M, mint(0)); return g;
}
bool sparseSqrt(Sparse f, int n, vector<mint>& out) {       // false if no square root
    f = normalizeSparse(f);
    out.assign(n, 0);
    if (f.empty()) return true;
    int s = f[0].first; if (s >= n) return true;           // 0 works modulo x^n
    if (s % 2) return false;
    mint r; if (!sqrtMint(f[0].second, r)) return false;
    for (auto& [i, x] : f) i -= s;
    vector<mint> g = sparsePowUnit(f, mint(2).inv(), r, n - s / 2);
    for (int i = 0; i + s / 2 < n; i++) out[i + s / 2] = g[i];
    return true;
}

Proof: the sparse recurrences

Compare coefficients of $x^{k-1}$ in each differential equation. Inverse: $[x^k](fg)=[k=0]$ gives $f_0{g}_k=[k=0]-{\sum }_{i\ge 1}{f}_i{g}_{k-i}$. Exp: $g=e^f$ satisfies $g^{\prime }=f^{\prime }g$: $k{g}_k={\sum }_{i}i{f}_i{g}_{k-i}$. Log: $g=\log f$ satisfies $f{g}^{\prime }=f^{\prime }$: ${\sum }_i{f}_i(k-i)g_{k-i}=k{f}_k$, isolate $i=0$. Power: $g=f^M$ satisfies $f{g}^{\prime }=M{f}^{\prime }g$: ${\sum }_i{f}_i(k-i)g_{k-i}=M{\sum }_{i}i{f}_i{g}_{k-i}$, giving $k{f}_0{g}_k={\sum }_{i\ge 1}{f}_i(Mi-(k-i))g_{k-i}$. This recurrence has a unique solution with prescribed $g_0$, and for $f_0=1$, $g_0=1$ the true power satisfies it, also for $M=\frac{1}{2}$ (then $g^2=f$ since $(g^2{)}^{\prime }f=2g{g}^{\prime }f=g^2{f}^{\prime }$ and both have constant term 1, so $g^2/f$ has zero derivative and equals 1). Scaling by $c=f_0$: $(cu{)}^M=c^M{u}^M$. The exponent enters the recurrence only modulo $p$, while $c^M$ uses the true $M$. Each coefficient touches $k$ terms: $O(nk)$. Verified against the dense NTT versions.

Sums of exponential times polynomial

${\sum }_{i<n}{a}^{i}f(i)$ for a polynomial $f$ of degree $\le d$ given by $f(\mathrm{0..}d)$ and $n\le 10^{18}$ (Library Checker “Sum of Exponential times Polynomial”, and its “Limit” version ${\sum }_{i\ge 0}{a}^{i}f(i)$ for $a\ne 1$), plus the binomial basis $f(x)={\sum }_k{c}_k\left(\genfrac{}{}{0ex}{}{x}{k}\right)$ with $c_k={\Delta }^{k}f(0)$ (Newton forward differences). For $a=1$ the sum is a polynomial of degree $d+1$ (interpolate); for $a\ne 1$, $S(N)={\sum }_{i\le N}{a}^{i}f(i)=a^{N}g(N)+c$ with $\mathrm{deg}g\le d$, $c=\frac{{\sum }_{i=0}^{d+1}\left(\genfrac{}{}{0ex}{}{d+1}{i}\right)(-a{)}^{d+1-i}S(i)}{(1-a{)}^{d+1}}$, $g(i)=(S(i)-c)a^{-i}$, and Lagrange gives $g(N)$.

mint sumExpTimesPoly(vector<mint> fv, mint a, ll n, const Comb& cb) {   // sum_{i<n} a^i f(i), fv = f(0..d); cb covers d + 3
    if (n <= 0) return 0;
    int d = sz(fv) - 1; fv.push_back(lagrangeConsecutive(fv, d + 1, cb));   // f(d + 1)
    vector<mint> S(d + 2); mint ap = 1;
    for (int i = 0; i <= d + 1; i++) S[i] = (i ? S[i - 1] : mint(0)) + ap * fv[i], ap *= a;
    ll N = n - 1;
    if (N <= d + 1) return S[N];
    if (a == mint(0)) return fv[0];
    if (a == mint(1)) return lagrangeConsecutive(S, N, cb);
    mint c = 0, ak = 1;
    for (int i = d + 1; i >= 0; i--) { mint t = cb.C(d + 1, i) * ak * S[i]; c += (d + 1 - i) % 2 ? -t : t; ak *= a; }
    c /= (mint(1) - a).pow(d + 1);
    vector<mint> g(d + 1); mint ainv = a.inv(), aip = 1;
    for (int i = 0; i <= d; i++) g[i] = (S[i] - c) * aip, aip *= ainv;
    return lagrangeConsecutive(g, N, cb) * a.pow(N) + c;
}
mint sumExpTimesPolyLimit(vector<mint> fv, mint a, const Comb& cb) {   // sum_{i>=0} a^i f(i) as the formal limit, a != 1
    if (fv.empty()) return 0;
    int d = sz(fv) - 1; fv.push_back(lagrangeConsecutive(fv, d + 1, cb));
    vector<mint> S(d + 2); mint ap = 1, s = 0;
    for (int i = 0; i <= d + 1; i++) s += ap * fv[i], S[i] = s, ap *= a;
    mint c = 0, ak = 1;
    for (int i = d + 1; i >= 0; i--) { mint t = cb.C(d + 1, i) * ak * S[i]; c += (d + 1 - i) % 2 ? -t : t; ak *= a; }
    return c / (mint(1) - a).pow(d + 1);
}
Poly toBinomialBasis(const vector<mint>& fv, const Comb& cb) {   // c_k = k! sum_j f(j)/j! (-1)^(k-j)/(k-j)!
    int n = sz(fv); Poly a(n), b(n);
    for (int i = 0; i < n; i++) a[i] = fv[i] * cb.fi[i], b[i] = i % 2 ? -cb.fi[i] : cb.fi[i];
    Poly c = polyMul(a, b); c.resize(n); for (int k = 0; k < n; k++) c[k] *= cb.f[k];
    return c;
}
Poly fromBinomialBasis(const vector<mint>& c, const Comb& cb) {  // f(i) = i! sum_k c_k/k! 1/(i-k)!
    int n = sz(c); Poly a(n), b(n);
    for (int i = 0; i < n; i++) a[i] = c[i] * cb.fi[i], b[i] = cb.fi[i];
    Poly f = polyMul(a, b); f.resize(n); for (int i = 0; i < n; i++) f[i] *= cb.f[i];
    return f;
}

Proof: the structure $S(N)=a^{N}g(N)+c$ and the binomial basis

We need a polynomial $g$ with $a^{N}g(N)-a^{N-1}g(N-1)=a^{N}f(N)$, i.e. $g(N)-a^{-1}g(N-1)=f(N)$. On polynomials of degree $\le d$ the map $g\mapsto g-a^{-1}g(\cdot -1)$ is triangular with respect to degree with diagonal $1-a^{-1}\ne 0$, hence bijective, so such $g$ exists and $S(N)-a^{N}g(N)$ is constant: $S(N)=a^{N}g(N)+c$. Since $\mathrm{deg}g\le d$, ${\Delta }^{d+1}g(0)={\sum }_i\left(\genfrac{}{}{0ex}{}{d+1}{i}\right)(-1{)}^{d+1-i}g(i)=0$; substituting $g(i)=(S(i)-c)a^{-i}$ and multiplying by $a^{d+1}$: ${\sum }_i\left(\genfrac{}{}{0ex}{}{d+1}{i}\right)(-a{)}^{d+1-i}S(i)=c{\sum }_i(\genfrac{}{}{0ex}{}{d+1}{i})(-a{)}^{d+1-i}=c(1-a{)}^{d+1}$. The “limit” is this $c$ (the value of ${\sum }_i{a}^{i}f(i)$ whenever it converges, $a^{N}g(N)\to 0$). For $a=1$, $S$ is a polynomial of degree $d+1$ (finite differences). Binomial basis: $\left(\genfrac{}{}{0ex}{}{x}{k}\right)$ is a basis of polynomials and $\Delta \left(\genfrac{}{}{0ex}{}{x}{k}\right)=\left(\genfrac{}{}{0ex}{}{x}{k-1}\right)$, so $c_k={\Delta }^{k}f(0)={\sum }_j(-1{)}^{k-j}(\genfrac{}{}{0ex}{}{k}{j})f(j)$, which is the convolution of $f(j)/j!$ with $(-1{)}^m/m!$ times $k!$; the inverse map $f(i)={\sum }_k{c}_k\left(\genfrac{}{}{0ex}{}{i}{k}\right)$ is the convolution with $1/m!$. Verified against direct summation (including $a=0,1$) and round trips.

Related notes

• Stirling, Bell and Partitions
• Generating Functions
• Knapsack (counting versions via exp / log)