Codechef Oct chanllenge Queries on Matrix-JIIT

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

Codechef Oct chanllenge Queries on Matrix-JIIT

首先发现矩阵的两个维度显然是互不相干的，假设最后操作后有 \(x\) 列被操作奇数次， \(y\) 行操作奇数次

那么最后为奇数的格子个数就是 \(x(m-y)+(n-x)y\)

考虑求出 \(q\) 操作后有 \(x\) 个位置被操作奇数次的方案数

考虑一个Naive的 \(dp\) ，令 \(dp_{i,j}\) 为操作 \(i\) 次后有 \(j\) 位置操作奇数的方案数，显然得到转移为

\(dp_{i,j}\cdot j\rightarrow dp_{i+1,j-1}\)

\(dp_{i,j}\cdot (n-j)\rightarrow dp_{i+1,j+1}\)

直接 \(dp\) 复杂度为 \(O(nq)\) ，用矩阵优化复杂度为 \(O(n^3\log q)\)

没有考虑过求向量列的现行递推式？说不定暴力求递推然后。。。

\(O(n^3)\)

~~由于笔者不会数学，所以~~考虑一个非常暴力非常直观的理解，可以完全抛开组合意义

每次是挑选一个位置异或上1，用形式幂级数可能比较蛋疼，不如直接搞成集合幂级数

即令集合 \(S\) 包含所有被操作奇数次的位置，用一个多项式 \(F=\sum_S a_S x^S\) 表示答案

那么转移多项式即为 \(F=\sum_{i=1}^{n} x^{\lbrace i \rbrace}\) ，转移运算为集合对称差运算~~(哎就是异或)~~

那么实际就是要求出 \(F^q\) ，可以直接用 \(\text{FWT}\) 优化，先 \(\text{FWT}\) ，然后求出每一项的 \(q\) 次幂，然后 \(\text{FWT}\) 回来

~~(这样不是 \(n2^n\) 的吗)~~

显然的可以发现， \(F^i\) 的任意位置系数 \([x^S]F^i\) 只与 \(|S|\) 有关，所以实际上只需要存下含有 \(0,1,2\cdots,n\) 个元素的项的系数即可

考虑在这样的多项式上模拟原先的 \(\text{FWT}\) 过程

按照快速沃尔什变换的式子，令 \(G=\text{FWT(F)}=\sum_{S}x^S \sum_{T}(-1)^{|S\cap T|}\cdot [x^T]F\)

考虑枚举 \(|S|,|T|,|S\cap T|\) ，然后组合数计算系数，令 \(F_i\) 表示 \([x^S]F(|S|=i)\)

则有 \(\begin{aligned} G_i=\sum_{j}F_j\sum_k C_i^k\cdot C_{n-i}^{j-k}(-1)^{k}\end{aligned}\) ，其中 \(C_{n-i}^{j-k}\) 表示从不相交的部分里选出 \(j\) 中剩下的元素

按照该式即可完成 \(O(n^3)\) 模拟 \(\text{FWT}\) ，注意 \(\text{IFWT}\) 时需要除掉 \(2^n\)

算上快速幂复杂度应为 \(O(n^3+n\log q)\)

\(O(n^2\log n)\)

用 \(\text{NTT}\) 优化上式

\(O(n^2)\)

依然考虑上面 \(\text{FWT}\) 的转移式，发现 \(i\) 项得到 \(j\) 项的贡献为一个常数，考虑直接计算这个常数 \(W_{i,j}=\sum_k C_i^k\cdot C_{n-i}^{j-k}(-1)^{k}\)

我们知道组合数就是二项展开的结果，所以发现实际就是 \(W_{i,j}=[x^j] (-x+1)^i\cdot (x+1)^{n-i}\)

该式可以 \(O(n^2)\) 按照 \(i\) 递推求出，每次乘上一个 \(\frac{1-x}{1+x}\) 即可

实际复杂度为 \(O(n^2+n\log q)\)

实际上应该还可以处理一些稍微复杂点的问题，比如每次可以操作若干个位置

不知道能否优化到 \(n^2\) 以下

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

typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef long double ldb;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

bool Mbe;

const int N=2050,P=998244353;

int n,m,Z; ll k;
int A[N],B[N];
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

int C[N][N];
int rev[N],T[N],IT[N];
void NTT(int n,int *a,int f){
	rep(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	static int e[N>>1];	
	for(reg int i=e[0]=1;i<n;i<<=1) {
		ll t=f==1?T[i]:IT[i];
		//qpow(f==1?3:(P+1)/3,(P-1)/i/2);
		for(reg int j=i-2;j>=0;j-=2) e[j+1]=t*(e[j]=e[j>>1])%P;
		for(reg int l=0;l<n;l+=i*2) {
			for(reg int j=l;j<l+i;++j) {
				reg int t=1ll*a[j+i]*e[j-l]%P;
				a[j+i]=a[j]-t,Mod2(a[j+i]);
				a[j]+=t,Mod1(a[j]);
			}
		}
	}
	if(f==-1) {
		ll base=qpow(n);
		rep(i,0,n-1) a[i]=a[i]*base%P;
	}
}

int Init(int n){
	int R=1,cc=-1;
	while(R<=n) R<<=1,cc++;
	rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
	return R;
}


int mk1[N],mk2[N];
void FWT(int n,int *A,int *B,int f,int *mk){
		static int X[N],Y[N];
		rep(i,0,n) {
			memset(X,0,sizeof X),memset(Y,0,sizeof Y);
			rep(k,0,i) X[i-k]=(k&1)?P-C[i][k]:C[i][k];
			rep(j,0,n) Y[j]=A[j];
			int R=Init(n+i+1);
			NTT(R,X,1),NTT(R,Y,1);
			rep(j,0,R-1) X[j]=1ll*X[j]*Y[j]%P;
			NTT(R,X,-1);
			rep(j,0,n-i) B[i]=(B[i]+1ll*C[n-i][j]*X[i+j])%P;
		}
}

void Solve(int n,int *A,int *mk) {
	static int X[N],Y[N];
	memset(X,0,sizeof X),memset(Y,0,sizeof Y);
	X[1]=1,FWT(n,X,Y,1,mk);
	rep(i,0,n) Y[i]=qpow(Y[i],k);
	memset(X,0,sizeof X);
	FWT(n,Y,X,2,mk);
	ll base=qpow(qpow(2,n),P-2);
	rep(i,0,n) A[i]=X[i]*base%P*C[n][i]%P;
}


bool Med;
int main(){
	//fprintf(stderr,"%.2lf\n",(&Med-&Mbe)/1024.0/1024.0);
	freopen("clone.in","r",stdin),freopen("clone.out","w",stdout);
	rep(i,0,N-1) rep(j,C[i][0]=1,i) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
	for(int i=1;i<N;i<<=1) {
		T[i]=qpow(3,(P-1)/i/2);
		IT[i]=qpow((P+1)/3,(P-1)/i/2);
	}
	scanf("%d%d%lld%d",&n,&m,&k,&Z);
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		mk1[i]=1,mk2[j]=1;
	}
	Solve(n,A,mk1),Solve(m,B,mk2);
	int ans=0;
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		ans=(ans+1ll*A[i]*B[j])%P;
	}
	printf("%d\n",ans);
}


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

typedef long long ll;
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

const int N=2050,P=998244353;

int n,m,Z; ll k;
int A[N],B[N];
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int C[N][N],W[N][N],mk1[N],mk2[N];

void Solve(int n,int *A,int *mk) {
	static int X[N],Y[N];
	rep(i,0,n) W[0][i]=C[n][i];
	rep(i,1,n) {
		rep(j,0,n) W[i][j]=W[i-1][j]-(j?W[i][j-1]:0),Mod2(W[i][j]);
		drep(j,n,0) W[i][j]=(j?-W[i][j-1]:0)+W[i][j],Mod2(W[i][j]);
	}

	memset(Y,0,sizeof Y);
	rep(i,0,n) X[i]=qpow(W[i][1],k);
	rep(i,0,n) if(mk[i]) rep(j,0,n) Y[i]=(Y[i]+1ll*W[i][j]*X[j])%P;
	ll base=qpow(qpow(2,n),P-2);
	rep(i,0,n) A[i]=Y[i]*base%P*C[n][i]%P;
}

int main(){
	freopen("clone.in","r",stdin),freopen("clone.out","w",stdout);
	scanf("%d%d%lld%d",&n,&m,&k,&Z);
	rep(i,0,max(n,m)) rep(j,C[i][0]=1,i) C[i][j]=C[i-1][j-1]+C[i-1][j],Mod1(C[i][j]);
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		mk1[i]=mk2[j]=1;
	}
	Solve(n,A,mk1),Solve(m,B,mk2);
	int ans=0;
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		ans=(ans+1ll*A[i]*B[j])%P;
	}
	printf("%d\n",ans);
}

