BZOJ4816 [Sdoi2017]数字表格 【莫比乌斯反演】

题目

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取模。

输入格式

有多组测试数据。

第一个一个数T,表示数据组数。
接下来T行,每行两个数n,m
T<=1000,1<=n,m<=10^6

输出格式

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

输入样例

3

2 3

4 5

6 7

输出样例

1

6

960

题解

一道满是套路的莫比乌斯反演题

我们要求:

\[\prod\limits_{i=1}^{N}\prod\limits_{j=1}^{M} f[gcd(i,j)] \]

根据莫比乌斯反演的套路,我们将gcd改为枚举d

\[\prod\limits_{d=1}^{N} f[d]^{\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M} 1 [gcd(i,j) == d]} \]

然后把指数中的\(d\)提掉

\[\prod\limits_{d=1}^{N} f[d]^{\sum\limits_{i=1}^{\lfloor \frac{N}{d} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{M}{d} \rfloor} 1 [gcd(i,j) == 1]} \]

然后指数部分就是经典的莫比乌斯反演

\[\prod\limits_{d=1}^{N} f[d]^{\sum\limits_{i=1}^{\lfloor \frac{N}{d} \rfloor} \mu(i) * \lfloor \frac{N}{i * d} \rfloor * \lfloor \frac{M}{i * d} \rfloor} \]

如果只有一组询问,这样接近\(O(n)\)可以过,但是多组询问,我们考虑继续优化
我们有一个枚举\(i * d\)的套路
我们记\(T = i * d\)

\[\prod\limits_{T=1}^{N} \prod\limits_{d|T} f[d]^{\mu(\frac{T}{d}) * \lfloor \frac{N}{T} \rfloor * \lfloor \frac{M}{T} \rfloor} \]

划分一下:

\[\prod\limits_{T=1}^{N} ( \prod\limits_{d|T} f[d]^{\mu(\frac{T}{d})} ) ^ {\lfloor \frac{N}{T} \rfloor * \lfloor \frac{M}{T} \rfloor} \]

容易发现,内层的东西就是关于\(T\)的因子的积,根据经验,这玩意可以通过枚举\(O(nlogn)\)预处理
外层是只与\(T\)有关的整除,可以\(O(\sqrt{N})\)分块

那我们就用\(O(nlogn + T\sqrt{N})\)的复杂度做完了

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
#define res register
using namespace std;
const int maxn = 1000005,maxm = 100005,INF = 1000000000,P = 1000000007,md = 1000000006;
int fv[maxn],f[maxn],g[maxn],gv[maxn];
int mu[maxn],p[maxn],fac[maxn],pi;
int isn[maxn];
int qpow(int a,LL b){
	b = (b % md + md) % md;
	int ans = 1;
	for (; b; b >>= 1,a = (LL)a * a % P)
		if (b & 1) ans = (LL)ans * a % P;
	return ans;
}
void init(){
	mu[1] = 1;
	for (res int i = 2; i < maxn; i++){
		if (!isn[i]) p[++pi] = i,mu[i] = -1;
		for (res int j = 1; j <= pi && i * p[j] < maxn; j++){
			isn[i * p[j]] = true;
			if (i % p[j] == 0){
				mu[i * p[j]] = 0;
				break;
			}
			mu[i * p[j]] = -mu[i];
		}
	}
	f[0] = 0; f[1] = 1;
	fv[0] = 1; fv[1] = 1;
	for (res int i = 2; i < maxn; i++){
		f[i] = (f[i - 1] + f[i - 2]) % P;
		fv[i] = qpow(f[i],P - 2);
	}
	for (res int i = 1; i < maxn; i++) g[i] = 1;
	for (res int i = 1; i < maxn; i++){
		for (int j = i; j < maxn; j += i){
			if (mu[j / i] == 1) g[j] = (LL)g[j] * f[i] % P;
			if (mu[j / i] == -1) g[j] = (LL)g[j] * fv[i] % P;
		}
	}
	g[0] = gv[0] = 1;
	gv[1] = qpow(g[1],P - 2);
	for (res int i = 2; i < maxn; i++){
		g[i] = (LL)g[i] * g[i - 1] % P;
		gv[i] = qpow(g[i],P - 2);
	}
}
int main(){
	init();
	int n,m,T;
	LL ans;
	cin >> T;
	while (T--){
		cin >> n >> m;
		ans = 1;
		if (n > m) swap(n,m);
		for (res int i = 1,nxt; i <= n; i = nxt + 1){
			nxt = min(n / (n / i),m / (m / i));
			ans = (LL)ans * qpow((LL)g[nxt] * gv[i - 1] % P,(LL)(n / i) * (m / i)) % P;
		}
		cout << ans << endl;
	}
	return 0;
}

posted @ 2018-04-10 16:57  Mychael  阅读(80)  评论(0编辑  收藏