Let's focus on the number of occurrences of each element: let $$$c_i$$$ be the number of time that element $$$i$$$ occurs in the array. To keep a subsequence is to decrease each $$$c_i$$$ by a non-negative integer. We need to make all the non-zero elements in $$$c$$$ the same while maximizing the sum of $$$c$$$.

For each integer $$$x$$$: suppose the value of all non-zero elements in $$$c$$$ is $$$x$$$, what is the length of the longest subsequence?

By Caylex.

def solve():
    n = int(input())
    ls = list(map(int, input().split()))
    a = [0] * n
    for x in ls:
        a[x - 1] += 1
    a = sorted(a)
    s = sum(a)
    ans = s
    for i in range(n):
        ans = min(ans, s - (n - i) * a[i])
    print(n - ans)
    return
 
 
t = int(input())
for _ in range(t):
    solve()

How to check if there is at least one way to choose the sets?

How to greedily find two more solutions to choose the sets?

By Caylex.

def solve():
    n, m = map(int, input().split())
    cnt = [0] * m
    t = 0
    v = []
    for i in range(n):
        v.append([])
        inp = list(map(int, input().split()))
        l = inp[0]
        for j in range(1, l + 1):
            x = inp[j]
            x -= 1
            t += not cnt[x]
            cnt[x] += 1
            v[i].append(x)
    ans = t == m
    for i in range(n):
        for x in v[i]:
            cnt[x] -= 1
            t -= not cnt[x]
        ans += t == m
        for x in v[i]:
            t += not cnt[x]
            cnt[x] += 1
    print("YES" if ans >= 3 else "NO")
    return


t = int(input())
for _ in range(t):
    solve()

How to check if an element $$$x$$$ is stable? For an element $$$p_i$$$, we only care about whether $$$p_i \lt x$$$, $$$p_i=x$$$ or $$$p_i \gt x$$$.

If an element is not stable, we can choose $$$m$$$ correctly in the first step to exclude $$$x$$$ from interval $$$[l,r]$$$.

By Caylex.

def solve():
    n = int(input())
    s = str(input())
    l = 0
    p = []
    for i in range(n):
        p.append(i + 1)
    while l < n:
        if s[l] == "1":
            l += 1
            continue
        r = l
        while r + 1 < n and s[r + 1] == "0":
            r += 1
        if r - l + 1 == 1:
            print("NO")
            return
        ok = 1
        for i in range(l, r):
            p[i] = i + 2
        p[r] = l + 1
        l = r + 1
    print("YES")
    for i in range(n):
        print(p[i], end=" ")
    print()
    return


t = int(input())
for _ in range(t):
    solve()

Try a few cases by hand or guess by the samples, what's the maximized answer for $$$l=0,r=n$$$?

We need to make the bitwise AND of $$$a_i$$$ and $$$b_i$$$ all equal to $$$0$$$. Let's remind you that the bitwise AND of $$$x$$$ and $$$2^k-1-x$$$ is $$$0$$$ for $$$0\le x \lt 2^k$$$. We can try to pair integers $$$x$$$ and $$$y$$$ by letting $$$a_x=y,a_y=x$$$.

Construct from the back.

By Caylex.

def solve():
    l, r = map(int, input().split())
    pw = 1
    ans = 0
    while pw < r:
        pw = pw << 1 | 1
    s = set()
    for i in range(l, r + 1):
        s.add(i)
    a = [0] * (r + 1)
    for i in range(r, l - 1, -1):
        while pw - i not in s:
            pw >>= 1
        a[i] = pw - i
        ans += i | a[i]
        s.remove(pw - i)
    print(ans)
    print(*a[l : r + 1])
    return


t = int(input())
for _ in range(t):
    solve()

Let's consider from the highest bits to the lowest bits. If the highest bits of $$$l$$$ and $$$r$$$ are the same, we can just ignore this bit and go on to a lower bit. We can still try to pair integers $$$x$$$ and $$$y$$$ by letting $$$a_x=y,a_y=x$$$.

Consider the highest bit that $$$l$$$ and $$$r$$$ differs. There exist an integer $$$t$$$ that this bit of $$$l\ldots t-1$$$ are $$$0$$$ and this bit of $$$t\ldots r$$$ are $$$1$$$. We need to pair a $$$0$$$ with a $$$1$$$ as more as possible while keeping the bitwise AND of lower bits small. Can we turn it into a smaller subproblem?

In D1, it is optimal to pair $$$x$$$ with $$$2^k-1-x$$$. Here, for interval $$$[l,r]$$$ and integer $$$t$$$ defined in Hint2, it is optimal to pair $$$x$$$ with $$$2t-1-x$$$.

By Caylex.

def g(x, j):
    return x >> j & 1


def solve():
    L, R = map(int, input().split())
    a = [0] * (R - L + 1)

    def solve(l, r, j):
        if l > r:
            return
        if l == r:
            a[l - L] = r
            return
        mid = l
        while mid + 1 <= r and g(mid + 1, j) == g(l, j):
            mid += 1
        if mid == r:
            solve(l, r, j - 1)
            return
        tl = mid + 1
        tr = mid
        while tl - 1 >= l and tr + 1 <= r:
            tl -= 1
            tr += 1
            a[tl - L] = tr
            a[tr - L] = tl
        solve(l, tl - 1, j - 1)
        solve(tr + 1, r, j - 1)
        return

    solve(L, R, 29)
    ans = 0
    for i in range(L, R + 1):
        ans += a[i - L] | i
    print(ans)
    print(*a)
    return


t = int(input())
for _ in range(t):
    solve()

It turns out that we can solve it greedily from the lowest bit to the highest bit. For each element in $$$[l,r]$$$, let's reverse its binary representation and insert it into a 01-trie. Then for each index $$$i$$$, we can do dfs on the trie while trying to keep every bit different with $$$i$$$. Formally, each step we try to go to the opposite side of each bit of $$$i$$$, and go on the other side if that subtree is empty.

We can prove that this gives the same construction with the other method constructing from the higher bits. You can read the implementation for better understanding.