CodeChef 2020 November Challenge - Scalar Product Tree (莫队)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

CodeChef 2020 November Challenge - Scalar Product Tree (莫队)

题目大意：给定一棵根为1的树，每个点有权值 \(A_i\) ，每个点按照其从根开始的路径记录下来一串数 \((A_1,\cdots,A_u)\) 构成一个向量 \(v_u\)

每次查询两个点 \((x,y)\) ，查询 \(v_x\cdot v_y\)

\[ \ \]

由于向量点积是对位相乘的，不好用数据结构维护

考虑用一个序列描述 \(\text{dfs}\) 遍历树时每个点入栈出栈的过程，扫描一段前缀即可得到遍历到每个点时 \(\text{dfs}\) 栈的情况，也就得到了题目指定的向量

每次查询两个点 \(x,y\) ，那么就是查询了两段前缀，用莫队维护两个前缀指针的移动，同时维护每个深度上两个前缀对应的值以及这些值的乘积即可

复杂度为 \(O(n\sqrt n)\)

#include<bits/stdc++.h>
using namespace std;
typedef pair <int,int> Pii;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

char IO;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) f|=IO=='-';
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=6e5+10,P=1e9+7;

int n,m;
struct Edge{
	int to,nxt;
}e[N];
int head[N],ecnt;
void AddEdge(int u,int v) {
	e[++ecnt]=(Edge){v,head[u]};
	head[u]=ecnt;
}
int id[N],dfn;

typedef unsigned U;
U A[N];
int L[N],K[N]; // L,K维护括号序列，L为编号,K为左括号还是右括号
U Sum,Ans[N];
int dep[N];

void dfs(int u,int f) {
	L[id[u]=++dfn]=u,K[dfn]=1;
	dep[u]=dep[f]+1;
	for(int i=head[u];i;i=e[i].nxt) {
		int v=e[i].to;
		if(v==f) continue;
		dfs(v,u);
	}
	L[++dfn]=u,K[dfn]=-1;
}

int len;
struct Que{
	int l,r,id;
	bool operator < (const Que __) const {
		if(l/len!=__.l/len) return l/len<__.l/len;
		return ((l/len)&1)?r<__.r:r>__.r;
	}
} Q[N];

U S[2][N];
void Add(int d,int x,int k) {
	int p=dep[x];
	Sum-=S[d][p]*S[!d][p];
	S[d][p]=k==1?A[x]:0;
	Sum+=S[d][p]*S[!d][p];
}

int main() {
	n=rd(),m=rd();
	rep(i,1,n) A[i]=rd();
	rep(i,2,n) {
		int u=rd(),v=rd();
		AddEdge(u,v),AddEdge(v,u);
	}
	dfs(1,0),len=sqrt(dfn);
	rep(i,1,m) {
		int l=rd(),r=rd();
		l=id[l],r=id[r];
		if(l>r) swap(l,r);
		Q[i]=(Que){l,r,i};
	}
	sort(Q+1,Q+m+1);
	int l=1,r=1;
	Add(0,1,1),Add(1,1,1);
	rep(i,1,m) {
		while(l<Q[i].l) ++l,Add(0,L[l],K[l]);
		while(l>Q[i].l) Add(0,L[l],-K[l]),l--;
		
		while(r<Q[i].r) ++r,Add(1,L[r],K[r]);
		while(r>Q[i].r) Add(1,L[r],-K[r]),r--;
		Ans[Q[i].id]=Sum;
	}
	rep(i,1,m) printf("%u\n",Ans[i]);
}

Codechef November Chanllenge 2019 Div1 PrettyBox (贪心，线段树)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

Codechef November Chanllenge 2019 Div1 PrettyBox (贪心，线段树)

原题链接

前言：这篇文章主要讲如何用线段树优化贪心，关于贪心的证明建议看官方题解

贪心思路：

首先肯定要按照 \((S_i,P_i)\) 递增的顺序排序

每次选取两个点，一个标记为左括号，权值为 \(-P_i\) ，一个标记为右括号，权值为 \(P_i\) ，显然只要是一个合法的括号序列即可

题解证明了在不断增加括号时，不会出现一个位置的括号情况改变

现在我们的贪心问题就在于怎样找到一对最优的括号，注意每次选出的 两个括号之间并不一定匹配

为了便于描述，把左括号看做1，右括号看做-1，一个合法括号序列满足任何一个前缀和 \(\ge 0\)

考虑什么样的情况可以放置左右括号，设分别放在 \(x,y\)

\(x<y\) 显然合法
\(x>y\) 时，如果存在一个括号对，将 \((x,y)\) 包含在一起，即 \((y,x)\) 这一段区间不跨过一个前缀和为 \(0\) 的位置

如果把序列看做由一段段 前缀和为0的位置 分割开来的一个个联通块，似乎比较好理解

也就是块内随意选，之间只能由小到大匹配

接下来考虑用线段树维护这样的块的信息，下面只讨论 \(x>y\) 的情况

由于线段树每个结点统计区间 \([l,r]\) 的信息，所以实际上块之间的间隔并不为0

设 \([l,r]\) 中最小的前缀和为 \(Min\) (是指从 \(l\) 开始的前缀和)

不妨统计 \([l,r]\) 中不跨过一个前缀和为 \(Min\) 的位置的答案 \(Ans\) ，以及跨过的答案 \(Ans2\)

合并两个区间时，需要找到

左区间中右边连续的一段不跨过最小值的最大权值 \(R\)

右区间中左边连续的一段不跨过最小值的最小权值 \(L\)

以及任意的最小值最大值 \(mi,ma\)

然后按照 \(Min\) 的权值大小关系，判断这4种权值的合并应该被分配到 \(Ans\) 还是 \(Ans2\)

合并 \(L,R\) 时注意 \(L\) 优先看左儿子， \(R\) 优先看右儿子，具体实现看代码中的 \(Up\) 函数

每次存下答案找到最优配对后，在序列上对应放置-1,1单点修改即可，复杂度为 \(O(n\log n)\)

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

typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ a=min(a,b); }
template <class T> inline void cmax(T &a,T b){ a=max(a,b); }

char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=2e5+10,M=N<<2;

int n;
Pii T[N];
int A[N];

int S[M]; // 区间和
int Min[M]; // 前缀最小值
struct Node{
	int x;
	Node(){}
	Node(int x):x(x){}
	int operator < (const Node __) const { return A[x]<A[__.x]; }
} L[M],R[M]; // 左最小，右最大  , 记录的是不跨过最小值的权值
Node mi[M],ma[M]; // 区间最大最小，没有限制

struct Pair{
	int x,y;
	Pair(){}
	Pair(int x,int y):x(x),y(y){}
	Pair(Node x,Node y):x(x.x),y(y.x){}
	int Val() const { return A[y]-A[x]; }
	int operator < (const Pair __) const { return Val()<__.Val(); }
} Ans[M],Ans2[M]; // 不跨过最小值的答案以及x<y的答案，包含最小值的答案
// 区间答案

void Up(int p){
	S[p]=S[p<<1]+S[p<<1|1];
	Min[p]=min(Min[p<<1],S[p<<1]+Min[p<<1|1]);
	mi[p]=min(mi[p<<1],mi[p<<1|1]),ma[p]=max(ma[p<<1],ma[p<<1|1]);

	Ans[p]=max(Ans[p<<1],Ans[p<<1|1]);
	Ans2[p]=Pair(mi[p<<1|1],ma[p<<1]);
	cmax(Ans2[p],Ans2[p<<1]);
	cmax(Ans2[p],Ans2[p<<1|1]);

	cmax(Ans[p],Pair(mi[p<<1],ma[p<<1|1]));
	Ans[p]=max(Ans[p],Pair(L[p<<1|1],R[p<<1]));

	if(Min[p<<1]!=Min[p]) {
		cmax(Ans[p],Ans2[p<<1]);
		cmax(Ans[p],Pair(L[p<<1|1],ma[p<<1]));
		L[p]=min(mi[p<<1],L[p<<1|1]);
	} else {
		L[p]=L[p<<1];
	}

	if(S[p<<1]+Min[p<<1|1]!=Min[p]) {
		cmax(Ans[p],Ans2[p<<1|1]);
		cmax(Ans[p],Pair(mi[p<<1|1],R[p<<1]));
		R[p]=max(R[p<<1],ma[p<<1|1]);
	} else {
		R[p]=R[p<<1|1];
	}

}

