数字表格(product)

Portal -->broken qwq

Description

​　　 求(prodlimits_{i=1}^nprodlimits_{j=1}^m f[gcd(i,j)](mod 10^9+7))，其中(f[0]=0,f[1]=1,f[i]=f[i-1]+f[i-2])

​　　 数据范围：多组数据，数据组数(T<=1000,1<=n,m<=10^6)

Solution

​​　　 这题。。一开始我觉得这题好像直接求个范围内(gcd(i,j)=d)的对数就做完了qwq

​​　　 式子长这个样子

[egin{aligned} Ans&=prodlimits_{d=1}^{min(n,m)}f(d)^{s(d)}\ s(d)&=sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=d]\ &=sumlimits_{k=1}^{lfloorfrac{min(n,m)}{d} floor}mu(k)lfloorfrac{n}{kd} floorlfloorfrac{m}{kd} floor end{aligned} ]

​　　 然后发现单组(O(nsqrt n logn))然后数据组数(1000)==（不过还是有60）

​　　（然而实际上也是可以做的==详见后面的update）

​​　　

​​　　 正解并没有题面那么简单qwq实际上正解很妙qwq

​​　　 我们考虑构造一个函数(g(x))，满足：

[f[n]=prodlimits_{d|n}g(d) ]

​​　　 这样一来我们就可以把原来的式子化一下变成：

[egin{aligned} &prodlimits_{i=1}^nprodlimits_{j=1}^mf[gcd(i,j)]\ =&prodlimits_{i=1}^nprodlimits_{j=1}^mprod_{d|i,d|j}g(d)\ =&prodlimits_{d=1}^{min(n,m)}g(d)^{lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor} end{aligned} ]

​​　　 然后这样我们就不用再在外面枚举一次了，单组复杂度(O(sqrt nlogn))

​​　　 前提是我们知道(g(d))的前缀积

​　　 这个东西长得跟反演式子差不多，考虑反演式子的容斥解释，本质上就是容斥的加加减减，换成乘法的话就是乘除就好了，所以我们直接将(mu)搬到指数上面就好了

[g(n)=prodlimits_{d|n}f(d)^{mu(frac{n}{d})} ]

​　　 这个可以大力枚举约数什么的直接预处理

​　　 然后这题就做完啦ovo

​　　

​　　 重大update！！

​　　 这题其实用最上面一开始推的那个式子也是可以做的！只不过。。我们可以把指数的式子少化一步变成：

[s(d)=sumlimits_{i=1}^{lfloorfrac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}[gcd(i,j)=1] ]

​　　 然后我们会发现底数也是可以分块的！我们预处理一下(f)的前缀积，然后直接大力分块一下就好了（注意(f[0]=0)的情况在求逆元的时候返回(1)）

​　　 这种做法的话。。很暴力但是可以压线过==算是一种水法吧。。不过因为式子非常好推场上比较容易想到所以还是记录一下

​　　 复杂度的话给给全说因为和不预处理的杜教筛形式类似所以是。。(O(n^{frac{3}{4}}))的然后一定要算的话。。再加上一个(sqrt n logn)（快速幂）
​　　

