[bzoj4816] [Sdoi2017]数字表格

Description

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项,那么

f[0]=0

f[1]=1

f[n]=f[n-1]+f[n-2],n>=2

Doris用老师的超级计算机生成了一个n×m的表格,第i行第j列的格子中的数是f[gcd(i,j)],其中gcd(i,j)表示i,

j的最大公约数。Doris的表格中共有n×m个数,她想知道这些数的乘积是多少。答案对10^9+7取模。

Input

有多组测试数据。

第一个一个数T,表示数据组数。

接下来T行,每行两个数n,m

T<=1000,1<=n,m<=10^6

Output

输出T行,第i行的数是第i组数据的结果

Sample Input

3
2 3
4 5
6 7

Sample Output

1
6
960

solution

前置知识:莫比乌斯反演

题目让求的是:

[ans=prod_{i=1}^{n}prod_{j=1}^mf(gcd(i,j)) ]

然后枚举(gcd)的结果:

[ans=prod_{d=1}^{min(n,m)}pow(f(d),{sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}}[gcd(i,j)=1]) ]

其中,(pow(a,x))表示(a^x)其实只是因为美观

整理一下可得:

[egin{align} ans&=prod_{d=1}^{min(n,m)}pow(f(d),sum_{d_1}{mu(d_1)lfloorfrac{n}{dd_1} floorlfloorfrac{m}{dd_1} floor})\ &=prod_{T=1}^{min(n,m)}prod_{d|T}pow(f(d),mu(frac{T}{d})lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor)\ &=prod_{T=1}^{min(n,m)}pow(prod_{d|T}pow(f(d),mu(frac{T}{d})),lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor)\ end{align} ]

定义(g)为:

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

[egin{align} ans&=prod_{T=1}^{min(n,m)}pow(g(T),lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor) end{align} ]

然后(g(n))(貌似?)并不满足积性函数的性质,由于(n)较小,我们可以大力预处理(g),然后其他的函数筛一筛,在数论分块下就做完了。

时间复杂度(O(n*log(n)+q*sqrt{n}*log(n)))

后面那个log是快速幂的复杂度。

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

#define ll long long 

void read(int &x) {
	x=0;int f=1;char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
	for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
	if(x<0) x=-x,putchar('-');
	if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

const int maxn = 1e6+1;
const int mod = 1e9+7;

int pri[maxn],f[maxn],mu[maxn],g[maxn],tot,vis[maxn],invf[maxn],invg[maxn];

int qpow(int a,int x) {
	int res=1;
	for(;x;x>>=1,a=1ll*a*a%mod) if(x&1) res=1ll*res*a%mod;
	return res;
}

void sieve() {
	f[0]=0,f[1]=1;invf[1]=invf[0]=1;mu[1]=1;
	for(int i=2;i<maxn;i++) {
		f[i]=(f[i-1]+f[i-2])%mod;
		invf[i]=qpow(f[i],mod-2);
		if(!vis[i]) pri[++tot]=i,mu[i]=-1;
		for(int t,j=1;j<=tot&&i*pri[j]<maxn;j++) {
			vis[t=i*pri[j]]=1;
			if(!(i%pri[j])) {mu[t]=0;break;}
			mu[t]=-mu[i];
		}
	}
	for(int i=0;i<maxn;i++) g[i]=1;
	for(int i=1;i<maxn;i++) {
		if(!mu[i])  continue;
		for(int j=i;j<maxn;j+=i) 
			g[j]=1ll*g[j]*(mu[i]==1?f[j/i]:invf[j/i])%mod;
	}
	for(int i=1;i<maxn;i++) g[i]=1ll*g[i]*g[i-1]%mod;
	//for(int i=1;i<=10;i++) write(g[i]);
}

int main() {
	int t;read(t);
	//int PRE=clock();
	sieve();
	//cerr << (double)(clock()-PRE)/1e6 << endl;
	while(t--) {
		int n,m;read(n),read(m);
		int T=1,ans=1;if(n>m) swap(n,m);
		while(T<=n) {
			int pre=T;T=min(n/(n/T),m/(m/T));
			ans=1ll*ans*(qpow(qpow(1ll*g[T]*qpow(g[pre-1],mod-2)%mod,n/T),m/T))%mod;
			T++;
		}
		write(ans);
	}
	return 0;
}

原文地址:https://www.cnblogs.com/hbyer/p/10065498.html