2223FUnreviewedunratedGPT-5.4 (high)

Zhily and Colorful Strings

Progressive hints first, then the full explanation and implementation when you're ready to cash out.

divide and conquer fft math

Original problem

Review status

AI-generated and still unreviewed. Double-check the details before internalizing them.

Hints

Progressive nudges

Open only as much as you need to keep the solve alive.

First, stop thinking about deletions as an algorithm. Think algebraically: for each type-color pair, you have a generator $g$ with relation $g^{d_i}=1$ . A colored string is valid iff the corresponding word becomes the identity.

For a single generator of order $d$ , the only valid positive lengths are multiples of $d$ , so its moment series is $M_d(z)=\frac{1}{1-z^d}$ . The useful transform here is the noncrossing-partition / free-cumulant series $C$ defined by $M(z)=C(zM(z))$ .

Solve that transform for one generator: if $C_d$ is the cumulant series, then $C_d(w)^d-C_d(w)^{d-1}=w^d.$ Now set $A_d(x)=C_d(x^{1/d})-1$ . Then $A_d$ is exactly the inverse series of $u(1+u)^{d-1}$ , i.e. $A_d(x)(1+A_d(x))^{d-1}=x.$

Let $F(\mathbf z)$ be the generating function of all valid strings, where $z_i$ marks one character of type $i$ . Since different types/colors are free, their cumulants add, and you get $F=1+\sum_i c_i\,A_{d_i}(z_i^{d_i}F^{d_i}).$ Immediate cheap shot: if some $d_i\nmid n_i$ , the answer is $0$ .

Now do the killer substitution $z_i^{d_i}=p_i z^{d_i}$ and $U=zF$ . Then $U=z\Phi(U),\qquad \Phi(u)=1+\sum_i c_i A_{d_i}(p_i u^{d_i}),$ so ordinary Lagrange inversion applies: $\text{Ans}=\frac1{D+1}[p_1^{k_1}\cdots p_m^{k_m}]\Bigl(1+\sum_i c_iA_{d_i}(p_i)\Bigr)^{D+1},$ where $k_i=n_i/d_i$ and $D=\sum_i n_i$ . Use $(1+\sum B_i)^{D+1}=(D+1)!\,[t^{D+1}]e^t\prod_i e^{tB_i}$ and define $E_i(t)=[q^{k_i}]e^{t c_i A_{d_i}(q)}$ . Then multiply all $E_i$ by NTT. Their coefficients are $[t^j]E_i(t)=\frac{c_i^j}{k_i(j-1)!}(-1)^{k_i-j}\binom{d_i k_i-j-1}{k_i-j},$ and in code it is best generated by the backward recurrence from $e_{i,k_i}=c_i^{k_i}/k_i!$ .

Zhily disguised a free-probability problem as a string problem. Rude. Also elegant.

1. Turn the string into a group word

For every type-color pair $(i,\alpha)$ , create a generator $g_{i,\alpha}$ with relation $g_{i,\alpha}^{d_i}=e.$ So each color of type $i$ is a cyclic generator of order $d_i$ .

Now read a colored string as a positive word in the free product $\ast_{i=1}^m \ast_{\alpha=1}^{c_i} C_{d_i}.$ Deleting $d_i$ consecutive equal letters is exactly using the relation $g_{i,\alpha}^{d_i}=e$ .

So a colored string is valid iff the corresponding word equals the identity. In other words, we are counting positive words with prescribed numbers of letters of each type whose reduced normal form is empty.

A first immediate consequence:

if $d_i \nmid n_i$ for some $i$ , the answer is $0$ .

Because the total number of letters of type $i$ must be a multiple of $d_i$ .

Let $k_i=\frac{n_i}{d_i}, \qquad D=\sum_i n_i=\sum_i d_i k_i.$

2. One generator: moments and cumulants

Take one generator $g$ of order $d$ . Among positive words over just $g$ , the identity appears exactly for lengths divisible by $d$ , so the moment series is $M_d(z)=1+z^d+z^{2d}+\cdots=\frac{1}{1-z^d}.$

Now introduce the standard noncrossing-partition transform $C_d$ , defined by $M_d(z)=C_d\bigl(zM_d(z)\bigr).$ If you know free probability, this is the free-cumulant series. If you do not, that is fine; we only use this identity and the fact that these cumulant series add for free sums.

