iv. Online Convolution

Online (relaxed) convolution returns $h_i={\sum }_{j\le i}{f}_j{g}_{i-j}$ as soon as $f_i$ and $g_i$ are known, so $f$ and $g$ may depend on $h_0..h_{i-1}$: the tool for sequences defined by a convolution with their own prefix ($f_n={\sum }_{j<n}{f}_j{g}_{n-j}$, rooted-tree counts, $\exp$ / $\log$ of something defined implicitly). Cost $O(n{\log }^{2}n)$ ($n=5\cdot 10^5$ in about 0.6 s), versus $O(n^2)$ naively. Uses ntt/convolution/Poly from NTT.

Setting	Method	Cost
$f_n={\sum }_{j<n}{f}_j{g}_{n-j}$, $g$ known in advance (semi-relaxed)	CDQ divide and conquer: finish the left half, add its contribution to the right half with one convolution, recurse	$O(n{\log }^{2}n)$
$f$ and $g$ both depend on the prefix of $h$ (fully relaxed)	RelaxedConvolution (blocks of size $2^z$ along the anti-diagonals)	$O(n{\log }^{2}n)$
the equation is algebraic in a known series ($\exp$, inverse, log of a known $A$)	Newton (Power Series)	$O(n\log n)$

Poly cdqSelfConvolution(const Poly& g, mint f0, int n) {    // f_0 = f0, f_i = sum_{j<i} f_j g_{i-j}  (i.e. f = f0 / (1 - g) with g_0 ignored)
    Poly f(n, 0); if (n == 0) return f; f[0] = f0;
    auto rec = [&](auto&& self, int l, int r) -> void {      // on entry f[l] is final; contributions from [0, l) already added to f[l, r)
        if (r - l == 1) return;
        int m = (l + r) / 2; self(self, l, m);
        Poly a(f.begin() + l, f.begin() + m), b(min(r - l, sz(g)), 0);
        for (int i = 1; i < sz(b); i++) b[i] = g[i];
        Poly c = polyMul(a, b);
        for (int i = m; i < r; i++) if (i - l < sz(c)) f[i] += c[i - l];
        self(self, m, r);
    };
    rec(rec, 0, n);
    return f;
}
 
template <int P> struct RelaxedConvolution {
    using Z = Mint<P>; static constexpr int NAIVE = 32;
    vector<Z> f, g, h; vector<vector<Z>> fBand, gBand;     // NTT (size 2s) of f[s-1, 2s-1) and g[s-1, 2s-1), s = 2^z
    Z push(Z a, Z b) {                                       // h_p = sum_{j + k = p} f_j g_k, with f_p = a, g_p = b
        int p = sz(f); f.push_back(a); g.push_back(b);
        if (sz(h) < 2 * p + 2) h.resize(2 * p + 2);
        for (int z = 0, s = 1; 2 * s <= p + 2 && (p + 2) % s == 0; z++, s *= 2) {
            int x = p + 1 - s; bool diag = x == s - 1;       // the square [s-1, 2s-1)^2, counted once
            if (s <= NAIVE) {
                for (int j = x; j < x + s; j++) for (int k = s - 1; k < 2 * s - 1; k++) h[j + k] += diag ? f[j] * g[k] : f[j] * g[k] + g[j] * f[k];
                continue;
            }
            vector<Z> u(f.begin() + x, f.begin() + x + s), v(g.begin() + x, g.begin() + x + s);
            u.resize(2 * s); v.resize(2 * s); ntt(u); ntt(v);
            if (diag) { fBand.resize(z + 1); gBand.resize(z + 1); fBand[z] = u; gBand[z] = v; for (int i = 0; i < 2 * s; i++) u[i] *= v[i]; }
            else for (int i = 0; i < 2 * s; i++) u[i] = u[i] * gBand[z][i] + v[i] * fBand[z][i];
            ntt(u, true);
            for (int i = 0; i < 2 * s - 1; i++) h[p + i] += u[i];
        }
        return h[p];
    }
};
 
vector<mint> rootedTrees(int n) {                           // unlabelled rooted trees a_1..a_n (index = size), relaxed
    vector<mint> a(n + 1, 0), b(n + 1, 0), iv(n + 1, 1); if (n >= 1) a[1] = 1;
    for (int i = 2; i <= n; i++) iv[i] = -mint(998244353 / i) * iv[998244353 % i];
    RelaxedConvolution<998244353> rc;                       // h_m = sum_{k=1}^{m} b_k a_{m+1-k}: push (b_{i+1}, a_{i+1})
    for (int m = 1; m < n; m++) {
        for (int d = 1; d * d <= m; d++) if (m % d == 0) { b[m] += mint(d) * a[d]; if (d * d != m) b[m] += mint(m / d) * a[m / d]; }
        mint h = rc.push(b[m], a[m]);                        // after push number m - 1 (0-indexed): sum_{j+k = m-1} b_{j+1} a_{k+1}
        a[m + 1] = h * iv[m];
    }
    return a;
}

Rooted trees. $r(x)=x{\prod }_k(1-x^k{)}^{-a_k}$ gives $a_{m+1}=\frac{1}{m}{\sum }_{k=1}^m{b}_k{a}_{m+1-k}$ with $b_k={\sum }_{d\mid k}d{a}_d$ (Graph Counting): $b_m$ needs $a_m$, and $a_{m+1}$ needs the convolution up to $m$, a relaxed product of the shifted sequences $(b_1,b_2,\dots )$ and $(a_1,a_2,\dots )$.

Proof: CDQ and the relaxed block schedule

CDQ. Invariant: when rec(l, r) is entered, f[l] is final and every contribution $f_j{g}_{i-j}$ with $j<l$, $i\in [l,r)$ has been added. The left recursion keeps the invariant for $[l,m)$ (contributions from $j<l$ are already in, those from $j\in [l,i)$ are added inside), after which f[l..m) are final. The convolution then adds exactly the contributions with $j\in [l,m)$ to every $i\in [m,r)$ (the kernel uses $g_1..g_{r-l-1}$, which covers $i-j\le r-l-1$), restoring the invariant for $[m,r)$. Each level does convolutions of total length $O(n)$: $O(n{\log }^{2}n)$. With $g_0$ excluded this computes $f=f_0+f\cdot (g-g_0)$, i.e. $f=f_0/(1-g+g_0)$.

Relaxed schedule. Partition the index pairs $(j,k)$ into squares: the diagonal squares $[s-1,2s-1{)}^2$ and, for every $s=2^z$, the off-diagonal blocks $[x,x+s)\times [s-1,2s-1)$ (and their mirror images) with $x+1\equiv 0(\mathrm{mod}s)$, $x\ge 2s-1$. These tile the quadrant $\{j,k\ge 0\}$ exactly once (a dyadic decomposition: each pair lies in the block determined by the largest power of two for which both coordinates fit). The block is processed at step $p=x+s-1$, when $f,g$ are known up to index $p$: its smallest sum is $x+s-1=p$, so its contribution is added before $h_p$ is returned, and it uses only indices $\le p$. The band transforms of $f,g$ on $[s-1,2s-1)$ are computed at the diagonal step and reused. Per $s$ there are $O(n/s)$ blocks costing $O(s\log s)$, so $O(n\log n)$ per level and $O(n{\log }^{2}n)$ in total. Both methods were checked against the naive $O(n^2)$ recurrences (including $f$ depending on $h$), against polyInv for $f=1/(1-g)$, and the rooted-tree counts against OEIS A000081.

Related notes

• CDQ Divide and Conquer
• Counting with Power Series
• Power Series