void Build(int p,int l,int r){
	if(l==r) {
		S[p]=Min[p]=0;
		L[p]=n+1,R[p]=n+2;
		mi[p]=ma[p]=l;
		Ans[p]=Ans2[p]=Pair(n+1,n+2);
		return;
	}
	int mid=(l+r)>>1;
	Build(p<<1,l,mid),Build(p<<1|1,mid+1,r);
	Up(p);
}
void Upd(int p,int l,int r,int x,int k) {
	if(l==r) {
		S[p]=Min[p]=k;
		L[p]=n+1,R[p]=n+2;
		mi[p]=n+1,ma[p]=n+2;
		Ans[p]=Ans2[p]=Pair(n+1,n+2);
		return;
	}
	int mid=(l+r)>>1;
	x<=mid?Upd(p<<1,l,mid,x,k):Upd(p<<1|1,mid+1,r,x,k);
	Up(p);
}

int main(){
	rep(i,1,n=rd()) T[i].first=rd(),T[i].second=rd();
	sort(T+1,T+n+1);
	rep(i,1,n) A[i]=T[i].second;
	A[n+1]=1e9+10,A[n+2]=-1e9-10;
	ll ans=0;
	int i=1;
	Build(1,1,n);
	while(i<=n/2) {
		Pair res=Ans[1];
		if(res.Val()<=0) break;
		printf("%lld\n",ans+=res.Val());
		Upd(1,1,n,res.x,1),Upd(1,1,n,res.y,-1);
		i++;
	}
	while(i<=n/2) printf("%lld\n",ans),i++;
}

CodeChef 2020 November - Challenge Chef and the Combination Lock (多项式)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

CodeChef 2020 November - Challenge Chef and the Combination Lock (多项式)

题目大意：给定了 \(n\) 个随机变量 \(x_i\in{0,1,\cdots,A_i}\) ，令 \(\Chi=\min_i\lbrace x_i,A_i-x_i\rbrace\) ，求 \(E(\Chi)\)

我们知道 \(E(\Chi)=\sum_{i=0}^{\infty} P(\Chi>i)\)

不妨考虑计算 \(P(\Chi>i)\) ，先计算方案数，发现方案数可以用一个多项式来表示

令 \(F(x)\) 为 \(\Chi>x\) 的方案数，则 \(\begin{aligned}F(x)=\prod_{i=1}^n (A_i-1-2x)　\end{aligned}\)

显然 \(\Chi \leq \min\lbrace \frac{A_i}{2} \rbrace\) ，不妨设这个上界为 \(U\) ，也就是说我们要求 \(\sum_{i=0}^U F(i)\)

常识：一个 \(n\) 次多项式前缀和可以用一个不超过 \(n+1\) 次的多项式来表示

如果暴力求出 \(F(x)\) 在 \(x=0,1,\cdots,n+1\) 处的值，累前缀和，然后用拉格朗日插值法求出解

暴力求值复杂度为 \(O(n^2)\) ，拉格朗日插值复杂度为 \(O(n)\)

可以用分治 \(\text{NTT}\) 优化 \(F(x)\) 的求解，然后用多项式多点求值求得点值

复杂度为 \(O(n\log ^2n)\) ，实际在CodeChef上的运行时间为0.53s

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

typedef long long ll;
typedef double db;
typedef unsigned long long ull;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=1<<18,P=998244353;


ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

typedef vector <int> Poly;
void Show(Poly a,int k=0){ 
	if(!k){ for(int i:a) printf("%d ",i); puts(""); }
	else for(int i:a) printf("%d\n",i);
}
int rev[N],w[N];
int Inv[N+1],Fac[N+1],FInv[N+1];

void Init() { 
	int t=qpow(3,(P-1)/N);
	w[N>>1]=1;
	rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
	drep(i,(N>>1)-1,1) w[i]=w[i<<1];
	Inv[0]=Inv[1]=Fac[0]=Fac[1]=FInv[0]=FInv[1]=1;
	rep(i,2,N) {
		Inv[i]=1ll*(P-P/i)*Inv[P%i]%P; 
		FInv[i]=1ll*FInv[i-1]*Inv[i]%P;
		Fac[i]=1ll*Fac[i-1]*i%P;
	}

}
int Init(int n){
	int R=1,c=-1;
	while(R<n) R<<=1,c++;
	rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
	return R;
}


void NTT(int n,int *a,int f){
	rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
	for(int i=1;i<n;i<<=1) {
		int *e=w+i;
		for(int l=0;l<n;l+=i*2) {
			for(int j=l;j<l+i;++j) {
				int t=1ll*a[j+i]*e[j-l]%P;
				a[j+i]=a[j]-t,Mod2(a[j+i]);
				a[j]+=t,Mod1(a[j]);
			}
		}
	}
	if(f==-1) {
		reverse(a+1,a+n);
		ll base=Inv[n];
		rep(i,0,n-1) a[i]=a[i]*base%P;
	}
}
void NTT(int n,Poly &a,int f){
	static int A[N];
	if((int)a.size()<n) a.resize(n);
	rep(i,0,n-1) A[i]=a[i];
	NTT(n,A,f);
	rep(i,0,n-1) a[i]=A[i];
}

Poly operator * (Poly a,Poly b){
	int n=a.size()+b.size()-1;
	int R=Init(n);
	a.resize(R),b.resize(R);
	NTT(R,a,1),NTT(R,b,1);
	rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
	NTT(R,a,-1);
	a.resize(n);
	return a;
}

Poly operator + (Poly a,Poly b) { 
	int n=max(a.size(),b.size());
	a.resize(n),b.resize(n);
	rep(i,0,n-1) a[i]+=b[i],Mod1(a[i]);
	return a; 
}
Poly operator - (Poly a,Poly b) { 
	int n=max(a.size(),b.size());
	a.resize(n),b.resize(n);
	rep(i,0,n-1) a[i]-=b[i],Mod2(a[i]);
	return a; 
}

Poly Poly_Inv(Poly a) { 
	int n=a.size();
	if(n==1) return Poly{(int)qpow(a[0],P-2)};
	Poly b=a; b.resize((n+1)/2); b=Poly_Inv(b);
	int R=Init(n<<1);
	a.resize(R),b.resize(R);
	NTT(R,a,1),NTT(R,b,1);
	rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
	NTT(R,a,-1);
	a.resize(n);
	return a;
}

// 应用转置原理优化的多项式多点求值
Poly Evaluate(Poly F,Poly X){
	static int ls[N<<1],rs[N<<1],cnt;
	static Poly T[N<<1];
	static auto TMul=[&] (Poly F,Poly G){
		int n=F.size(),m=G.size();
		if(n<=20 && m<=20){
			rep(i,0,n-m) {
				int t=0;
				rep(j,0,m-1) t=(t+1ll*F[i+j]*G[j])%P;
				F[i]=t;
			} 
			F.resize(n-m+1);
			return F;
		}
		reverse(G.begin(),G.end());
		int R=Init(n);
		NTT(R,F,1),NTT(R,G,1);
		rep(i,0,R-1) F[i]=1ll*F[i]*G[i]%P;
		NTT(R,F,-1); Poly T(n-m+1);
		rep(i,0,n-m) T[i]=F[i+m-1];
		return T;
	};
	static function <int(int,int)> Build=[&](int l,int r) {
		int u=++cnt; ls[u]=rs[u]=0;
		if(l==r) {
			T[u]=Poly{1,P-X[l]};
			return u;
		}
		int mid=(l+r)>>1;
		ls[u]=Build(l,mid),rs[u]=Build(mid+1,r);
		T[u]=T[ls[u]]*T[rs[u]];
		return u;
	};

	int n=F.size(),m=X.size();
	cmax(n,m),F.resize(n),X.resize(n);
	cnt=0,Build(0,n-1);
	F.resize(n*2+1),T[1]=TMul(F,Poly_Inv(T[1]));
	int p=0;
	rep(i,1,cnt) if(ls[i]) {
		swap(T[ls[i]],T[rs[i]]);

		int R=Init(T[i].size()),n=T[i].size(),m1=T[ls[i]].size(),m2=T[rs[i]].size();
		NTT(R,T[i],1);
		reverse(T[ls[i]].begin(),T[ls[i]].end()); reverse(T[rs[i]].begin(),T[rs[i]].end());
		NTT(R,T[ls[i]],1); NTT(R,T[rs[i]],1);
		rep(j,0,R-1) {
			T[ls[i]][j]=1ll*T[ls[i]][j]*T[i][j]%P;
			T[rs[i]][j]=1ll*T[rs[i]][j]*T[i][j]%P;
		}
		NTT(R,T[ls[i]],-1); NTT(R,T[rs[i]],-1);
		rep(j,0,n-m1) T[ls[i]][j]=T[ls[i]][j+m1-1];
		T[ls[i]].resize(n-m1+1);
		rep(j,0,n-m2) T[rs[i]][j]=T[rs[i]][j+m2-1];
		T[rs[i]].resize(n-m2+1);
	} else X[p++]=T[i][0];
	X.resize(m);
	return X;
}