Let $w=zM_d(z)$ . Then $M_d(z)=C_d(w)$ and $z=w/C_d(w)$ . Plugging into $M_d(z)=1/(1-z^d)$ gives $C_d(w)^d-C_d(w)^{d-1}=w^d.$

Define $A_d(x)=C_d\bigl(x^{1/d}\bigr)-1.$ Then $A_d(0)=0$ , and the previous equation becomes $A_d(x)(1+A_d(x))^{d-1}=x.$ So $A_d$ is the compositional inverse of $u(1+u)^{d-1}$ .

That inverse series is the whole game.

3. All types together

For type $i$ , let $X_i=\sum_{\alpha=1}^{c_i} g_{i,\alpha}.$ Expanding $X_i^{n_i}$ chooses colors for the letters of type $i$ .

Now consider the multivariate generating series $F(z_1,\dots,z_m)=\sum a_{n_1,\dots,n_m} z_1^{n_1}\cdots z_m^{n_m},$ where $a_{n_1,\dots,n_m}$ is exactly the answer we want.

Because different types/colors live in different free factors, their cumulants add. Therefore the cumulant series of the sum $\sum_i z_i X_i$ is just the sum of the type contributions, and the moment-cumulant relation gives $F=1+\sum_{i=1}^m c_i\,A_{d_i}(z_i^{d_i}F^{d_i}).$

This is already a correct functional equation, but in this form it looks like something that wants to punch you in the face.

4. The substitution that saves the day

We only need the coefficient for fixed $k_i=n_i/d_i$ . So write $z_i^{d_i}=p_i z^{d_i}.$ Then the coefficient of $p_1^{k_1}\cdots p_m^{k_m} z^D$ in $F$ is exactly the desired answer.

Now set $U=zF.$ Then the equation becomes $U=z\Phi(U), \qquad \Phi(u)=1+\sum_{i=1}^m c_i A_{d_i}(p_i u^{d_i}).$

Beautiful. Now it is an ordinary one-variable Lagrange inversion setup.

By Lagrange:

=\frac{1}{D+1}[u^D p_1^{k_1}\cdots p_m^{k_m}]\Phi(u)^{D+1}.$$ Since $U=zF$, this is exactly our answer. But every monomial in $A_{d_i}(p_i u^{d_i})$ is of the form $(p_i u^{d_i})^r$. So once the degree in each $p_i$ is fixed to $k_i$, the degree in $u$ is **automatically** $$\sum_i d_i k_i=D.$$ Therefore the $u^D$ extraction is forced and we can drop it: $$\boxed{\text{Ans}=\frac{1}{D+1}[p_1^{k_1}\cdots p_m^{k_m}]\left(1+\sum_{i=1}^m c_i A_{d_i}(p_i)\right)^{D+1}}.$$ That is the key formula. ## 5. Split the multivariate coefficient into one polynomial per type Use the exponential trick $$X^{D+1}=(D+1)!\,[t^{D+1}]e^{tX}.$$ Apply it with $$X=1+\sum_i c_i A_{d_i}(p_i).$$ Then $$\text{Ans}=D!\,[t^{D+1}]e^t\prod_{i=1}^m E_i(t),$$ where $$E_i(t)=[q^{k_i}]\exp\bigl(t c_i A_{d_i}(q)\bigr).$$ Now the monster has been chopped into $m$ ordinary one-variable polynomials. That is where the FFT/NTT tag finally stops trolling and starts making sense. ## 6. Explicit formula for $E_i(t)$ We still need coefficients of $E_i$. Since $A_d$ satisfies $$A_d(q)=q(1+A_d(q))^{1-d},$$ Lagrange again gives, for any series $G$, $$[q^k]G(A_d(q))=\frac1k[u^{k-1}]G'(u)(1+u)^{-(d-1)k}.$$ Take $$G(u)=e^{tcu}.$$ Then $$E_i(t)=\frac{t c_i}{k_i}[u^{k_i-1}]e^{t c_i u}(1+u)^{-(d_i-1)k_i}.$$ Extracting the coefficient of $t^j$ gives $$[t^j]E_i(t)=\frac{c_i^j}{k_i(j-1)!}(-1)^{k_i-j}\binom{d_i k_i-j-1}{k_i-j} \qquad (1\le j\le k_i).$$ In code, do **not** recompute binomials from scratch like a maniac. Use the backward recurrence: $$e_{i,k_i}=\frac{c_i^{k_i}}{k_i!},$$ $$e_{i,j}=-e_{i,j+1}\cdot \frac{j\,(d_i k_i-j-1)}{c_i\,(k_i-j)} \qquad (j=k_i-1,\dots,1).$$ This builds each polynomial in $O(k_i)$. ## 7. Final combination Let $$P(t)=\prod_{i=1}^m E_i(t)=\sum_{j\ge 0} P_j t^j.$$ Then $$\text{Ans}=D!\,[t^{D+1}]e^t P(t) =\sum_{j\ge 0} P_j\,\frac{D!}{(D+1-j)!}.$$ Since $\deg P=\sum_i k_i=:K \le 2\cdot 10^5$, we only need $j\le K$. So compute the falling factorials iteratively: $$w_1=1, \qquad w_{j+1}=w_j\,(D-j+1).$$ Then $$\text{Ans}=\sum_{j=1}^{K} P_j\,w_j.$$ Even though $D$ itself can be huge, $K<\text{MOD}$, so this product can be done directly modulo $998244353$; if one factor hits $0$, all later ones stay $0$ and life gets simpler. ## 8. Complexity Let $$K=\sum_i k_i=\sum_i \frac{n_i}{d_i}.$$ The statement guarantees total $K$ over all test cases is at most $2\cdot 10^5$. - Building all $E_i$: $O(K)$. - Multiplying them with balanced NTT merges: $O(K\log^2 K)$. - Final weighted sum: $O(K)$. So the whole thing is comfortably fast. Yeah, the problem is basically: **free probability + two Lagrange inversions + NTT**. Completely normal Codeforces behavior. Totally not cursed.

