The problem requires counting unlabeled unrooted colored trees with maximum independent set size falling in a given range, which is an extended variant of classic unlabeled unrooted tree counting, so we can adapt standard techniques for that problem.
For unlabeled rooted trees, the generating function $F$ satisfies $F = x\mathcal{E}(F)$, where $\mathcal{E}(F)$ is the Euler transform of $F$. Euler transform can be expresssed as: $$\mathcal{E}(F) = \prod_{i=1}^{\infty} (1-x^i)^{-[x^i]F} = \exp\left(\sum_{i=1}^{\infty} \frac{F(x^i)}{i}\right)$$ Its combinatorial meaning corresponds to unbounded knapsack for combining subtrees. To compute Euler transform, we can explicitly maintain the sum term and use incremental online convolution to compute the exponential.
For unrooted trees, the standard counting approach counts each tree exactly once when the root is set to its centroid, and handles the special case of trees with two centroids separately.
Adding vertex coloring to the problem only introduces a constant multiplicative factor of $c$ (the number of colors) for the root node, so the base equation becomes $F = cx\mathcal{E}(F)$.
To incorporate the maximum independent set constraint, we use the property that the greedy algorithm (select a node if none of its children are selected, working bottom-up) correctly computes the maximum independent set for trees. For a rooted tree, we split into two generating functions: $F$ tracks trees where the root is selected in the maximum independent set, $G$ tracks trees where the root is not selected. We only track the size of the maximum independent set as the exponent, dropping the total node count dimension to keep complexity feasible.
The system of equations for rooted trees becomes: $$ \begin{cases} F = cx\mathcal{E}(G) \ G = c\mathcal{E}(F + G) - c\mathcal{E}(G) \end{cases} $$
When computing incrementally, note that $[x^i]G$ contributes to both $[x^i]\mathcal{E}(F+G)$ and $[x^i]\mathcal{E}(G)$, and the contributions cancel out, eliminating circular dependency when updating coefficients.
The core challenge of this problem is converting unrooted tree counting to rooted tree counting, since we do not track total tree size we cannot use the standard centroid definition directly. We need a unique anchor per unrooted tree that can be defined based on maximum independent set properties, with special handling for the small number of cases where multiple anchors exist.
A natural approach extends the centroid definition to a weighted centroid weighted by the maximum independent set size of subtrees. Define a black rooted tree as one where the root is selected in the greedy maximum independent set, and a white tree as one where the root is not selected. Multiple weighted centroids can only occur in two scenarios: two white trees with equal maximum independent set size connected directly by an edge, or two black trees of equal size connected through an intermediate white node. Adjusting the centroid definition to avoid three centroids still leaves an edge case: a tree of two nodes connected by an edge will have no node counted as a centroid, because the greedy maximum independent set changes when the root is changed, making the original definition ill-posed.
This issue only affects two cases: two white trees directly connected (two valid centroids) and two black trees directly connected (no valid centroid). We can leave the general definition as-is and correct the total count for these special cases separately with inclusion-exclusion.
A key insight simplifies the inclusion-exclusion: since we only need the sum of counts for maximum independent set size in a range, we do not need a full extra divide-and-conquer NTT pass to correct for non-centroid roots.
Brute-force Euler transform implementation:
#include <cstdio>
using namespace std;
const int MOD = 998244353;
const int MAX_N = 500003;
typedef long long ll;
int read(){ /* input implementation omitted */ }
int LOWER, n, c;
int sel[MAX_N], nsel[MAX_N], tsel[MAX_N], tnsel[MAX_N], ans_arr[MAX_N], inv[MAX_N];
inline void add(int &x, int v){ if((x += v) >= MOD) x -= MOD; }
inline void sub(int &x, int v){ if((x -= v) < 0) x += MOD; }
void update_transform(int *w, int step, int coeff) {
for(int t = n; ~t; --t) {
for(int i = 1, pos = step, cur = 1; t + pos <= n; ++i, pos += step) {
cur = (ll)cur * inv[i] % MOD * (coeff + i - 1) % MOD;
w[t + pos] = (w[t + pos] + (ll)cur * w[t]) % MOD;
}
}
}
int main(){
LOWER = read(); n = read(); c = read(); inv[1] = 1;
tsel[0] = tnsel[0] = 1;
for(int i = 2; i <= n; ++i) inv[i] = (ll)inv[MOD % i] * (MOD - MOD/i) % MOD;
for(int i = 1; i <= n; ++i) {
sel[i] = (ll)tsel[i-1] * c % MOD;
add(nsel[i], (ll)(tnsel[i] + sel[i]) * c % MOD);
sub(nsel[i], (ll)tsel[i] * c % MOD);
update_transform(tsel, i, nsel[i]);
update_transform(tnsel, i, (sel[i] + nsel[i]) % MOD);
}
for(int i = 1; i <= n; ++i) {
add(ans_arr[i] = sel[i], nsel[i]);
for(int j = (i+1)/2; j < i; ++j) sub(ans_arr[i], (ll)nsel[j] * sel[i - j] % MOD);
for(int j = (i+2)/2; j <= i; ++j) {
sub(ans_arr[i], (ll)(sel[j] + nsel[j]) * nsel[i - j] % MOD);
sub(ans_arr[i], (ll)sel[j] * sel[i - j + 1] % MOD);
}
}
for(int i = 1; 2*i <= n; ++i) if(nsel[i]) {
sub(ans_arr[2*i], ((ll)nsel[i] * (nsel[i] - 1) / 2) % MOD);
}
for(int i = 1; 2*i - 1 <= n; ++i) if(sel[i]) {
add(ans_arr[2*i - 1], ((ll)sel[i] * (sel[i] + 1) / 2) % MOD);
}
int total_ans = 0;
for(int i = LOWER; i <= n; ++i) add(total_ans, ans_arr[i]);
printf("%d\n", total_ans);
return 0;
}
Online convolution implementation:
#include <cstdio>
#include <algorithm>
#define ADD_ARR(arr) for(int i = mid+1; i <= r; ++i) add((arr)[i], C[i - l])
using namespace std;
const int MOD = 998244353;
const int MAX_N = 500003;
const int NTT_SIZE = 1 << 20;
typedef long long ll;
int read(){ /* input implementation omitted */ }
int qp(int a, int b = MOD-2){ /* fast exponentiation implementation omitted */ }
int LOWER, n, c, ntt_len, inv_len, bit_level;
int sel[MAX_N], nsel[MAX_N], af[MAX_N], ag[MAX_N], tsel[MAX_N], tnsel[MAX_N], inv[MAX_N], pre_sel[MAX_N], pre_nsel[MAX_N];
int rev[NTT_SIZE], cw[NTT_SIZE + 1];
inline void add(int &x, int v){ if((x += v) >= MOD) x -= MOD; }
inline void sub(int &x, int v){ if((x -= v) < 0) x += MOD; }
void NTT(int *E){ /* DFT implementation omitted */ }
void INTT(int *E){ /* inverse DFT implementation omitted */ }
void init_ntt(int cur_len){ /* NTT initialization omitted */ }
int A[NTT_SIZE], B[NTT_SIZE], C[NTT_SIZE];
void convolve(int *a_l, int *a_r, int *b_l, int *b_r) {
++a_r; ++b_r;
int ap = 0, bp = 0;
do A[ap++] = *a_l++; while(a_l != a_r);
do B[bp++] = *b_l++; while(b_l != b_r);
NTT(A); NTT(B);
for(int i = 0; i < ntt_len; ++i) {
C[i] = (ll)A[i] * B[i] % MOD;
A[i] = B[i] = 0;
}
INTT(C);
}
void solve(int l, int r) {
if(l == r) {
tsel[r] = (ll)(tsel[r] + af[r]) * inv[r] % MOD;
tnsel[r] = (ll)(tnsel[r] + ag[r]) * inv[r] % MOD;
sel[r] = (ll)tsel[r-1] * c % MOD;
nsel[r] = ((ll)sel[r] + tnsel[r] + MOD - tsel[r]) * c % MOD;
int vf = (ll)r * nsel[r] % MOD;
int vg = (ll)r * (sel[r] + nsel[r]) % MOD;
for(int t = r; t <= n; t += r) {
add(af[t], vf); add(ag[t], vg);
}
add(tsel[r], nsel[r]);
add(tnsel[r], sel[r]); add(tnsel[r], nsel[r]);
return;
}
int mid = (l + r) >> 1;
solve(l, mid);
init_ntt(r - l);
convolve(tsel + l, tsel + mid, af, af + (min(r - l, mid))); ADD_ARR(tsel);
convolve(af + l, af + mid, tsel, tsel + (min(r - l, l - 1))); ADD_ARR(tsel);
convolve(tnsel + l, tnsel + mid, ag, ag + (min(r - l, mid))); ADD_ARR(tnsel);
convolve(ag + l, ag + mid, tnsel, tnsel + (min(r - l, l - 1))); ADD_ARR(tnsel);
solve(mid + 1, r);
}
int main(){
LOWER = read(); n = read(); c = read(); inv[1] = 1;
tsel[0] = tnsel[0] = 1;
for(int i = 2; i <= n; ++i) inv[i] = (ll)inv[MOD % i] * (MOD - MOD/i) % MOD;
solve(1, n);
for(int i = 1; i <= n; ++i) {
pre_sel[i] = pre_sel[i-1]; add(pre_sel[i], sel[i]);
pre_nsel[i] = pre_nsel[i-1]; add(pre_nsel[i], nsel[i]);
}
int total_ans = 0;
for(int i = LOWER; i <= n; ++i) {
add(total_ans, sel[i]); add(total_ans, nsel[i]);
}
for(int i = (LOWER + 1)/2; 2*i <= n; ++i) if(nsel[i]) {
sub(total_ans, ((ll)nsel[i] * (nsel[i] - 1) / 2) % MOD);
}
for(int i = (LOWER + 2)/2; 2*i - 1 <= n; ++i) if(sel[i]) {
add(total_ans, ((ll)sel[i] * (sel[i] + 1) / 2) % MOD);
}
for(int i = 1; i <= n; ++i) {
int l = max(LOWER - i, 1);
int r = min(n - i, i);
if(l <= r) {
total_ans = (total_ans + (ll)(pre_sel[l-1] - pre_sel[r] + MOD) * nsel[i]) % MOD;
}
l = max(LOWER - i, 0);
r = min(n - i, i - 1);
if(l <= r) {
total_ans = (total_ans + (ll)(pre_nsel[l-1] - pre_nsel[r] + MOD) * (sel[i] + nsel[i])) % MOD;
total_ans = (total_ans + (ll)(pre_sel[l] - pre_sel[r+1] + MOD) * sel[i]) % MOD;
}
}
printf("%d\n", total_ans);
return 0;
}