int n;
int A[N];
int I[N],J[N];
int F[N],L[N],R[N];

Poly Solve(int l,int r) {
	if(l==r) return Poly{A[l]-1,P-2};
	int mid=(l+r)>>1;
	return Solve(l,mid)*Solve(mid+1,r);
}


int main() {
	Init();
	rep(i,J[0]=1,N-1) J[i]=1ll*J[i-1]*i%P;
	I[N-1]=qpow(J[N-1]);
	drep(i,N-1,1) I[i-1]=1ll*I[i]*i%P;

	rep(kase,1,rd()) {
		int x=P,All=1; 
		n=rd();
		rep(i,1,n+1) F[i]=0;
		F[0]=1;
		int f=0;
		rep(i,1,n) {
			A[i]=rd();
			f|=!A[i];
			cmin(x,(A[i]-1)/2),All=1ll*All*(A[i]+1)%P;
		}
		if(f){ puts("0"); continue; }
		Poly Y=Solve(1,n),X(n+2);
		rep(i,0,n+1) X[i]=i;
		Y=Evaluate(Y,X);
		rep(i,1,n+1) Y[i]=(Y[i]+Y[i-1])%P;
		int ans=0;
		if(x<=n+1) ans=Y[x];
		else {
            // 拉格朗日插值
			L[0]=x;
			rep(i,1,n+1) L[i]=1ll*L[i-1]*(x-i)%P;
			R[n+2]=1;
			drep(i,n+1,0) R[i]=1ll*R[i+1]*(x-i)%P;
			rep(i,0,n+1) {
				int t=1ll*Y[i]*(i?L[i-1]:1)%P*R[i+1]%P*I[i]%P*I[n+1-i]%P;
				if((n+1-i)&1) t=P-t;
				ans=(ans+t)%P;
			}
		}
		ans=ans*qpow(All)%P;
		printf("%d\n",ans);
	}
}

Codechef March Challenge 2021 Div2 Consecutive Adding(CONSADD)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

Codechef March Challenge 2021 Div2 Consecutive Adding(CONSADD)

题目大意：

给定两个 \(n\times m\) 矩阵 \(A\) ， \(B\) 和一个常数 \(x\)

现在对于 \(A\) 操作，每次可以选择一行或者一列连续的 \(x\) 个，一起改变同一个数值 \(v\in \Z\)

判断是否可以由 \(A\) 变成 \(B\)

显然可以先将 \(A,B\) 作差，转化为操作成0矩阵

进一步，我们将 \(A\) 矩阵行内差分，使得每次行操作变为一个单点 \(A_{i,j}+v\) ，一个单点 \(A_{i,j+x}-v\)

在此基础上，继续差分即可将行列操作都转化为单点操作

此时容易发现， \(A_{i,j}\) 的数值有关联的部分都是 \(A_{i,j},A_{i+x,j},A_{i,j+x}\cdots A_{i+ax,j+bx}\)

也就是相差 \(x\) 的，考虑可以将这一部分子矩形提取出来，这样问题变成了

每次操作一个数 \(A_{i,j}+v\) ，可以选择相邻一个数 \(A_{i,j+1}\) 或 \(A_{i+1,j}\) 去 \(-v\)

对于每个这样的子问题，容易发现有解的充要条件：子矩阵元素和为0

（可以依次考虑每个元素贪心构造方案）

如此可以 \(O(nm)\) 判定


const int N=1010,INF=1e9+10;

int n,m,k;
ll A[N][N],B[N][N];
int V[N][N];

int main(){
	rep(kase,1,rd()) {
		n=rd(),m=rd(),k=rd();
		rep(i,1,n+1) rep(j,1,m+1) A[i][j]=V[i][j]=0;
		rep(i,1,n) rep(j,1,m) A[i][j]=rd();
		rep(i,1,n) rep(j,1,m) A[i][j]-=rd();
		rep(i,1,n+1) drep(j,m+1,1) A[i][j]-=A[i][j-1];
		drep(i,n+1,1) rep(j,1,m+1) A[i][j]-=A[i-1][j];
        // 3 次作差
		int f=1;
		rep(i,1,n+1) rep(j,1,m+1) if(!V[i][j]) {
			ll s=0;
            // 子问题判定
			for(int a=i;a<=n+1;a+=k) for(int b=j;b<=m+1;b+=k) {
				V[a][b]=1;
				s+=A[a][b];
			}
			f&=s==0;
		}
		puts(f?"Yes":"No");
	}
}