#include <bits/stdc++.h>
using namespace std;

using ll = long long;

void setIO() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
}

const int MOD = 998244353;
const int G = 3;
const int MAXK = 200000 + 5;

ll modpow(ll a, ll e) {
    ll r = 1;
    while (e > 0) {
        if (e & 1) r = r * a % MOD;
        a = a * a % MOD;
        e >>= 1;
    }
    return r;
}

vector<int> invv(MAXK), fact(MAXK), invfact(MAXK);

void initComb() {
    fact[0] = 1;
    for (int i = 1; i < MAXK; ++i) fact[i] = (ll)fact[i - 1] * i % MOD;
    invfact[MAXK - 1] = modpow(fact[MAXK - 1], MOD - 2);
    for (int i = MAXK - 1; i >= 1; --i) invfact[i - 1] = (ll)invfact[i] * i % MOD;

    invv[1] = 1;
    for (int i = 2; i < MAXK; ++i) {
        invv[i] = MOD - (ll)(MOD / i) * invv[MOD % i] % MOD;
    }
}

void ntt(vector<int>& a, bool invert) {
    int n = (int)a.size();
    if (n == 1) return;

    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 len = 2; len <= n; len <<= 1) {
        ll wlen = modpow(G, (MOD - 1) / len);
        if (invert) wlen = modpow(wlen, MOD - 2);
        for (int i = 0; i < n; i += len) {
            ll w = 1;
            for (int j = 0; j < len / 2; ++j) {
                int u = a[i + j];
                int v = (ll)a[i + j + len / 2] * w % MOD;
                a[i + j] = u + v;
                if (a[i + j] >= MOD) a[i + j] -= MOD;
                a[i + j + len / 2] = u - v;
                if (a[i + j + len / 2] < 0) a[i + j + len / 2] += MOD;
                w = w * wlen % MOD;
            }
        }
    }

    if (invert) {
        ll inv_n = modpow(n, MOD - 2);
        for (int &x : a) x = (ll)x * inv_n % MOD;
    }
}