#include<bits/stdc++.h>
using namespace std;
#define int         long long
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
int l,r;
int tot;
int son[6000005][2],cnt[6000005];
void Main() {
	cin>>l>>r;
	tot=1;
	son[0][0]=son[0][1]=-1;
	REP(i,l,r+1){
		int x=0;
		REP(j,0,30){
			int &p=son[x][(i>>j)&1];
			if(p==-1)p=tot++,cnt[p]=0,son[p][0]=son[p][1]=-1;
			x=p;++cnt[x];
		}
	}
	int ans=0;
	vector<int>res;
	REP(i,l,r+1){
		int x=0,t=0;
		REP(j,0,30){
			int o=!((i>>j)&1);
			if(son[x][o]==-1||!cnt[son[x][o]])o=!o;
			t|=o<<j;
			x=son[x][o];--cnt[x];
		}
		res.pb(t);ans+=i|t;
	}
	cout<<ans<<'\n';
	for(auto i:res)cout<<i<<' ';
	cout<<'\n';
}
void TC() {
    int tc=1;
    cin>>tc;
	while(tc--){
		Main();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}
/*
1. CLEAR the arrays (ESPECIALLY multitests)
2. DELETE useless output
 */

For an array $$$b$$$ and an integer $$$x$$$ that does not appear in $$$b$$$, let $$$g(b,x)$$$ be the number of elements in $$$b$$$ that is greater than $$$x$$$. Then the weight of $$$b$$$ equals to $$$g(b,\rm{mex}(b))$$$. However, when $$$x\neq \rm{mex}(b)$$$, what will happen to $$$g(b,x)$$$?

For a fixed right endpoint $$$r$$$, let's enumerate over all the $$$x$$$ and maximize $$$g(a[l\dots r],x)$$$ while keeping $$$x$$$ no occurred in $$$a[l\dots r]$$$. It's clear that the value is maximized for $$$x$$$ if $$$l$$$ equals the last occurrence of $$$x$$$ plus $$$1$$$.

Use a segment tree to maintain the answer for each $$$x$$$.

By LMydd0225.

#include<bits/stdc++.h>
using namespace std;
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
int n;
int a[300005];
struct ds{
	int seg[1200005],tag[1200005];
	void build(int l,int r,int p){
		seg[p]=tag[p]=0;
		if(l==r)return;
		int m=(l+r)>>1;
		build(l,m,p*2+1);build(m+1,r,p*2+2);
	}
	void update(int l,int r,int s,int t,int p){
		if(l<=s&&t<=r){
			++seg[p];++tag[p];
			return;
		}
		int m=(s+t)>>1;
		if(m>=l)update(l,r,s,m,p*2+1);
		if(m<r)update(l,r,m+1,t,p*2+2);
		seg[p]=max(seg[p*2+1],seg[p*2+2])+tag[p];
	}
	void point(int pos,int l,int r,int p,int val=0){
		if(l==r)return seg[p]=-val,void();
		int m=(l+r)>>1;
		val+=tag[p];
		if(m>=pos)point(pos,l,m,p*2+1,val);
		else point(pos,m+1,r,p*2+2,val);
		seg[p]=max(seg[p*2+1],seg[p*2+2])+tag[p];
	}
}seg;
void Main() {
	cin>>n;
	REP(i,0,n)cin>>a[i];
	seg.build(0,n,0);
	REP(i,0,n){
		seg.update(0,a[i],0,n,0);
		seg.point(a[i],0,n,0);
		cout<<seg.seg[0]<<' ';
	}
	cout<<'\n';
}
void TC() {
    int tc=1;
    cin>>tc;
	while(tc--){
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}
/*
1. CLEAR the arrays (ESPECIALLY multitests)
2. DELETE useless output
 */

There's a harder version of the problem. Can you solve it?

There're $$$q$$$ queries: $$$x,y$$$. Find $$$\max_{x\le l\le r\le y}w(l,r)$$$. $$$n,q\le 3\cdot 10^5$$$.

Let's focus on bubble sort first. For an element $$$p_x$$$, how many rounds does it take to make all the numbers before the position small than $$$p_x$$$?

Try to count something else instead of permutations $$$p$$$. Can you come up with an $$$O(n^2)$$$ dynamic programming method?

The number of restrictions is small ($$$m\le 1000$$$). Can we decrease the number of states from $$$O(n^2)$$$ to $$$O(nm)$$$ and $$$O(m^2)$$$? How to calculate the coefficients then?

By LMydd0225.

#include<bits/stdc++.h>
using namespace std;
#define MOD         998244353
#define speMOD      2933256077ll
#define int         long long
#define pii         pair<int,int>
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
#define over(x)     {cout<<(x)<<endl;return;}
#define lowbit(x)   ((x)&(-(x)))
#define cntbit(x)   __builtin_popcount(x)
#define rf(v)       shuffle(all(v),sd);
#define deal(v)     sort(all(v));v.erase(unique(v.begin(),v.end()),v.end())
#define lbound(v,x) lower_bound(all(v),x)-v.begin()
int qpow(int a,int b,int m=MOD,int res=1){
	a%=m;
	while(b>0)res=(b&1)?(res*a%m):(res),a=a*a%m,b>>=1;
	return res;
}
int fac[2000005],inv[2000005];
void init(int n){
	fac[0]=inv[0]=1;
	REP(i,1,n+1)fac[i]=fac[i-1]*i%MOD;
	inv[n]=qpow(fac[n],MOD-2);
	for(int i=n-1;i>=1;--i)inv[i]=inv[i+1]*(i+1)%MOD;
}
int binom(int x,int y){
	return x<y||y<0? 0:fac[x]*inv[y]%MOD*inv[x-y]%MOD;
}
struct restr{
	int l,r,k;
};
int n,m;
int r1[1000005],r2[1000005];
int dp[10005],sum[10005];
int qry1(int l,int r,int k){
	if(k<l)return qpow(k+1,r-l+1);
	if(k>r)return fac[r+1]*inv[l]%MOD;
	return fac[k+1]*inv[l]%MOD*qpow(k+1,r-k)%MOD;
}
int qry2(int l,int r,int k,int x){
	return (qry1(l,r,k)+MOD-qry1(l,r,x))%MOD;
}
void Main() {
	cin>>n>>m;
	vector<int>v;vector<restr>a;
	vector<int>t;
	REP(i,0,m){
		int k,l,r;
		cin>>k>>l>>r;
		--l;
		if(l>=0)v.pb(l);
		if(r<n)v.pb(r);
		t.pb(k);
		a.pb({l,r,k});
	}
	v.pb(n-1);t.pb(-1);t.pb(n-1);
	deal(v);deal(t);
	for(auto i:v)r1[i]=t.size()-1,r2[i]=0;
	for(auto i:a){
		if(i.k==n-1&&i.r<n-1)over(0)
		if(i.l>=0)r1[i.l]=min(r1[i.l],lbound(t,i.k));
		if(i.r<n)r2[i.r]=max(r2[i.r],lbound(t,i.k));
	}
	int lst=0,m=t.size();dp[1]=1;REP(i,2,m+1)dp[i]=0;
	for(auto i:v){
		REP(j,1,m)sum[j]=dp[j];
		REP(j,2,m)(sum[j]+=sum[j-1])%=MOD;
		REP(j,1,m)if(j>r2[i]&&j<=r1[i]){
			dp[j]=(dp[j]*qry1(lst,i,t[j])+sum[j-1]*qry2(lst,i,t[j],t[j-1]))%MOD;
		}else dp[j]=0;
		lst=i+1;
	}
	int res=0;
	REP(j,1,m+1)(res+=dp[j])%=MOD;
	cout<<res<<endl;
}
void TC() {
	init(1000000);
    int tc=1;
    cin>>tc;
	while(tc--){
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}
/*
1. CLEAR the arrays (ESPECIALLY multitests)
2. DELETE useless output
 */

How many odd integers are there in a pair $$$(a, b, c)$$$?

It's always better to put exactly $$$2$$$ odd integers and exactly $$$1$$$ even integer in one pair of $$$(a, b, c)$$$. Is there a way to always choose optimal $$$(a, b, c)$$$?

Try consecutive integers.

T = int(input())
 
for _ in range(T) :
	l, r = map(int, input().split())
	print(((r + 1) // 2 - l // 2) // 2)

Consider the position of the maximum value in one array. If there's multiple maximum values, choose any. What can you observe?

T = int(input())
 
for _ in range(T) :
	n, m = map(int, input().split())
	v = max(map(int, input().split()))
	for __ in range(m) :
		c, l, r = input().split()
		l = int(l)
		r = int(r)
		if (l <= v <= r) :
			if (c == '+') :
				v = v + 1
			else :
				v = v - 1
		print(v, end = ' ' if __ != m - 1 else '\n')

If you can perform an operation that decreases $$$a$$$ or $$$b$$$ from an element, will it be different?

If $$$a = b$$$, can you solve the problem?

import math
 
T = int(input())
for _ in range(T) :
	n, a, b = map(int, input().split())
	d = math.gcd(a, b)
	a = list(map(int, input().split()))
	a = [x % d for x in a]
	a.sort()
	res = a[n - 1] - a[0]
	for i in range(1, n) :
		res = min(res, a[i - 1] + d - a[i])
	print(res)

How to calulate the weight of a leaf quickly?

Consider a formed string. Let's delete useless part that won't contribute to the count of $$$\tt{01}$$$ and $$$\tt{10}$$$ substrings. What will be the rest?

The root is somehow important. When to fill $$$s_1$$$?

T = int(input())
 
for _ in range(T) :
    n = int(input())
    x, y, z, w = 0, 0, 0, 0
    deg = [0] * (n + 1)
    for __ in range(n - 1) :
        u, v = map(int, input().split())
        deg[u] += 1
        deg[v] += 1
    s = " " + input()
    for i in range(2, n + 1) :
        if deg[i] == 1 :
            if s[i] == '?' :
                z += 1
            elif s[i] == '0' :
                x += 1
            else :
                y += 1
        elif s[i] == '?' :
            w += 1
    if s[1] == '0' :
        print(y + (z + 1) // 2)
    elif s[1] == '1' :
        print(x + (z + 1) // 2)
    else :
        print(max(x, y) + (z + (w % 2 if x == y else 0)) // 2)

The previous statement of the problem did not directly tell you that the vertices are numbered by dfs order. It was changed to make the statement more readable. Can you prove that the following property is the same as telling you the vertices are numbered by dfs order?

On a tree with nodes numbered with dfs order, we have to find the distance between two consecutive nodes. What can you think of?

How to calculate the maximum value of $$$\text{dist}(i,(i\bmod n)+1)$$$? Is it different when the weight of all of its edges are known?

#include <bits/stdc++.h>
 
namespace FastIO {
	template <typename T> inline T read() { T x = 0, w = 0; char ch = getchar(); while (ch < '0' || ch > '9') w |= (ch == '-'), ch = getchar(); while ('0' <= ch && ch <= '9') x = x * 10 + (ch ^ '0'), ch = getchar(); return w ? -x : x; }
	template <typename T> inline void write(T x) { if (!x) return; write<T>(x / 10), putchar(x % 10 ^ '0'); }
	template <typename T> inline void print(T x) { if (x < 0) putchar('-'), x = -x; else if (x == 0) putchar('0'); write<T>(x); }
	template <typename T> inline void print(T x, char en) { if (x < 0) putchar('-'), x = -x; else if (x == 0) putchar('0'); write<T>(x), putchar(en); }
}; using namespace FastIO;
 
#define MAXN 200001
int fa[MAXN], dep[MAXN];
int c1[MAXN], c2[MAXN], len[MAXN];
 
void solve() {
	int N = read<int>(); long long w = read<long long>();
	for (int i = 2; i <= N; ++i) fa[i] = read<int>();
	for (int i = 2; i <= N; ++i) dep[i] = dep[fa[i]] + 1;
	for (int i = 1; i <= N; ++i) len[i] = c1[i] = 0;
	for (int i = 1, x, y; i <= N; ++i) {
		x = i, y = (i == N ? 1 : i + 1);
		while (x != y) {
			if (dep[x] < dep[y]) std::swap(x, y);
			(c1[x] ? c2[x] : c1[x]) = i, x = fa[x], ++len[i];
		}
	}
	long long sum = 0, sur = N;
	for (int i = 1, x; i < N; ++i) {
		x = read<int>(), sum += read<long long>();
		if ((--len[c1[x]]) == 0) --sur;
		if ((--len[c2[x]]) == 0) --sur;
		print<long long>(sum * 2 + sur * (w - sum), " \n"[i == N - 1]);
	}
}
 
int main() { int T = read<int>(); while (T--) solve(); return 0; }

Can you solve the problem if the nodes are not numbered by dfs order?

Let's try to judge whether an array is brilliant, and try to sort it first. What happens if we minus $$$1$$$ from each element?

We only care about the relative positions between the elements. Let's consider the difference array $$$b$$$ of array $$$a$$$: $$$b_i=a_i-a_{i-1}(i\ge 2)$$$.

What will the difference array look like after no more operations can be done?

#include<bits/stdc++.h>
using namespace std;
#define int         long long
#define pii         pair<int,int>
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
#define over(x)     {cout<<x<<endl;return;}
#define cntbit(x)   __buitin_popcount(x)
int n;
int a[400005],lsteq[400005];
struct two_pointers{
	int l,r,mid;
	int val[400005],rt;
	void init(){
		l=r=mid=rt=0;
		val[0]=a[0];
	}
	inline void addr(){++r;rt=__gcd(rt,a[r]);}
	inline void addl(){if(l+1>r)addr();++l;if(l>mid){rt=0,mid=r;val[mid]=a[mid];
	for(int i=mid-1;i>=l;--i)val[i]=__gcd(val[i+1],a[i]);}}
	inline int query(){return __gcd(val[l],rt);}
};
void Main() {
	cin>>n;
	REP(i,0,n)cin>>a[i];
	REP(i,1,n)a[n-i]-=a[n-i-1];
	REP(i,1,n)a[i-1]=llabs(a[i]);
	--n;
	lsteq[n]=n;
	for(int i=n-1;i>=0;--i)lsteq[i]=(a[i]? i:lsteq[i+1]);
	two_pointers p;p.init();
	int ans=1,cnt=1;
	REP(i,0,n)if(lsteq[i]==i){
		while(p.l<i)p.addl();
		while(p.r<n&&cntbit(p.query())>1)p.addr();
		ans+=(n-p.r)*cnt+1;cnt=1;
	}else ans+=lsteq[i]-i+1,++cnt;
	cout<<ans<<endl;
}
void TC() {
	int tc=1;
	cin>>tc;
	while(tc--){
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}

#include<bits/stdc++.h>
using namespace std;
#define int         long long
#define pii         pair<int,int>
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
#define over(x)     {cout<<x<<endl;return;}
#define cntbit(x)   __buitin_popcount(x)
int n;
int a[400005],res[400005],lsteq[400005];
int st[400005][20];
int query(int l,int r){
	int s=__lg(r-l+1);
	return __gcd(st[l][s],st[r-(1<<s)+1][s]);
}
void Main() {
	cin>>n;
	REP(i,0,n)cin>>a[i];
	--n;
	if(!n)over(1)
	REP(i,0,n)st[i][0]=llabs(a[i]-a[i+1]);
	REP(j,0,__lg(n)){
		REP(i,0,n-(1<<(j+1))+1){
			st[i][j+1]=__gcd(st[i][j],st[i+(1<<j)][j]);
		}
	}
	lsteq[n]=n;
	for(int i=n-1;i>=0;--i)lsteq[i]=(a[i]==a[i+1]? lsteq[i+1]:i);
	int ans=1;
	int l=0,r=0;
	REP(i,0,n){
		l=max(l,lsteq[i]);r=max(r,l);
		while(r<n&&cntbit(query(l,r))>1)++r;
		ans+=n-r;ans+=lsteq[i]-i+1;
	}
	cout<<ans<<endl;
}
void TC() {
	int tc=1;
	cin>>tc;
	while(tc--){
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}

There's a harder but similar version of the problem. Can you solve it?

There're $$$q$$$ queries: $$$L,R,k$$$. Find the number of pairs $$$(l,r)$$$ that $$$L\le l\le r\le R$$$, and $$$(\lvert Q({a_l\dots a_r}) \rvert-1) \cdot k\ge \max(a_l\dots a_r) - \min(a_l\dots a_r) $$$, where $$$Q(array)$$$ means the result of putting the array into the set.

Let's use another way to describe the problem. Let's define an array good, that it can be reordered to satisfy the condition. How many numbers do we have to change to make it good?

To handle this, we need to know what is the optimal way to reorder an array.

Sort the array. Suppose $$$a_1$$$ is the minimum and $$$a_n$$$ is the maximum. Where to place the two elements? Can you prove this is optimal?

How to maintain the answer better using the properties?

#include<bits/stdc++.h>
using namespace std;
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
#define over(x)     {cout<<(x)<<endl;return;}
int n,q,r;
int a[300005];
int sum[300005];
int ql[300005],qr[300005];
int ans[300005];
int cnt[300005],cnt2[300005];
void Main() {
	cin>>n>>q>>r;
	REP(i,0,n)cin>>a[i];
	sum[0]=0;
	int B=sqrt(r);
	REP(i,0,n)sum[i+1]=sum[i]+(a[i]<=B);
	REP(i,0,q){
		int x,y;
		cin>>x>>y;
		--x,--y;
		ql[i]=x;qr[i]=y;
		ans[i]=0;
		if(sum[y+1]-sum[x]<(y-x+1)/2)ans[i]=(y-x+1)/2-sum[y+1]+sum[x];
	}
	REP(i,0,B){
		REP(j,0,n+1)cnt[j]=cnt2[j]=0;
		REP(j,0,n){
			cnt[j+1]=cnt[j]+(a[j]<=i);
			cnt2[j+1]=cnt2[j]+(a[j]<=r/(i+1));
		}
		REP(j,0,q){
			int x=ql[j],y=qr[j],l=y-x+1;
			int c1=cnt2[y+1]-cnt2[x],c2=cnt[y+1]-cnt[x];
			ans[j]=max(ans[j],min((l-c1-c2+1)/2,l/2-c2));
		}
	}
	REP(i,0,q)cout<<ans[i]<<' ';
	cout<<endl;
}
void TC() {
    int tc=1;
    cin>>tc;
	while(tc--){
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}
/*
1. CLEAR the arrays (ESPECIALLY multitests)
2. DELETE useless output
 */

Consider the degress of vertices in a binary tree. What can you observe?

Choose a root. It should have degree $$$\le 2$$$. Is it always possible to build a binary tree on the initial tree if all nodes have degree $$$\le 3$$$? What will be the depth of the binary tree with the given root?

Consider the binary depth of the full tree. We should to go over every possible root, but how to find the maximum distance from the root to other nodes? Are there any special nodes?

Diameter.

Is it possible to maintain the answer for all the roots when adding a leaf to the tree?

#include<bits/stdc++.h>
using namespace std;
#define int         long long
#define pii         pair<int,int>
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(int)(e);++i)
#define over(x)     {cout<<(x)<<endl;return;}
struct ds{
	int seg[2000005],tag[2000005];
	void build(int l,int r,int p){
		seg[p]=(l==0? 0:1e18),tag[p]=0;
		if(l==r)return;
		int m=(l+r)>>1;
		build(l,m,p*2+1);build(m+1,r,p*2+2);
	}
	void pushdown(int p){
		if(!tag[p])return;
		tag[p*2+1]+=tag[p];tag[p*2+2]+=tag[p];
		seg[p*2+1]+=tag[p];seg[p*2+2]+=tag[p];
		tag[p]=0;
	}
	void add(int l,int r,int s,int t,int p,int val){
		if(l<=s&&t<=r){
			seg[p]+=val;tag[p]+=val;
			return;
		}
		int m=(s+t)>>1;pushdown(p);
		if(m>=l)add(l,r,s,m,p*2+1,val);
		if(m<r)add(l,r,m+1,t,p*2+2,val);
		seg[p]=min(seg[p*2+1],seg[p*2+2]);
	}
	void update(int pos,int l,int r,int p,int val){
		if(l==r){
			seg[p]=val;
			return;
		}
		int m=(l+r)>>1;pushdown(p);
		if(m>=pos)update(pos,l,m,p*2+1,val);
		else update(pos,m+1,r,p*2+2,val);
		seg[p]=min(seg[p*2+1],seg[p*2+2]);
	}
}seg;
int n;
vector<int>v[500005];
int fa[500005],an[500005][21];
int dep[500005],dfn[500005],rev[500005];
int tot,deg[500005],ls[500005];
void dfs(int x,int pre,int d){
	dep[x]=d;fa[x]=pre;an[x][0]=fa[x];ls[x]=1;
	REP(i,0,__lg(n+1))if(an[x][i]==-1)an[x][i+1]=-1;else an[x][i+1]=an[an[x][i]][i];
	rev[tot]=x;dfn[x]=tot++;
	for(auto i:v[x])dfs(i,x,d+1),ls[x]+=ls[i];
}
int getlca(int x,int y){
	if(dep[x]<dep[y])swap(x,y);
	int d=dep[x]-dep[y];
	for(int i=__lg(d);i>=0;--i)if((d>>i)&1)x=an[x][i];
	if(x==y)return x;
	d=dep[x];
	for(int i=__lg(d);i>=0;--i)if((1<<i)<=dep[x]&&an[x][i]!=an[y][i])x=an[x][i],y=an[y][i];
	return fa[x];
}
int getan(int x,int y){
	if(dfn[y]>=dfn[x]&&dfn[y]<=ls[x]){
		int d=dep[y]-dep[x]-1;
		for(int i=0;(1<<i)<=d;++i)if((d>>i)&1)y=an[y][i];
		return y;
	}else return fa[x];
}
void updside(int x,int y){
	if(dfn[y]>=dfn[x]&&dfn[y]<=ls[x]){
		seg.add(0,n-1,0,n-1,0,1);
		x=getan(x,y);
		seg.add(dfn[x],ls[x],0,n-1,0,-1);
	}else seg.add(dfn[x],ls[x],0,n-1,0,1);
}
int getdist(int x,int y){return dep[x]+dep[y]-2*dep[getlca(x,y)];}
struct diameter{
	int x,y,len;
	int update(int z){
		if(x==y){
			int d=getdist(x,z);
			if(d<=len)return d+len;
			updside(y,z);
			y=getan(y,z);
			return d+len;
		}else{
			int d1=getdist(x,z),d2=getdist(y,z),d3=min(d1,d2);
			if(d3<=len)return d3+len+1;
			if(d1>d2)swap(x,y);
			updside(y,z);
			y=x;++len;return len*2;
		}
	}
	int query(int z){
		if(x==y)return len+getdist(x,z);
		else return len+min(getdist(x,z),getdist(y,z))+1;
	}
};
void Main() {
	cin>>n;
	REP(i,0,n)v[i].clear();
	REP(i,1,n){
		cin>>fa[i];--fa[i];
		v[fa[i]].pb(i);
	}
	tot=0;dfs(0,-1,0);seg.build(0,n-1,0);
	REP(i,0,n)deg[i]=0,ls[i]+=dfn[i]-1;
	diameter d={0,0,0};
	cout<<"1 ";
	REP(i,1,n){
		++deg[fa[i]];++deg[i];
		if(deg[fa[i]]==4){
			REP(j,i,n)cout<<-1<<' ';
			cout<<endl;return;
		}
		if(deg[fa[i]]==3)seg.update(dfn[fa[i]],0,n-1,0,1e18);
		seg.update(dfn[i],0,n-1,0,d.update(i));
		cout<<seg.seg[0]+1<<' ';
	}
	cout<<endl;
}
void TC() {
    int tc=1;
    cin>>tc;
	while(tc--){
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}
/*
1. CLEAR the arrays (ESPECIALLY multitests)
2. DELETE useless output
 */

Each column or row should be painted at most once. Suppose we paint each column and row exactly once. Then what does the number on one position tell us?

We can build a directed bipartite graph. How to judge whether it is possible to paint the graph? How to calculate the beauty of one single grid?

How to solve the original problem when the initial grid has positive beauty? How to optimize? Consider adjacent layers.

How to solve the original problem when the initial grid is impossible to be painted out? Are there too many ways to flip a block?

How to optimize the process of building the graph and doing topological sort?

#include<bits/stdc++.h>
using namespace std;
#define MOD         998244353
#define speMOD      2933256077ll
#define int         long long
#define pii         pair<int,int>
#define all(v)      v.begin(),v.end()
#define pb          push_back
#define REP(i,b,e)  for(int i=(b);i<(e);++i)
#define over(x)     {cout<<x<<endl;return;}
#define lowbit(x)   ((x)&(-(x)))
#define cntbit(x)   __builtin_popcount(x)
#define rf(v)       shuffle(all(v),sd);
mt19937 sd(std::random_device{}());
int qpow (int a,int b,int m=MOD,int res=1){
	a%=m;
	while(b>0)res=(b&1)?(res*a%m):(res),a=a*a%m,b>>=1;
	return res;
}
int fac[2000005],inv[2000005],si[2000005];
void init(int n){
	fac[0]=inv[0]=si[0]=1;
	REP(i,1,n+1)fac[i]=fac[i-1]*i%MOD;
	inv[n]=qpow(fac[n],MOD-2);
	for(int i=n-1;i>=1;--i)inv[i]=inv[i+1]*(i+1)%MOD;
	REP(i,1,n+1)si[i]=inv[i]*fac[i-1]%MOD;
}
int binom(int x,int y){
	return x<y||y<0? 0:fac[x]*inv[y]%MOD*inv[x-y]%MOD;
}
namespace fastIO {
#define BUF_SIZE 6000010
    bool IOerror = 0;
    inline char nc() {
        static char buf[BUF_SIZE], *p1 = buf + BUF_SIZE, *pend = buf + BUF_SIZE;
        if(p1 == pend) {
            p1 = buf;
            pend = buf + fread(buf, 1, BUF_SIZE, stdin);
            if(pend == p1) {
                IOerror = 1;
                return -1;
            }
        }
        return *p1++;
    }
    inline bool blank(char ch) {
        return ch == ' ' || ch == '\n' || ch == '\r' || ch == '\t';
    }
    inline void read(int &x) {
        char ch;
        while(blank(ch = nc()));
        if(IOerror)
            return;
        for(x = ch - '0'; (ch = nc()) >= '0' && ch <= '9'; x = x * 10 + ch - '0');
    }
#undef BUF_SIZE
};
using namespace fastIO;
int n,m;
int a[200005],b[200005];
vector<pii>c;
int d1[200005],d2[200005];
int visa[200005],visb[200005];
vector<int>h[200005],nodes[200005],buc[200005];
int X,Y;
int id1[200005],id2[200005],cnt[400005];
int solve(){
	REP(i,0,n+1)buc[i].clear();
	REP(i,0,n)buc[a[i]].pb(i);
	int num=0;
	REP(i,0,n+1)for(auto j:buc[i])id1[num++]=j;
	REP(i,0,n+1)buc[i].clear();
	REP(i,0,n)buc[b[i]].pb(i);
	num=0;
	REP(i,0,n+1)for(auto j:buc[i])id2[num++]=j;
	REP(i,0,n)d1[i]=d2[i]=-1;
	int ans=1,x=0,y=0,cur=0;
	while(x<n||y<n){
		if(x<n&&a[id1[x]]<=y){
			int sum=0;
			while(x<n&&a[id1[x]]<=y)++sum,d1[id1[x++]]=cur;
			if(cur)(ans*=fac[sum])%=MOD;
		}else if(y<n&&b[id2[y]]<=x){
			int sum=0;
			while(y<n&&b[id2[y]]<=x)++sum,d2[id2[y++]]=cur;
			if(cur)(ans*=fac[sum])%=MOD;
		}else return X=x,Y=y,0;
		++cur;
	}
	return ans;
}
int updater(int x,int y){
	int f=0,ret;
	for(auto i:c)if(i==make_pair(x,y))f=1;
	if(f)--b[y],++a[x],ret=solve(),++b[y],--a[x];
	else ++b[y],--a[x],ret=solve(),--b[y],++a[x];
	return ret;
}
void Main() {
	fastIO::read(n);fastIO::read(m);
	c.clear();
	REP(i,0,n)a[i]=b[i]=0;
	REP(i,0,n*2+2)cnt[i]=0;
	REP(i,0,n)visa[i]=visb[i]=0,h[i].clear(),nodes[i].clear();
	REP(i,0,m){
		int x,y;
		fastIO::read(x);fastIO::read(y);
		--x,--y;
		c.pb({x,y});
		++b[y];--a[x];
	}
	REP(i,0,n)a[i]+=n;
	int ans=solve();
	if(!ans){
		if(X==n-1||Y==n-1)over(0)
		int x=id1[X],y=id2[Y],f=0;
		for(auto i:c)if(i==make_pair(x,y)){f=1;break;}
		if(f){
			for(auto i:c)if(i.first==x)visb[i.second]=1;
			for(auto i:c)if(i.second==y)visa[i.first]=1;
			X=Y=-1;
			for(auto i:c)if(i.first!=x&&i.second!=y){
				if(!visa[i.first]&&!visb[i.second]){X=i.first,Y=i.second;break;}
			}
			if(X==-1)over(0);
		}else{
			X=Y=-1;
			for(auto i:c)h[i.first].pb(i.second);
			int sum=0;
			for(auto i:h[x])visa[i]=1,++sum;
			REP(i,0,n){
				int found=0;
				for(auto j:h[i])if(j==y){found=1;break;}
				if(!found)continue;
				vector<int>z;
				for(auto j:h[i])if(visa[j])visa[j]=0,--sum,z.pb(j);
				if(sum){
					X=i;
					REP(j,0,n)if(visa[j])Y=j;
					break;
				}
				for(auto j:z)++sum,visa[j]=1;
			}
			if(X==-1)over(0);
		}
		int x2=X,y2=Y;
		int ans=updater(x,y)+updater(x,y2)+updater(x2,y)+updater(x2,y2);
		cout<<(ans%MOD*qpow(n*n%MOD,MOD-2))%MOD<<endl;
		return;
	}
	REP(i,0,n)++cnt[d1[i]],++cnt[d2[i]];
	int res=0;
	(res+=ans*cnt[0]%MOD*(cnt[1]==1? cnt[2]+1:1))%=MOD;
	if(cnt[2]&&cnt[1])(res+=ans*(cnt[2]==1? cnt[3]+1:1))%=MOD;
	REP(i,2,2*n-1)if(cnt[i]&&cnt[i+1]){
		(res+=ans*(cnt[i+1]==1? cnt[i+2]+1:1)%MOD*(cnt[i]==1? cnt[i-1]+1:1))%=MOD;
	}
	cout<<res*qpow(n*n%MOD,MOD-2)%MOD<<endl;
}
void TC() {
	init(1000000);
	int tc=1;
	fastIO::read(tc);
	while(tc--) {
		Main();
		cout.flush();
	}
}
signed main() {
	return cin.tie(0),cout.tie(0),ios::sync_with_stdio(0),TC(),0;
}
/*
1. CLEAR the arrays (ESPECIALLY multitests)
2. DELETE useless output
 */

Try to find some way to "sort" the grid, for swapping different rows and different columns doesn't matter. How?

How to calculate the beauty?

How to solve the original problem when the initial grid has positive beauty? How to optimize?

How to solve the original problem when the initial grid is impossible to be painted out? Are there too many ways to flip a block?

Read the hints.

Let's sort the rows and columns in some way, which means you should swap two rows or two columns each time to make them ordered. Let's sort by the number of $$$1$$$ in each row and column. For example, the result of sorting the left part is shown in the right part:

We claim that the grid is solvealbe if and only if the later rows are "included" by previous rows, and so for the columns. For example, if a position in any row is filled with $$$1$$$, then the position in any previous row must also be $$$1$$$. This is correct because you'll notice that any submatrices of the following two types cannot appear in a solveable grid. To check that, you'll only need to check adjacent rows and columns.

After you sort all of that, the number of ways to paint it is easy to think about. It's similar to the product of factorials of the number of each kind of rows(for example, the last two rows in the right above). The only to notice is the first step of painting shouldn't be counted.

Consider modifying one block. If the initial grid is solveable, there're only $$$O(n)$$$ types of changing(in each row, all the $$$1$$$ and $$$2$$$ can be considered the same). Each type of changing can be done in $$$O(1)$$$ using simple discussion.

If the initial grid is unsolveable, we should find the two adjacent rows and columns that starts the problem. We claim that the block to change must belong to one of the four blocks in the intersection of the rows and columns. Just try every possible situation and do brute force on it.

For example, consider the following matrix: The first two columns and the last two rows went wrong. So let's try to change blocks $$$(3,1),(3,2),(4,1),(4,2)$$$ to judge whether it is possible.

The time complexity is $$$O((n+m)\log n)$$$, and if you implement carefully enough you can solve it in $$$O(n+m)$$$.

#if defined(LOCAL) or not defined(LUOGU)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast,unroll-loops")
#endif
#include<bits/stdc++.h>
using namespace std;
struct time_helper
{
#ifdef LOCAL
	clock_t time_last;
#endif
	time_helper()
	{
#ifdef LOCAL
		time_last=clock();
#endif
	}
	void test()
	{
#ifdef LOCAL
		auto time_now=clock();
		std::cerr<<"time:"<<1.*(time_now-time_last)/CLOCKS_PER_SEC<<";all_time:"<<1.*time_now/CLOCKS_PER_SEC<<std::endl;
		time_last=time_now;
#endif
	}
	~time_helper()
	{
		test();
	}
}time_helper;
#ifdef LOCAL
#include"dbg.h"
#else
#define dbg(...) (__VA_ARGS__)
#endif
namespace Fread{const int SIZE=1<<16;char buf[SIZE],*S,*T;inline char getchar(){if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return'\n';}return *S++;}}namespace Fwrite{const int SIZE=1<<16;char buf[SIZE],*S=buf,*T=buf+SIZE;inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;}inline void putchar(char c){*S++=c;if(S==T)flush();}struct NTR{~NTR(){flush();}}ztr;}
#define getchar Fread::getchar
#define putchar Fwrite::putchar
#define Setprecision 10
#define between '\n'
#define __int128 long long
template<typename T>struct is_char{static constexpr bool value=(std::is_same<T,char>::value||std::is_same<T,signed char>::value||std::is_same<T,unsigned char>::value);};template<typename T>struct is_integral_ex{static constexpr bool value=(std::is_integral<T>::value||std::is_same<T,__int128>::value)&&!is_char<T>::value;};template<typename T>struct is_floating_point_ex{static constexpr bool value=std::is_floating_point<T>::value||std::is_same<T,__float128>::value;};namespace Fastio{struct Reader{template<typename T>typename std::enable_if_t<std::is_class<T>::value,Reader&>operator>>(T&x){for(auto &y:x)*this>>y;return *this;}template<typename T>typename std::enable_if_t<is_integral_ex<T>::value,Reader&>operator>>(T&x){char c=getchar();short f=1;while(c<'0'||c>'9'){if(c=='-')f*=-1;c=getchar();}x=0;while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}x*=f;return *this;}template<typename T>typename std::enable_if_t<is_floating_point_ex<T>::value,Reader&>operator>>(T&x){char c=getchar();short f=1,s=0;x=0;T t=0;while((c<'0'||c>'9')&&c!='.'){if(c=='-')f*=-1;c=getchar();}while(c>='0'&&c<='9'&&c!='.')x=x*10+(c^48),c=getchar();if(c=='.')c=getchar();else return x*=f,*this;while(c>='0'&&c<='9')t=t*10+(c^48),s++,c=getchar();while(s--)t/=10.0;x=(x+t)*f;return*this;}template<typename T>typename std::enable_if_t<is_char<T>::value,Reader&>operator>>(T&c){c=getchar();while(c=='\n'||c==' '||c=='\r')c=getchar();return *this;}Reader&operator>>(char*str){int len=0;char c=getchar();while(c=='\n'||c==' '||c=='\r')c=getchar();while(c!='\n'&&c!=' '&&c!='\r')str[len++]=c,c=getchar();str[len]='\0';return*this;}Reader&operator>>(std::string&str){str.clear();char c=getchar();while(c=='\n'||c==' '||c=='\r')c=getchar();while(c!='\n'&&c!=' '&&c!='\r')str.push_back(c),c=getchar();return*this;}Reader(){}}cin;const char endl='\n';struct Writer{typedef __int128 mxdouble;template<typename T>typename std::enable_if_t<std::is_class<T>::value,Writer&>operator<<(T x){for(auto &y:x)*this<<y<<between;*this<<'\n';return *this;}template<typename T>typename std::enable_if_t<is_integral_ex<T>::value,Writer&>operator<<(T x){if(x==0)return putchar('0'),*this;if(x<0)putchar('-'),x=-x;static int sta[45];int top=0;while(x)sta[++top]=x%10,x/=10;while(top)putchar(sta[top]+'0'),--top;return*this;}template<typename T>typename std::enable_if_t<is_floating_point_ex<T>::value,Writer&>operator<<(T x){if(x<0)putchar('-'),x=-x;x+=pow(10,-Setprecision)/2;mxdouble _=x;x-=(T)_;static int sta[45];int top=0;while(_)sta[++top]=_%10,_/=10;if(!top)putchar('0');while(top)putchar(sta[top]+'0'),--top;putchar('.');for(int i=0;i<Setprecision;i++)x*=10;_=x;while(_)sta[++top]=_%10,_/=10;for(int i=0;i<Setprecision-top;i++)putchar('0');while(top)putchar(sta[top]+'0'),--top;return*this;}template<typename T>typename std::enable_if_t<is_char<T>::value,Writer&>operator<<(T c){putchar(c);return*this;}Writer&operator<<(char*str){int cur=0;while(str[cur])putchar(str[cur++]);return *this;}Writer&operator<<(const char*str){int cur=0;while(str[cur])putchar(str[cur++]);return*this;}Writer&operator<<(std::string str){int st=0,ed=str.size();while(st<ed)putchar(str[st++]);return*this;}Writer(){}}cout;}
#define cin Fastio::cin
#define cout Fastio::cout
#define endl Fastio::endl
 
void solve();
main()
{
	int t=1;
	cin>>t;
	while(t--)solve();
}
template <uint32_t mod>
struct LazyMontgomeryModInt {
	using mint = LazyMontgomeryModInt;
	using i32 = int32_t;
	using u32 = uint32_t;
	using u64 = uint64_t;
	static constexpr u32 get_r() {
		u32 ret = mod;
		for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
		return ret;
	}
	static constexpr u32 r = get_r();
	static constexpr u32 n2 = -u64(mod) % mod;
	static_assert(r * mod == 1, "invalid, r * mod != 1");
	static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30");
	static_assert((mod & 1) == 1, "invalid, mod % 2 == 0");
	u32 a;
	constexpr LazyMontgomeryModInt() : a(0) {}
	constexpr LazyMontgomeryModInt(const int64_t &b)
			: a(reduce(u64(b % mod + mod) * n2)){};
	static constexpr u32 reduce(const u64 &b) {
		return (b + u64(u32(b) * u32(-r)) * mod) >> 32;
	}
	constexpr mint &operator+=(const mint &b) {
		if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
		return *this;
	}
	constexpr mint &operator-=(const mint &b) {
		if (i32(a -= b.a) < 0) a += 2 * mod;
		return *this;
	}
	constexpr mint &operator*=(const mint &b) {
		a = reduce(u64(a) * b.a);
		return *this;
	}
	constexpr mint &operator/=(const mint &b) {
		*this *= b.inverse();
		return *this;
	}
	constexpr mint operator+(const mint &b) const { return mint(*this) += b; }
	constexpr mint operator-(const mint &b) const { return mint(*this) -= b; }
	constexpr mint operator*(const mint &b) const { return mint(*this) *= b; }
	constexpr mint operator/(const mint &b) const { return mint(*this) /= b; }
	constexpr bool operator==(const mint &b) const {
		return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
	}
	constexpr bool operator!=(const mint &b) const {
		return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
	}
	constexpr mint operator-() const { return mint() - mint(*this); }
	constexpr mint pow(u64 n) const {
		mint ret(1), mul(*this);
		while (n > 0) {
			if (n & 1) ret *= mul;
			mul *= mul;
			n >>= 1;
		}
		return ret;
	}
	constexpr mint inverse() const { return pow(mod - 2); }
	friend ostream &operator<<(ostream &os, const mint &b) {
		return os << b.get();
	}
	friend istream &operator>>(istream &is, mint &b) {
		int64_t t;
		is >> t;
		b = LazyMontgomeryModInt<mod>(t);
		return (is);
	}
	constexpr u32 get() const {
		u32 ret = reduce(a);
		return ret >= mod ? ret - mod : ret;
	}
	static constexpr u32 get_mod() { return mod; }
	explicit operator u32() const {return get();}
};
#define modint LazyMontgomeryModInt
#define mint modint<p>
constexpr int p=998244353;
std::pmr::monotonic_buffer_resource pool(100000000);
void solve()
{
	int n,m;
	cin>>n>>m;
	vector<std::pmr::vector<int>>a(n,std::pmr::vector<int>{&pool});
	time_helper.test();
	while(m--)
	{
		int u,v;
		cin>>u>>v;
		u--,v--;
		a[u].emplace_back(v);
	}
	time_helper.test();
	std::pmr::vector<mint>inv{&pool},fac{&pool},ifac{&pool};
	inv.resize(n+1),fac.resize(n+1),ifac.resize(n+1);
	std::pmr::vector<int>arr_cnt{&pool};
	arr_cnt.resize(n+2);
	auto arr_sort=[&](const auto&a)
	{
		for(int i=0;i<n;i++)
		arr_cnt[a[i].size()+1]++;
		for(int x=1;x<=n;x++)
		arr_cnt[x]+=arr_cnt[x-1];
		vector<std::pmr::vector<int>>b(n,std::pmr::vector<int>{&pool});
		for(int i=0;i<n;i++)
		b[arr_cnt[a[i].size()]++]=move(a[i]);
		for(int x=0;x<=n+1;x++)
		arr_cnt[x]=0;
		return b;
	};
	a=move(arr_sort(a));
	auto arr_eq=[&](const auto&a,const auto&b)
	{
		if(a.size()!=b.size())return false;
		bool res=1;
		for(auto q:a)arr_cnt[q]=1;
		for(auto q:b)if(!arr_cnt[q]){res=0;break;}
		for(auto q:a)arr_cnt[q]=0;
		return res;
	};
	auto arr_includes=[&](const auto&a,const auto&b)
	{
		bool res=1;
		for(auto q:a)arr_cnt[q]=1;
		for(auto q:b)if(!arr_cnt[q]){res=0;break;}
		for(auto q:a)arr_cnt[q]=0;
		return res;
	};
	auto arr_set_difference=[&](const auto&a,const auto&b,auto &diff)
	{
		for(auto q:b)arr_cnt[q]=1;
		for(auto q:a)
		if(!arr_cnt[q])
		{
			diff.emplace_back(q);
			if(diff.size()>2)break;
		}
		for(auto q:b)arr_cnt[q]=0;
	};
	{
		fac[0]=1;
		for(int x=1;x<=n;x++)
		fac[x]=fac[x-1]*x;
		ifac[n]=fac[n].inverse();
		for(int x=n;x>=1;x--)
		ifac[x-1]=ifac[x]*x,inv[x]=ifac[x]*fac[x-1];
	}
	int failcnt=0;
	for(int x=1;x<n;x++)
	failcnt+=!(arr_includes(a[x],a[x-1]));
	if(failcnt==0)
	{
		int nc=0;
		mint ansbas=fac[a[0].size()],ans=0;
		std::pmr::vector<int>tnc{&pool},tmc{&pool};
		tmc.emplace_back(a[0].size());
		int ct=0;
		for(int x=0;x<n;x++)
		if(x==n-1||!arr_eq(a[x],a[x+1]))
		{
			nc++;
			if(a[x].size()!=n)
			ansbas*=fac[nc];
			int mc=0;
			if(x==n-1)mc=n-a[x].size();
			else mc=a[x+1].size()-a[x].size();
			if(x!=n-1)
			ansbas*=fac[mc];
			tnc.emplace_back(nc);
			tmc.emplace_back(mc);
			nc=0;
			ct++;
		}
		else nc++;
		for(int x=0;x<ct;x++)
		{
			int nc=tnc[x],mc=tmc[x+1];
			if(mc==0)continue;
			mint fix=ansbas;
			if(x==ct-1)
			{
				fix*=mc;
			}
			if(x==ct-2&&mc==1&&tmc[x+2]==0)
			{
				fix*=inv[tnc[x+1]+1];
			}
			if(x==ct-1||mc!=1)
			{
				if(mc)
				{
					if(nc>1)ans+=fix;
					else ans+=fix*(tmc[x]+1);
				}
			}
			else if(nc>1)ans+=fix*(tnc[x+1]+1);
			else ans+=fix*(tnc[x+1]+1)*(tmc[x]+1);
		}
		for(int x=0;x<ct;x++)
		{
			int nc=tnc[x],mc=tmc[x];
			if(mc==0)continue;
			mint fix=ansbas;
			if(x==ct-1)
			{
			    if(tmc[x+1]==0)fix*=nc;
			    if(nc==1)fix*=inv[tmc[x+1]+1];
			}
			if(x==0||mc!=1)
			{
				if(mc)
				{
					if(nc>1)ans+=fix;
					else ans+=fix*(tmc[x+1]+1);
				}
			}
			else if(nc>1)ans+=fix*(tnc[x-1]+1);
			else ans+=fix*(tnc[x-1]+1)*(tmc[x+1]+1);
		}
		cout<<(ans*inv[n]*inv[n]).get()<<endl;
		return;
	}
	if(failcnt>2){cout<<0<<endl;return;}
	mint ans=0;
    auto calc=[&](vector<std::pmr::vector<int>>&a)->mint
	{
        a=move(arr_sort(a));
		for(int x=1;x<n;x++)
		if(!arr_includes(a[x],a[x-1]))return 0;
		int nc=0;
		mint ansbas=fac[a[0].size()];
		for(int x=0;x<n;x++)
		if(x==n-1||!arr_eq(a[x],a[x+1]))
		{
			nc++;
			if(a[x].size()!=n)
			ansbas=ansbas*fac[nc];
			int mc=0;
			if(x==n-1)mc=n-a[x].size();
			else mc=a[x+1].size()-a[x].size();
			if(x!=n-1)
			ansbas=ansbas*fac[mc];
			nc=0;
		}
		else nc++;
		return ansbas;
	};
	std::pmr::set<pair<int,int>>st{&pool};
	std::pmr::vector<int>diff{&pool};
	vector<std::pmr::vector<int>>b(n,std::pmr::vector<int>{&pool});
	time_helper.test();
	for(int x=1;x<n;x++)
	if(!arr_includes(a[x],a[x-1]))
	{
		for(int y=x+2;y<n;y++)
		if(!arr_includes(a[y],a[y-1]))goto fail;
		{
			diff.clear();
			arr_set_difference(a[x],a[x-1],diff);
			if(diff.size()<=1)
			{
				if(!st.count({x-1,diff[0]}))
				{
					b=a;
					b[x-1].emplace_back(diff[0]);
					ans+=calc(b),st.insert({x-1,diff[0]});
				}
				if(!st.count({x,diff[0]}))
				{
					b=a;
					b[x].erase(find(b[x].begin(),b[x].end(),diff[0]));
					ans+=calc(b),st.insert({x,diff[0]});
				}
			}
		}
		{
			diff.clear();
			arr_set_difference(a[x-1],a[x],diff);
			if(diff.size()<=1)
			{
				if(!st.count({x,diff[0]}))
				{
					b=a;
					b[x].emplace_back(diff[0]);
					ans+=calc(b),st.insert({x,diff[0]});
				}
				if(!st.count({x-1,diff[0]}))
				{
					b=a;
					b[x-1].erase(find(b[x-1].begin(),b[x-1].end(),diff[0]));
					ans+=calc(b),st.insert({x-1,diff[0]});
				}
			}
		}
		fail:
		break;
	}
	cout<<(ans*inv[n]*inv[n]).get()<<endl;
}