Search for a command to run...
Progressive hints first, then the full explanation and implementation when you're ready to cash out.
Review status
AI-generated and still unreviewed. Double-check the details before internalizing them.
Hints
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 with relation . A colored string is valid iff the corresponding word becomes the identity.
For a single generator of order , the only valid positive lengths are multiples of , so its moment series is . The useful transform here is the noncrossing-partition / free-cumulant series defined by .
Solve that transform for one generator: if is the cumulant series, then Now set . Then is exactly the inverse series of , i.e.
Let be the generating function of all valid strings, where marks one character of type . Since different types/colors are free, their cumulants add, and you get Immediate cheap shot: if some , the answer is .
Now do the killer substitution and . Then so ordinary Lagrange inversion applies: where and . Use and define . Then multiply all by NTT. Their coefficients are and in code it is best generated by the backward recurrence from .
Zhily disguised a free-probability problem as a string problem. Rude. Also elegant.
For every type-color pair , create a generator with relation So each color of type is a cyclic generator of order .
Now read a colored string as a positive word in the free product Deleting consecutive equal letters is exactly using the relation .
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:
Because the total number of letters of type must be a multiple of .
Let
Take one generator of order . Among positive words over just , the identity appears exactly for lengths divisible by , so the moment series is
Now introduce the standard noncrossing-partition transform , defined by 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 . Then and . Plugging into gives
Define Then , and the previous equation becomes So is the compositional inverse of .
That inverse series is the whole game.
For type , let Expanding chooses colors for the letters of type .
Now consider the multivariate generating series where 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 is just the sum of the type contributions, and the moment-cumulant relation gives
This is already a correct functional equation, but in this form it looks like something that wants to punch you in the face.
We only need the coefficient for fixed . So write Then the coefficient of in is exactly the desired answer.
Now set Then the equation becomes
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';
}
}