[Codechef March Challenge 2021 Random Walk Queries(RWALKS) (动态点分治) ](https://www.codechef.com/MARCH21B/problems/RWALKS)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

Codechef March Challenge 2021 Random Walk Queries(RWALKS) (动态点分治)

题目大意：

对于给定的无根树 \(T\) ，要求强制在线维护两种操作

1.游走 \((u,d)\) ，以 \(u\) 为根在树上游走，从 \(u\) 开始，最多走 \(d\) 步，每次随机从儿子中选择一个点

2.查询 \(u\) ，当前 \(u\) 被遍历的期望次数

\[ \ \]

~~灵光一闪想到这么个憨批树上结构~~

对于更新 \((u,d)\) ，考虑 \(u\) 跨过当前点分根 到达其他点分子树里的贡献

一个点由当前点分根到达的概率是一个定值，可以预处理出来，并在查询时计算

因此更新贡献时，可以描述为 \(dep\leq d\) 的点接受到以 \(x\) 的概率访问当前点分根

可以简单用树状数组维护

为了剔除对于自己所在子树的非法贡献，需要额外开一些树状数组来维护

一个节点有 \(\log n\) 个点分父节点，每次需要两次树状数组查询

因此查询部分复杂度为 \(O(m\log ^2n)\) ，预处理以及空间复杂度为 \(O(n\log n)\)

const int N=2e5+10,K=19,P=1e9+7;

int n,m,I[N];
struct Edge{
	int to,nxt;
}e[N<<1];
int head[N],ecnt,deg[N];
void AddEdge(int u,int v){
	e[++ecnt]=(Edge){v,head[u]};
	head[u]=ecnt,deg[v]++;
}
#define erep(u) for(int i=head[u],v=e[i].to;i;i=e[i].nxt,v=e[i].to)

struct BIT{
	int n;
	vector <int> s;
	BIT(){};
	BIT(int n):n(n){ s.resize(n+1); }
	void Add(int p,int x){ 
		for(cmin(p,n);p;p-=p&-p) s[p]+=x,Mod1(s[p]);
	}
	int Que(int p){
		int res=0;
		while(p<=n) res+=s[p],Mod1(res),p+=p&-p;
		return res;
	}
} T[N];
vector <BIT> G[N];
//  Dep:点分树上的dep，id:节点在每层的编号， dep:节点在每层的dep，s:节点在每层由根到达的系数
int Dep[N],id[K][N],dep[K][N],s[K][N],vis[N],sz[N],fa[N],Root;

int mi,rt;
void FindRt(int n,int u,int f){
	int ma=0; sz[u]=1;
	erep(u) if(v!=f && !vis[v]) {
		FindRt(n,v,u);
		sz[u]+=sz[v],cmax(ma,sz[v]);
	}
	cmax(ma,n-sz[u]);
	if(mi>ma) mi=ma,rt=u;
}

int D,maxd;
void dfs(int u,int f,int id){
	cmax(maxd,dep[D][u]=dep[D][f]+1),::id[D][u]=id;
	erep(u) if(v!=f && !vis[v]) {
		s[D][v]=1ll*s[D][u]*I[deg[u]-1]%P;
		dfs(v,u,id);
	}
}

// 预处理点分治，开树状数组
int Divide(int n,int u){
	mi=1e9,FindRt(n,u,0),u=rt;
	int sonc=0;
	vis[u]=s[Dep[u]=D][u]=1,id[D][u]=-1;
	int t=0;
	erep(u) if(!vis[v]) {
		maxd=0;
		s[D][v]=1,dfs(v,u,sonc);
		G[u].pb(BIT(maxd));
		sonc++;
		cmax(t,maxd);
	}
	T[u]=BIT(t);
	erep(u) if(!vis[v]) {
		if(sz[v]>sz[u]) sz[v]=n-sz[u];
		D++,fa[Divide(sz[v],v)]=u,D--;
	}
	return u;
}

int sum[N];
int Que(int u){
	ll ans=sum[u];
	for(int v=u,d=Dep[v];(d--,v=fa[v]);) 
		ans=(ans+ 1ll* (T[v].Que(dep[d][u])+G[v][id[d][u]].Que(dep[d][u])) *s[d][u])%P;
	return (ans%P+P)%P;
}
void Upd(int u,int d){
	sum[u]++,Mod1(sum[u]),T[u].Add(d,I[deg[u]]);
	for(int v=fa[u],D=Dep[u]-1;v;v=fa[v],D--) {
		if(d<dep[D][u]) continue;
		int x=1ll*I[deg[u]]*s[D][u]%P;
		sum[v]+=x,Mod1(sum[v]);
		x=1ll*x*I[deg[v]-1]%P;
		T[v].Add(d-dep[D][u],x),G[v][id[D][u]].Add(d-dep[D][u],P-x);
	}
}

int lst;
int Get() { return (rd()+lst)%n+1; }

int main(){
	I[0]=I[1]=1;
	rep(i,2,N-1) I[i]=1ll*(P-P/i)*I[P%i]%P;
	n=rd(),m=rd();
	rep(i,2,n){
		int u=rd(),v=rd();
		AddEdge(u,v),AddEdge(v,u);
	}
	Root=Divide(n,1);
	while(m--) {
		int opt=rd();
		if(opt==1) {
			int u=Get(),d=Get();
			Upd(u,d);
		} else printf("%d\n",lst=Que(Get()));
	}
}

CodeChef 2020 November Challenge - Red-Black Boolean Expression

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

CodeChef 2020 November Challenge - Red-Black Boolean Expression

吐槽：这题很蠢，很套路

题目大意：

给定 \(n\) 个布尔变量 \(x_i\) ，每个变量有其反变量 \(\overline {x_i}\)

有 \(n\) 组关系 \(a_i,b_i\) ，要求 \(a_i\lor b_i\) 为真

并且保证所有 \(a_i,b_i\) 关系构成一张二分图，其中 \(x_i\) 与 \(\overline{x_i}\) 有一条边相连

给定每个变量的初始值 \(s_i\) ，以及翻转其所需的代价 \(C_i\) ，求最小满足条件的代价

\[ \ \]

\(a_i\lor b_i\) 为真即不存在 \(a_i,b_i\) 均为假的情况

如果是2-sat上的理解，即可以由 \(a_i\) 假推 \(b_i\) 真， \(b_i\) 假推 \(a_i\) 真，但是 \(2-sat\) 没法带权

由于题目保证了关系的二分图性质，不妨把所有变量分成两个集合 \(A,B\)

这个问题令人联想到网络流最小割模型，我们用一条边 \((u,v)\) 限制 \((u,v)\) 不同时为假的情况

对于 \(A\) 中的点，我们令源点 \(S\) 向 \(u\) 连的边 \((S,u,w)\) 表示 \(u\) 变成 \(0\) 所需代价，令 \((u,T,w)\) 表示 \(u\) 变成1的代价

对于 \(B\) 中的点，采取相反的连接方式

任意一个关系的两点不在同一集合，不妨对于 \(u\in A\) 的情况考虑，实际上可以分为两类考虑

\((u,v)\) 不同时为0，那么连接一条边 \((v,u,\infty)\) ，表示如果合法必然有一条让 \(u\) 或 \(v\) 变成1的边被割掉
\((u,v)\) 不同时为1，连接一条边 \((v,u,\infty)\) ，原理类似

然后就可以跑网络流最小割了

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

typedef long long ll;
typedef double db;
typedef unsigned long long ull;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
template <class T=int> T rd(){
	T s=0; int f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=1e5+10,INF=1e9+10;

int n,m,S,T;
char s[N];
struct Edge{
	int to,nxt,w;
}e[N<<1];
int head[N],ecnt=1;
void AddEdge(int u,int v,int w) {
	e[++ecnt]=(Edge){v,head[u],w};
	head[u]=ecnt;
}
void Link(int u,int v,int w){ 
	AddEdge(u,v,w),AddEdge(v,u,0);
}

int F[N],X[N],Y[N],col[N],A[N];
// 这里我用带权并查集实现了二分图
int Find(int x){ 
	if(F[x]==x) return x;
	int f=F[x]; F[x]=Find(F[x]);
	col[x]^=col[f];
	return F[x];
}

int dis[N];
int Bfs() {
	rep(i,1,T) dis[i]=INF;
	static queue <int> que;
	dis[S]=0; que.push(S);
	while(!que.empty()) {
		int u=que.front(); que.pop();
		for(int i=head[u];i;i=e[i].nxt) {
			int v=e[i].to,w=e[i].w;
			if(!w || dis[v]<=dis[u]+1) continue;
			dis[v]=dis[u]+1,que.push(v);
		}
	}
	return dis[T]<INF;
}

int Dfs(int u,int in) {
	if(u==T) return in;
	int out=0;
	for(int i=head[u];i;i=e[i].nxt) {
		int v=e[i].to,w=e[i].w;
		if(!w || dis[v]!=dis[u]+1) continue;
		int t=Dfs(v,min(in-out,w));
		e[i].w-=t,e[i^1].w+=t,out+=t;
		if(in==out) break;
	}
	if(!out) dis[u]=0;
	return out;
}

int Dinic() {
	int ans=0; 
	while(Bfs()) ans+=Dfs(S,INF);
	return ans;
}


int main() {
	rep(kase,1,rd()) {
		n=rd(),m=rd();
		rep(i,1,n) F[i]=i,col[i]=0;
		scanf("%s",s+1);
		S=n+1,T=n+2;
		rep(i,1,n) A[i]=rd();
		rep(i,1,m) {
			X[i]=rd(),Y[i]=rd();
			int x=abs(X[i]),y=abs(Y[i]);
			int u=Find(x),v=Find(y);
			if(u==v) continue;
			F[u]=v,col[u]=col[x]^col[y]^(X[i]<0)^(Y[i]<0)^1;
		}
		rep(i,1,n) Find(i);
		rep(i,1,n) {
			if(col[i]^s[i]^'0') Link(S,i,A[i]);
			else Link(i,T,A[i]);
		}
		rep(i,1,m) {
			int t=col[abs(X[i])]^(X[i]<0);
			assert(col[abs(X[i])]^col[abs(Y[i])]^(X[i]<0)^(Y[i]<0));
			if(t) Link(abs(X[i]),abs(Y[i]),INF);
			else Link(abs(Y[i]),abs(X[i]),INF);
		}
		printf("%d\n",Dinic());
		rep(i,1,T) head[i]=0; ecnt=1;
	}
}

CodeChef November Challenge2019 Winning Ways (3-FWT)

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

CodeChef November Challenge2019 Winning Ways (3-FWT)

显然每个把每个数换成其因子个数-1，就能转为一个扩展的 \(\text{Nim}\) 游戏

每次操作 \(1,2,\cdots,k\) 堆的 \(\text{Nim}\) 游戏，其判定方法是：

将每个数二进制分解，对于每个二进制位上分别计算1的个数 \(\mod (k+1)\) ，如果均为0则先手必败

对于这道题 \(k=3\) ，我们考虑将其转为二进制之后的形式累成3进制，然后就能进行3进制按位不进位加法，即类异或

然后问题实际上有非常多的部分需要考虑

Part1 如何求因子个数

一个简单的方法是枚举 \([1,\sqrt n]\) 判断是否整除，复杂度过高

对于 \(n=\prod p_i^{c_i}\) ( \(p_i\) 为质数)，其因子个数为 \(\prod (c_i+1)\)

由这个式子对于 \(n\) 进行质因数分解，枚举 \([1,\sqrt n]\) 中的质数，复杂度为 \(O(\pi(\sqrt n))=O(\frac{\sqrt n}{\log n})\) ，这个应该够了？

然后是一个常规套路型的分解方法：

先对于 \([1,n^{\frac{1}{3}}]\) 的质数筛 \(n\) ，剩余的部分只有3种情况

\(n\) 被筛成1了
\(n\) 被筛到只剩一个质数，可以用 \(\text{Miller_Rabin}\) 算法快速判断，可以参考
\(n\) 仍然是若干质数的乘积，此时质因子必然 \(>n^{\frac{1}{3}}\) ，因此最多只有两个

那么只需要判断 \(n\) 是否是完全平方数即可

总复杂度为 \(O(w\cdot \log n+\frac{n^{\frac{1}{3}}}{\log n})\) ，其中 \(w\) 为 \(\text{Miller_Rabin}\) 筛选次数

const int N=2e5+10;
int notpri[N],pri[N],pc;

ll qpow(ll x,ll k=P-2,int P=::P) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}
int Trans(int x){ 
    static int buf[20],l=0;
    l=0;
    while(x) buf[++l]=x%2,x/=2;
    int s=0;
    drep(i,l,1) s=s*3+buf[i];
    cmax(ma,s);
    return s;
}

int Miller_Rabin(int x) {
    if(x<N) return !notpri[x];
    if(x==2) return 1;
    if(x<=1 || ~x&1) return 0;
    ll s=0,t=x-1;
    while(~t&1) s++,t>>=1;
    rep(i,1,4) {
        ll a=pri[rand()%pc+1],b=qpow(a,t,x),c;
        rep(j,1,s) {
            c=1ll*b*b%x;
            if(c==1 && b!=1 && b!=x-1) return 0;
            b=c;
        }
        if(b!=1) return 0;
    }
    return 1;
}
int CheckSqr(int n){
    int y=round(sqrt(n));
    return y*y==n;
}

int Count(int n){
    int res=1;
    for(int i=1;i<=pc && pri[i]*pri[i]*pri[i]<=n;++i) if(n%pri[i]==0) {
        int c=0;
        while(n%pri[i]==0) n/=pri[i],c++;
        res*=c+1;
    }
    if(n==1) return res;
    if(CheckSqr(n)) return res*3;
    if(Miller_Rabin(n)) return res*2;
    return res*4;
}

void Init(){
    rep(i,2,N-1) if(!notpri[i]) {
        pri[++pc]=i;
        for(int j=i+i;j<N;j+=i) notpri[j]=1;
    }
}

Part2 快速计算答案

\(10^9\) 以内的数，最大因子个数为 \(1334\) ，这个数为 \(931170240\)

转成二进制之后最多包含 \(11\) 位，三进制下最大为 \(3^{11}-1=177146\) ，令这个上界为 \(M\)

一种非常暴力的方法就是直接枚举， \(NM\) 计算每次选择一个数，复杂度为 \(O(NMK)\) ，应该可以通过 \(N\leq 12\) 的数据

一个比较浅显的优化可以用快速幂维护乘法，复杂度为 \(O(M^2\log K)\)

由于是3进制类异或，接下来考虑用 \(\text{3-FWT}\) 优化乘法，可以参考

模数为 \(P=10^9+7\) ，不存在整数 \(3\) 阶单位根，因此要用类似拆系数 \(\text{FFT}\) 方法做

复杂度为 \(O(M\log M\log K)\) ，似乎已经比较小了，但是常数非常大，应该难以通过

int R;// R为上界
struct Cp{
    db x,y;
    Cp(){}
    Cp(db x,db y):x(x),y(y) {}
    Cp operator + (const Cp t){ return Cp(x+t.x,y+t.y); }
    Cp operator - (const Cp t){ return Cp(x-t.x,y-t.y); }
    Cp operator * (const Cp t){ return Cp(x*t.x-y*t.y,x*t.y+y*t.x); }
} A[N],B[N],C[N],D[N];
Cp w[30];

int Add(int x,int y) {
    static int A[20],B[20];
    rep(i,0,19) A[i]=x%3,x/=3;
    rep(i,0,19) B[i]=y%3,y/=3;
    int ans=0;
    drep(i,19,0) ans=ans*3+((A[i]+B[i])%3);
    return ans;
}

void FWT(Cp *a,int f) {
    for(int i=1;i<R;i*=3) {
        for(int l=0;l<R;l+=i*3) {
            for(int j=l;j<l+i;++j) {
                static Cp t[3];
                if(f==1) {
                    t[0]=a[j]+a[j+i]+a[j+i*2];
                    t[1]=a[j]+w[1]*a[j+i]+w[2]*a[j+i*2];
                    t[2]=a[j]+w[2]*a[j+i]+w[1]*a[j+i*2];
                } else {
                    t[0]=a[j]+a[j+i]+a[j+i*2];
                    t[1]=a[j]+w[2]*a[j+i]+w[1]*a[j+i*2];
                    t[2]=a[j]+w[1]*a[j+i]+w[2]*a[j+i*2];
                }
                a[j]=t[0],a[j+i]=t[1],a[j+i*2]=t[2];
                if(f==-1) {
                    rep(d,0,2) {
                        a[j+i*d].x/=3,a[j+i*d].y/=3;
                    }
                }
            }
        }
    }
}

const int S=(1<<15)-1;

#define FWTs

struct Poly{ 
    int a[N];
    Poly operator * (const Poly __) const {
        Poly res;
        // 拆系数，任意模数3-FWT
#ifdef FWTs
        rep(i,0,R-1) A[i]=Cp((a[i]&S),(a[i]>>15)),B[i]=Cp((a[i]&S),0);
        rep(i,0,R-1) C[i]=Cp((__.a[i]&S),(__.a[i]>>15));
        FWT(A,1),FWT(B,1),FWT(C,1);
#define E(x) ((ll)(x+0.5))%P
        rep(i,0,R-1) A[i]=A[i]*C[i];
        rep(i,0,R-1) B[i]=B[i]*C[i];
        FWT(A,-1),FWT(B,-1);
        rep(i,0,R-1) {
            ll a=E(B[i].x),b=E(A[i].y),c=E(B[i].x-A[i].x);
            res.a[i]=(a+1ll*b*(S+1)+1ll*c*(S+1)*(S+1))%P;
        }
#else
        rep(i,0,R-1) res.a[i]=0;
        rep(i,0,R-1) if(a[i]) rep(j,0,R-1) if(__.a[j]) {
            int k=Add(i,j);
            res.a[k]=(res.a[k]+1ll*a[i]*__.a[j])%P;
        }
#endif
        return res;
    }
} x,res;



ll qpow(ll x,ll k=P-2) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}