​​　　 正解的代码大概长这个样子

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
const int N=1e6+10,MOD=1e9+7;
int f[N],g[N],p[N],miu[N],invf[N],invg[N];
int vis[N];
int n,m,T,ans;
int add(int x,int y){return (1LL*x+MOD+y)%MOD;}
int mul(int x,int y){return 1LL*x*y%MOD;}
int ksm(int x,ll y){
	int ret=1,base=x;
	for (;y;y>>=1,base=mul(base,base))
		if (y&1) ret=mul(ret,base);
	return ret;
}
int inv(int x){return ksm(x,MOD-2);}
void prework(int n){
	int cnt=0;
	miu[1]=1;
	for (int i=2;i<=n;++i){
		if (!vis[i])
			p[++cnt]=i,miu[i]=-1;
		for (int j=1;j<=cnt&&p[j]*i<=n;++j){
			vis[i*p[j]]=true;
			if (i%p[j]==0){
				miu[i*p[j]]=0;
				break;
			}
			else
				miu[i*p[j]]=-miu[i];
		}
	}
	f[0]=0; f[1]=1; g[1]=1;
	for (int i=2;i<=n;++i) f[i]=add(f[i-1],f[i-2]),g[i]=1;
	for (int i=1;i<=n;++i) invf[i]=inv(f[i]);
	for (int i=1;i<=n;++i)
		for (int j=i;j<=n;j+=i){
			if (miu[j/i]==-1) 
				g[j]=mul(g[j],invf[i]);
			else if (miu[j/i]==1)
				g[j]=mul(g[j],f[i]);
		}
	for (int i=2;i<=n;++i) g[i]=mul(g[i],g[i-1]);
	invg[0]=1;
	for (int i=1;i<=n;++i) invg[i]=inv(g[i]);
}
int solve(){
	int ret=1;
	if (n>m) swap(n,m);
	for (int i=1,pos=0;i<=n;i=pos+1){
		pos=min(n/(n/i),m/(m/i));
		ret=mul(ret,ksm(mul(g[pos],invg[i-1]),1LL*(n/i)*(m/i)));
	}
	return ret;
}

int main(){
#ifndef ONLINE_JUDGE
	freopen("a.in","r",stdin);
#endif
	scanf("%d",&T);
	prework(N-10);
	for (int o=1;o<=T;++o){
		scanf("%d%d",&n,&m);
		ans=solve();		
		printf("%d
",ans);
	}
}

​　　
​　　 水法的代码大概长这个样子

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
const int N=1e6+10,MOD=1e9+7,inf=2147483647;
ll f[N],miu[N],p[N];
int vis[N];
int n,m,T;
int mul(int x,int y){return 1LL*x*y%MOD;}
int add(int x,int y){return (1LL*x+y)%MOD;}
void prework(int n){
	int cnt=0;
	miu[1]=1;
	for (int i=2;i<=n;++i){
		if (!vis[i]){
			p[++cnt]=i; miu[i]=-1;
		}
		for (int j=1;j<=cnt&&i*p[j]<=n;++j){
			vis[i*p[j]]=1;
			if (i%p[j])
				miu[i*p[j]]=-miu[i];
			else{
				miu[i*p[j]]=0; break;
			}
		}
	}
	for (int i=2;i<=n;++i) miu[i]+=miu[i-1];
	f[0]=0; f[1]=1;
	for (int i=2;i<=n;++i) f[i]=add(f[i-1],f[i-2]);
	for (int i=2;i<=n;++i) f[i]=mul(f[i],f[i-1]);
}
int ksm(int x,ll y){
	int ret=1,base=x;
	for (;y;y>>=1,base=mul(base,base))
		if (y&1) ret=mul(ret,base);
	return ret;
}
int inv(int x){return x==0?1:ksm(x,MOD-2);}
ll s(int n,int m){
	ll ret=0;
	for (int i=1,pos=0;i<=n;i=pos+1){
		pos=min(n/(n/i),m/(m/i));
		//ret=add(ret,mul(add(miu[pos],-miu[i-1]),mul(n/i,m/i)));
		ret+=1LL*(n/i)*(m/i)*(miu[pos]-miu[i-1]);
	}
	return ret;
}
int solve(int n,int m){
	int ret=1;
	for (int i=1,pos=0;i<=n;i=pos+1){
		pos=min(n/(n/i),m/(m/i));
		ret=mul(ret,ksm(mul(f[pos],inv(f[i-1])),s(n/i,m/i)));
	}
	return ret;
}

int main(){
#ifndef ONLINE_JUDGE
	freopen("a.in","r",stdin);
#endif
	scanf("%d",&T);
	prework(N-10);
	for (int i=1;i<=T;++i){
		scanf("%d%d",&n,&m);
		if (n>m) swap(n,m);
		printf("%d
",solve(n,m));
	}
}