vector<int> multiply(const vector<int>& a, const vector<int>& b) {
    if (a.empty() || b.empty()) return {};
    int n = (int)a.size(), m = (int)b.size();
    if ((ll)n * m <= 4000) {
        vector<int> res(n + m - 1);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                res[i + j] = (res[i + j] + (ll)a[i] * b[j]) % MOD;
            }
        }
        return res;
    }

    int sz = 1;
    while (sz < n + m - 1) sz <<= 1;
    vector<int> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    fa.resize(sz);
    fb.resize(sz);

    ntt(fa, false);
    ntt(fb, false);
    for (int i = 0; i < sz; ++i) fa[i] = (ll)fa[i] * fb[i] % MOD;
    ntt(fa, true);
    fa.resize(n + m - 1);
    return fa;
}

vector<int> build_poly(int n, int c, int k) {
    vector<int> poly(k + 1);
    ll invc = modpow(c, MOD - 2);
    poly[k] = (ll)modpow(c, k) * invfact[k] % MOD;
    for (int j = k - 1; j >= 1; --j) {
        ll val = poly[j + 1];
        val = (MOD - val) % MOD;
        val = val * j % MOD;
        val = val * (n - j - 1LL) % MOD;
        val = val * invc % MOD;
        val = val * invv[k - j] % MOD;
        poly[j] = (int)val;
    }
    return poly;
}

int main() {
    setIO();
    initComb();

    int T;
    cin >> T;
    while (T--) {
        int m;
        cin >> m;

        vector<tuple<int, int, int>> items;
        items.reserve(m);

        bool bad = false;
        ll D = 0;
        int K = 0;

        for (int i = 0; i < m; ++i) {
            int n, c, d;
            cin >> n >> c >> d;
            items.emplace_back(n, c, d);
            D += n;
            if (n % d != 0) bad = true;
            else K += n / d;
        }

        if (bad) {
            cout << 0 << '\n';
            continue;
        }

        vector<vector<int>> polys;
        polys.reserve(m * 2 + 5);
        priority_queue<pair<int, int>, vector<pair<int, int>>, greater<pair<int, int>>> pq;

        for (auto [n, c, d] : items) {
            int k = n / d;
            polys.push_back(build_poly(n, c, k));
            int id = (int)polys.size() - 1;
            pq.push({(int)polys[id].size(), id});
        }

        while (pq.size() > 1) {
            auto [sa, ia] = pq.top(); pq.pop();
            auto [sb, ib] = pq.top(); pq.pop();
            polys.push_back(multiply(polys[ia], polys[ib]));
            int id = (int)polys.size() - 1;
            pq.push({(int)polys[id].size(), id});
        }

        vector<int> P = polys[pq.top().second];

        ll ans = 0;
        ll w = 1;              // w_j = D! / (D + 1 - j)!, starts with w_1 = 1
        ll Dmod = D % MOD;

        for (int j = 1; j < (int)P.size(); ++j) {
            ans += w * P[j] % MOD;
            if (ans >= MOD) ans -= MOD;

            ll factor = Dmod - (j - 1);
            factor %= MOD;
            if (factor < 0) factor += MOD;
            w = w * factor % MOD;
            if (w == 0) break;
        }

        cout << ans % MOD << '\n';
    }
}

#include <bits/stdc++.h>
using namespace std;

using ll = long long;

void setIO() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
}

const int MOD = 998244353;
const int G = 3;
const int MAXK = 200000 + 5;

ll modpow(ll a, ll e) {
    ll r = 1;
    while (e > 0) {
        if (e & 1) r = r * a % MOD;
        a = a * a % MOD;
        e >>= 1;
    }
    return r;
}

vector<int> invv(MAXK), fact(MAXK), invfact(MAXK);

void initComb() {
    fact[0] = 1;
    for (int i = 1; i < MAXK; ++i) fact[i] = (ll)fact[i - 1] * i % MOD;
    invfact[MAXK - 1] = modpow(fact[MAXK - 1], MOD - 2);
    for (int i = MAXK - 1; i >= 1; --i) invfact[i - 1] = (ll)invfact[i] * i % MOD;

    invv[1] = 1;
    for (int i = 2; i < MAXK; ++i) {
        invv[i] = MOD - (ll)(MOD / i) * invv[MOD % i] % MOD;
    }
}