int main(){
    rep(i,2,N-1) if(!notpri[i]) {
        pri[++pc]=i;
        for(int j=i+i;j<N;j+=i) notpri[j]=1;
    }
    w[0]=Cp(1,0);
    rep(i,1,29) w[i]=w[i-1]*Cp(cos(2*Pi/3),sin(2*Pi/3));
    // 复平面3阶单位根

    n=rd(),m=rd();
    rep(i,1,n) x.a[Count(rd())]++;
    for(R=1;R<=ma;R*=3);

    res.a[0]++;
    for(;m;m>>=1) {
        if(m&1) res=res*x;
        x=x*x;
    }
    int ans=0;
    rep(i,1,R-1) ans+=res.a[i],Mod1(ans);
    printf("%d\n",ans);
}

Further

对于形式幂级数多项式，我们知道 \(K\) 次幂的循环卷积可以直接 \(\text{DFT}\) 一次，每一位快速幂，然后 \(\text{IDFT}\)

同理的，如果你学习了 \(\text{K-FWT}\) 就知道这就是一个按 \(K\) 进制位，每一位分别进行循环卷积，因此也可以用类似的方法做

但是遇到一个非常大的问题就是无法找到模意义下的 \(3\) 阶单位根(指 \(3\not |(P-1)\) )

如果用复平面单位根 \(\omega_n=cos(\frac{2\pi}{n})+sin(\frac{2\pi}{n})\cdot i\) ( \(i=\sqrt {-1})\) ，无法在计算时保证值域精度

这里由于 \(n=3\) 比较特殊，发现 \(\omega_3=cos(\frac{2\pi}{3})+sin(\frac{2\pi}{3})\cdot i=-\frac{1}{2}+\frac{\sqrt 3}{2}\cdot i\)

而 \(3\) 在 \(\mod 10^9+7\) 下存在二次剩余，因此可以用一个模意义下的复数描述复平面单位根

应该是有通行的单位根求法，会根据 \(n\) 不同要用更复杂的高维复数描述，但是我并不会.jpg

总复杂度为 \(O(M(\log M+\log K))\) ，分别为进行 \(\text{3-FWT}\) 以及快速幂的复杂度

Code总览:

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

typedef long long ll;
typedef unsigned long long ull;
typedef long double db;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
int rd(){
    int s=0,f=0;
    while(!isdigit(IO=getchar())) if(IO=='-') f=1;
    do s=(s<<1)+(s<<3)+(IO^'0');
    while(isdigit(IO=getchar()));
    return f?-s:s;
}
bool Mbe;
const int N=2e5,P=1e9+7;
const int Quad3=82062379; // 3在Mod P意义下的二次剩余
const db Pi=acos((db)-1);

const int MaxX=931170240;


int n,m,ma,R;
int notpri[N],pri[N],pc;

ll qpow(ll x,ll k=P-2,int P=::P) {
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}
int Trans(int x){ 
    static int buf[20],l=0;
    l=0;
    while(x) buf[++l]=x%2,x/=2;
    int s=0;
    drep(i,l,1) s=s*3+buf[i];
    cmax(ma,s);
    return s;
}

int Miller_Rabin(int x) {
    if(x<N) return !notpri[x];
    if(x==2) return 1;
    if(x<=1 || ~x&1) return 0;
    ll s=0,t=x-1;
    while(~t&1) s++,t>>=1;
    rep(i,1,4) {
        ll a=pri[rand()%pc+1],b=qpow(a,t,x),c;
        rep(j,1,s) {
            c=1ll*b*b%x;
            if(c==1 && b!=1 && b!=x-1) return 0;
            b=c;
        }
        if(b!=1) return 0;
    }
    return 1;
}
int CheckSqr(int n){
    int y=round(sqrt(n));
    return y*y==n;
}

int Count(int n){
    int res=1;
    for(int i=1;i<=pc && pri[i]*pri[i]*pri[i]<=n;++i) if(n%pri[i]==0) {
        int c=0;
        while(n%pri[i]==0) n/=pri[i],c++;
        res*=c+1;
    }
    if(n==1) return res;
    if(CheckSqr(n)) return res*3;
    if(Miller_Rabin(n)) return res*2;
    return res*4;
}

struct Cp{
    int x,y;
    Cp(){}
    Cp(int x,int y):x(x),y(y) {}
    Cp operator + (const Cp t){ 
        Cp res(x+t.x,y+t.y);
        Mod1(res.x),Mod1(res.y);
        return res;
    }
    Cp operator * (const Cp t){ return Cp((1ll*x*t.x+1ll*(P-y)*t.y)%P,(1ll*x*t.y+1ll*y*t.x)%P); }
} A[N]; // 模意义下 模拟复平面单位根
Cp w1,w2; // 3阶单位根及其平方

Cp qpow(Cp x,ll k=P-2) {
    Cp res(1,0);
    for(;k;k>>=1,x=x*x) if(k&1) res=res*x;
    return res;
}

// 下面是展开的FWT式子
void FWT() {
    for(int i=1;i<R;i*=3) {
        for(int l=0;l<R;l+=i*3) {
            for(int j=l;j<l+i;++j) {
                Cp a=A[j]+A[j+i]+A[j+i*2];
                Cp b=A[j]+w1*A[j+i]+w2*A[j+i*2];
                Cp c=A[j]+w2*A[j+i]+w1*A[j+i*2];
                A[j]=a,A[j+i]=b,A[j+i*2]=c;
            }
        }
    }
}

void IFWT() {
    for(int i=1;i<R;i*=3) {
        for(int l=0;l<R;l+=i*3) {
            for(int j=l;j<l+i;++j) {
                Cp a=A[j]+A[j+i]+A[j+i*2];
                Cp b=A[j]+w2*A[j+i]+w1*A[j+i*2];
                Cp c=A[j]+w1*A[j+i]+w2*A[j+i*2];
                A[j]=a,A[j+i]=b,A[j+i*2]=c;
            }
        }
    }
    ll base=qpow(R);
    rep(i,0,R-1) A[i].x=A[i].x*base%P;
}

int main(){
    rep(i,2,N-1) if(!notpri[i]) {
        pri[++pc]=i;
        for(int j=i+i;j<N;j+=i) notpri[j]=1;
    }
    w1=Cp(P-(P+1)/2,1ll*Quad3*(P+1)/2%P);
    w2=w1*w1;

    rep(kase,1,rd()) {
        n=rd(),m=rd();
        rep(i,1,n) A[Trans(Count(rd())-1)].x++;
        for(R=1;R<=ma;R*=3);
        FWT();
        rep(i,0,R-1) A[i]=qpow(A[i],m);
        IFWT();
        int ans=0;
        rep(i,1,R-1) ans+=A[i].x,Mod1(ans);
        rep(i,0,R-1) A[i].x=A[i].y=0;
        printf("%d\n",ans);
    }
}

2014多校6 Another Letter Tree

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

2014多校6 Another Letter Tree

点分治做法

就裸地离个线，放到点分治上，从每个根开始，维护 \(dp_{u,l,r}\) 表示这条链匹配了序列中 \([l,r]\) 的部分

注意dp数组要一正一反，俩家伙一个含根一个不含

查询要合并两个dp数组，但是只需要知道 \(dp_{1,|s_0|}\) ，因此合并复杂度是 \(O(|s_0|)\) 的

最终复杂度，处理为 \(O(n\log n|s_0|^2+q|s_0|)\)

树剖线段树做法

类似上面的dp，线段树维护即可

问题1

需要存正反！

然后你发现内存从中间裂开！！

正反分两次，离线跑两次就可以了啊啊啊啊

问题2

如果直接查询合并，合并两个dp数组复杂度为 \(O(|s_0|^3)\)

查询复杂度为 \(O(q\log ^2n|s_0|^3)\)

妙啊!!比 \(n^2\) 还大!!

所以最后不能合并dp数组，而应该直接累加到答案数组上

问题3

没错现在我们的复杂度为 \(O(q\log ^2n|s_0|^2)\)

依然大得令人无法忍受

~~但是没想到吧，数据全部都是链，树剖是 \(O(1)\) 的~~

优化：查询重链时，只有最后依次是在链上查询 \([l,r]\) 都在中间的，而对于直接跳到top的部分，可以预处理出来

算上线段树的预处理，这样总复杂度就是 \(O(n|s_0|^3+q\log n |s_0|^2)\)

并查集做法

把问题拆成查询两条 \(u\) 到它的祖先 \(v\) 的答案

每个节点存储一个dp矩阵，用带权并查集维护

具体方法是:将询问按照 \(\text{LCA}\) 深度逆序排序后，每次查询一直将 \(u\) 合并到 \(v\) 为止

复杂度为 \(O(n\alpha(n)|s_0|^3)\) ，理论上来说，这个转移矩阵应当很稀疏，乘法应该很快，但是实际常数比较大

#include<bits/stdc++.h>
using namespace std;
#pragma GCC optimize(2)
#define reg register
#define Mod1(x) ((x>=P)&&(x-=P))
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)
char IO;
int rd(){
	int s=0;
	while(!isdigit(IO=getchar()));
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return s;
}

const int N=5e4+10,P=10007;

int n,m,len;
char str[N],t[N];
struct Node{
	int a[31][31];
	void clear(){ memset(a,0,sizeof a); }
	void Init(int u) { clear(); rep(i,1,len) if(t[i]==str[u]) a[i][i]=1; }
	Node operator + (const Node &x) const {
		Node res; res.clear();
		rep(i,1,len) for(int j=i;j<len && a[i][j];++j) for(int k=j+1;k<=len && x.a[j+1][k];++k) res.a[i][k]=(res.a[i][k]+a[i][j]*x.a[j+1][k])%P;
		rep(i,1,len) rep(j,i,len) {
			res.a[i][j]+=a[i][j],Mod1(res.a[i][j]);
			res.a[i][j]+=x.a[i][j],Mod1(res.a[i][j]);
		}
		return res;
	}
} s[N];

struct Edge{
	int to,nxt;
}e[N<<1];
int head[N],ecnt;
void AddEdge(int u,int v) {
	e[++ecnt]=(Edge){v,head[u]};
	head[u]=ecnt; 
}
#define erep(u,i) for(int i=head[u];i;i=e[i].nxt)

int QX[N],QY[N],QL[N];
int dep[N],id[N],fa[N][18];
void pre_dfs(int u,int f) {
	dep[u]=dep[fa[u][0]=f]+1;
	rep(i,1,17) fa[u][i]=fa[fa[u][i-1]][i-1];
	erep(u,i) {
		int v=e[i].to;
		if(v==f) continue;
		pre_dfs(v,u);
	}
}
int LCA(int x,int y) {
	if(dep[x]<dep[y]) swap(x,y);
	for(int i=0,del=dep[x]-dep[y];(1<<i)<=del;++i) if(del&(1<<i)) x=fa[x][i];
	if(x==y) return x;
	drep(i,17,0) if(fa[x][i]!=fa[y][i]) x=fa[x][i],y=fa[y][i];
	return fa[x][0];
}

int F[N],ans[N][31],Ans[N];
int Find(int x,int k=0) {
	if(F[x]==x) return x;
	int f=F[x]; F[x]=Find(f,k);
	if(F[f]!=f){
		if(!k) s[x]=s[x]+s[f];
		else s[x]=s[f]+s[x];
	}
	return F[x];
}