void ntt(vector<int>& a, bool invert) {
    int n = (int)a.size();
    if (n == 1) return;

    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 len = 2; len <= n; len <<= 1) {
        ll wlen = modpow(G, (MOD - 1) / len);
        if (invert) wlen = modpow(wlen, MOD - 2);
        for (int i = 0; i < n; i += len) {
            ll w = 1;
            for (int j = 0; j < len / 2; ++j) {
                int u = a[i + j];
                int v = (ll)a[i + j + len / 2] * w % MOD;
                a[i + j] = u + v;
                if (a[i + j] >= MOD) a[i + j] -= MOD;
                a[i + j + len / 2] = u - v;
                if (a[i + j + len / 2] < 0) a[i + j + len / 2] += MOD;
                w = w * wlen % MOD;
            }
        }
    }

    if (invert) {
        ll inv_n = modpow(n, MOD - 2);
        for (int &x : a) x = (ll)x * inv_n % MOD;
    }
}

vector<int> multiply(const vector<int>& a, const vector<int>& b) {
    if (a.empty() || b.empty()) return {};
    int n = (int)a.size(), m = (int)b.size();
    if ((ll)n * m <= 4000) {
        vector<int> res(n + m - 1);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                res[i + j] = (res[i + j] + (ll)a[i] * b[j]) % MOD;
            }
        }
        return res;
    }

    int sz = 1;
    while (sz < n + m - 1) sz <<= 1;
    vector<int> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    fa.resize(sz);
    fb.resize(sz);

    ntt(fa, false);
    ntt(fb, false);
    for (int i = 0; i < sz; ++i) fa[i] = (ll)fa[i] * fb[i] % MOD;
    ntt(fa, true);
    fa.resize(n + m - 1);
    return fa;
}

vector<int> build_poly(int n, int c, int k) {
    vector<int> poly(k + 1);
    ll invc = modpow(c, MOD - 2);
    poly[k] = (ll)modpow(c, k) * invfact[k] % MOD;
    for (int j = k - 1; j >= 1; --j) {
        ll val = poly[j + 1];
        val = (MOD - val) % MOD;
        val = val * j % MOD;
        val = val * (n - j - 1LL) % MOD;
        val = val * invc % MOD;
        val = val * invv[k - j] % MOD;
        poly[j] = (int)val;
    }
    return poly;
}

int main() {
    setIO();
    initComb();

    int T;
    cin >> T;
    while (T--) {
        int m;
        cin >> m;

        vector<tuple<int, int, int>> items;
        items.reserve(m);

        bool bad = false;
        ll D = 0;
        int K = 0;

        for (int i = 0; i < m; ++i) {
            int n, c, d;
            cin >> n >> c >> d;
            items.emplace_back(n, c, d);
            D += n;
            if (n % d != 0) bad = true;
            else K += n / d;
        }

        if (bad) {
            cout << 0 << '\n';
            continue;
        }

        vector<vector<int>> polys;
        polys.reserve(m * 2 + 5);
        priority_queue<pair<int, int>, vector<pair<int, int>>, greater<pair<int, int>>> pq;

        for (auto [n, c, d] : items) {
            int k = n / d;
            polys.push_back(build_poly(n, c, k));
            int id = (int)polys.size() - 1;
            pq.push({(int)polys[id].size(), id});
        }

        while (pq.size() > 1) {
            auto [sa, ia] = pq.top(); pq.pop();
            auto [sb, ib] = pq.top(); pq.pop();
            polys.push_back(multiply(polys[ia], polys[ib]));
            int id = (int)polys.size() - 1;
            pq.push({(int)polys[id].size(), id});
        }

        vector<int> P = polys[pq.top().second];

        ll ans = 0;
        ll w = 1;              // w_j = D! / (D + 1 - j)!, starts with w_1 = 1
        ll Dmod = D % MOD;

        for (int j = 1; j < (int)P.size(); ++j) {
            ans += w * P[j] % MOD;
            if (ans >= MOD) ans -= MOD;

            ll factor = Dmod - (j - 1);
            factor %= MOD;
            if (factor < 0) factor += MOD;
            w = w * factor % MOD;
            if (w == 0) break;
        }

        cout << ans % MOD << '\n';
    }
}

View on Codeforces Back to the list