int main(){
	rep(kase,1,rd()) {
		n=rd(),m=rd();
		rep(i,1,n) head[i]=ecnt=0;
		rep(i,2,n) {
			int u=rd(),v=rd();
			AddEdge(u,v),AddEdge(v,u);
		}
		scanf("%s%s",str+1,t+1),len=strlen(t+1);
		pre_dfs(1,0);
		rep(i,1,m) id[i]=i,QX[i]=rd(),QY[i]=rd(),QL[i]=LCA(QX[i],QY[i]);
		sort(id+1,id+m+1,[&](int x,int y){ return dep[QL[x]]>dep[QL[y]]; });

		rep(i,1,n) F[i]=i,s[i].clear();
		rep(k,1,m){
			int i=id[k],x=QX[i];
			while(1){
				int y=Find(x);
				if(y==QL[i]) break;
				F[y]=fa[y][0],s[y].Init(y);
			}
			ans[i][0]=1;
			rep(j,1,len) ans[i][j]=s[x].a[1][j];
			drep(j,len,1) if(str[QL[i]]==t[j]) ans[i][j]+=ans[i][j-1],Mod1(ans[i][j]);
		}

		rep(i,1,n) F[i]=i,s[i].clear();
		rep(k,1,m){
			int i=id[k],x=QY[i]; Ans[i]=0;
			while(1){
				int y=Find(x,1);
				if(y==QL[i]) break;
				F[y]=fa[y][0],s[y].Init(y);
			}
			rep(j,0,len-1) Ans[i]=(Ans[i]+ans[i][j]*s[x].a[j+1][len])%P;
			Ans[i]=(Ans[i]+ans[i][len])%P;
		}
		rep(i,1,m) printf("%d\n",Ans[i]);
	}
}

\[ \ \]

伪矩阵求逆做法

同样的，把问题拆成查询两条 \(u\) 到它的祖先 \(v\) 的答案(不包含v)

以从 \(v\) 到 \(u\) 的字符串为例，设 \(dp_u\) 为 \(u\) 的祖先链的dp矩阵，我们要求的部分答案是 \(x\)

则 \(dp_v\cdot x=dp_u, x=\frac{dp_u}{dp_v}\)

一般来说，矩阵求逆是一个很难实现的东西

但是发现对于一种 \(dp\) ，它的矩阵一定是一个上/下对角的矩阵

我们需要求出矩阵第一维为1或者第二维为 \(|s_0|\) 的部分

如果暴力求，可以看做求解一个 \(|s_0|\) 元的线性方程组，可以用高斯消元在 \(O(|s_0|^3)\) 时间内求解

而实际上，这个线性方程组是含拓扑序关系的，任何一个含拓扑关系的线性方程组求解是不需要高斯消元的

而且这个问题列出的方程矩阵就已经是上对角矩阵了

所以写出来就是容斥吧

tips: 预处理部分一次只插入一个字符，复杂度为 \(O(n|s_0|^2)\) (也可以认为是稀疏矩阵乘法)

查询部分求解线性方程复杂度为 \(O(|s_0|^2)\) ，合并答案复杂度为 \(O(|s_0|)\)

因此复杂度为 \(O((n+q)|s_0|^2)\)

Code: 注意两种dp共用了一个数组

#include<bits/stdc++.h>
using namespace std;
#define reg register
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)
char IO;
int rd(){
	int s=0;
	while(!isdigit(IO=getchar()));
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return s;
}

const int N=5e4+10,P=10007;

int n,m,len;
char str[N],t[N];
struct Edge{
	int to,nxt;
}e[N<<1];
int head[N],ecnt;
void AddEdge(int u,int v) {
	e[++ecnt]=(Edge){v,head[u]};
	head[u]=ecnt; 
}
#define erep(u,i) for(int i=head[u];i;i=e[i].nxt)

int dep[N],id[N],fa[N][18];
int dp[N][31][31];
int f1[31],f2[31];

void pre_dfs(int u,int f) {
	dep[u]=dep[fa[u][0]=f]+1;
	rep(i,1,17) fa[u][i]=fa[fa[u][i-1]][i-1];
	rep(i,1,len) rep(j,1,len) {
		dp[u][i][j]=dp[f][i][j];
		if(str[u]==t[j]) {
			if(i==j) dp[u][i][j]++,Mod1(dp[u][i][j]);
			if(i<j) dp[u][i][j]+=dp[f][i][j-1],Mod1(dp[u][i][j]);
			if(i>j) dp[u][i][j]+=dp[f][i][j+1],Mod1(dp[u][i][j]);
		}
	}
	erep(u,i) {
		int v=e[i].to;
		if(v==f) continue;
		pre_dfs(v,u);
	}
}
int LCA(int x,int y) {
	if(dep[x]<dep[y]) swap(x,y);
	for(int i=0,d=dep[x]-dep[y];(1<<i)<=d;++i) if(d&(1<<i)) x=fa[x][i];
	if(x==y) return x;
	drep(i,17,0) if(fa[x][i]!=fa[y][i]) x=fa[x][i],y=fa[y][i];
	return fa[x][0];
}

void Calcdp1(int u,int f) {
	rep(i,f1[0]=1,len) {
		f1[i]=dp[u][i][1];
		rep(j,0,i-1) f1[i]=(f1[i]-f1[j]*dp[f][i][j+1])%P;
		Mod2(f1[i]);
	}
}

void Calcdp2(int u,int f) {
	drep(i,len,f2[len+1]=1) {
		f2[i]=dp[u][i][len];
		rep(j,i+1,len+1) f2[i]=(f2[i]-dp[f][i][j-1]*f2[j])%P;
		Mod2(f2[i]);
	}
}
int Que(int x,int y) {
	int lca=LCA(x,y);
	Calcdp1(x,fa[lca][0]),Calcdp2(y,lca);
	int ans=0;
	rep(i,0,len) ans=(ans+f1[i]*f2[i+1])%P;
	return ans;
}

int main(){
	rep(kase,1,rd()) {
		n=rd(),m=rd();
		rep(i,1,n) head[i]=ecnt=0;
		rep(i,2,n) {
			int u=rd(),v=rd();
			AddEdge(u,v),AddEdge(v,u);
		}
		scanf("%s%s",str+1,t+1),len=strlen(t+1);
		pre_dfs(1,0);
		rep(i,1,m) {
			int x=rd(),y=rd();
			printf("%d\n",Que(x,y));
		}
	}
}

HDU-5608

Posted on 2023-08-12 Edited on 2024-04-19 In CP 题解

HDU-5608

题意： \(G(n)=n^2−3n+2=\sum_{d|n}F(d)\) ，求 \(\sum_1^nF(i)\)

反演得到： \(F(n)=\sum_{d|n}\mu(d)G(\frac{n}{d})\)

则 \(\sum_1^nF(i)=\sum_i\sum_{d|i}\mu(d)G(\frac{i}{d})\)

\(\sum_1^nF(i)=\sum_{i=1}^{n}G(i)\sum_{j=1}^{\lfloor \frac{n}{i}\rfloor }\mu(j)\)

问题就是要快速求 \(G(i)\) 前缀和和 \(\mu(i)\) 前缀和

第一个是 \(O(1)\) 求，第二个是杜教筛


const int N=5e6+10,P=1e9+7;

ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

const int Inv6=qpow(6);

int n;
char notpri[N],w[N];
int pri[N/4],pc,Sw[N];
map <int,int> M;

int SumG(int n){ // O(1)求出G函数的前缀和
	ll ans=1ll*n*(n+1)%P*(2*n+1)%P*Inv6%P;
	ans=(ans-3ll*n*(n+1)/2%P+2*n)%P;
	ans=(ans%P+P)%P;
	return ans;
}

int Sumw(int n){ // 杜教筛求mobius函数前缀和
	if(n<N) return Sw[n];
	if(M.count(n)) return M[n];
	int ans=1;
	for(int i=2,j;i<=n;i=j+1) {
		j=n/(n/i);
		ans-=(j-i+1)*Sumw(n/i);
	}

	return M[n]=ans;
}

int SumF(int n){ // 答案
	int ans=0;
	for(int i=1,j;i<=n;i=j+1) {
		j=n/(n/i);
		ans=(ans+1ll*(SumG(j)-SumG(i-1))*Sumw(n/i))%P;
	}
	ans=(ans%P+P)%P;
	return ans;
}

int main(){
	w[1]=1;
	rep(i,2,N-1) {
		if(!notpri[i]) pri[++pc]=i,w[i]=-1;
		for(int j=1;j<=pc && 1ll*i*pri[j]<N;++j) {
			notpri[i*pri[j]]=1;
			if(i%pri[j]==0) {
				w[i*pri[j]]=0;
				break;
			}
			w[i*pri[j]]=-w[i];
		}
	}
	rep(i,1,N-1) Sw[i]=Sw[i-1]+w[i];
	rep(kase,1,rd()) printf("%d\n",SumF(rd()));
}