bzoj4816/sdoi 2017 数字表格(莫比乌斯函数

在这里插入图片描述

分析

f ( i ) f(i) f(i) F i b o n a c c i Fibonacci Fibonacci 数列的第 i i i
题目本质上要求
a n s = ∏ i = 1 n ∏ j = 1 m f ( g c d ( i , j ) ) ans = \prod\limits_{i=1}^{n}\prod\limits_{j=1}^{m}f(gcd(i,j)) ans=i=1nj=1mf(gcd(i,j))
枚举 d = g c d ( i , j ) d = gcd(i,j) d=gcd(i,j)
a n s = ∏ d = 1 n f ( d ) ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = d ] ans = \prod\limits_{d=1}^{n}f(d)^{\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(i,j)==d]} ans=d=1nf(d)i=1nj=1m[gcd(i,j)==d]
对于右上角那部分,老生常谈了。
化简后要求的即是
a n s = ∏ d = 1 n f ( d ) ∑ x = 1 ⌊ n d ⌋ μ ( x ) ⌊ n d x ⌋ ⌊ m d x ⌋ ans = \prod\limits_{d=1}^{n}f(d)^{\sum\limits_{x=1}^{{\tiny \left\lfloor\dfrac{n}{d}\right\rfloor}}\mu(x){{\tiny\left\lfloor\dfrac{n}{dx}\right\rfloor}}{{\tiny\left\lfloor\dfrac{m}{dx}\right\rfloor}}} ans=d=1nf(d)x=1dnμ(x)dxndxm
然后套路枚举 d d d x x x 的乘积
a n s = ∏ T = 1 n ∏ d ∣ T f ( d ) μ ( T d ) ⌊ m T ⌋ ⌊ m T ⌋ ans = \prod\limits_{T=1}^{n}\prod\limits_{d|T}f(d)^{\mu({\tiny\dfrac{T}{d}}){{\tiny\left\lfloor\dfrac{m}{T}\right\rfloor}}{{\tiny\left\lfloor\dfrac{m}{T}\right\rfloor}}} ans=T=1ndTf(d)μ(dT)TmTm
预处理 h ( T ) = ∏ d ∣ T f ( d ) μ ( T d ) h(T) = \prod\limits_{d|T}f(d)^{\mu({\tiny\dfrac{T}{d}})} h(T)=dTf(d)μ(dT) 及其前缀积即可。每次询问分块+快速幂即可。
特别一提的是,当 μ \mu μ − 1 -1 1 时,乘上的应是 f ( d ) f(d) f(d) 的逆元。

代码如下

#include <bits/stdc++.h>
#define mod 1000000007
#define N 1000007
#define LL long long 
using namespace std;

LL z = 1;
int f[N], mu[N], x[N], p[N], g[N], s[N], inv[N], cnt, maxn = N - 7;
LL t; 
int ksm(int a, int b, int p){
	int s = 1;
	while(b){
		if(b & 1) s = z * s * a % mod;
		a = z * a * a % mod;
		b >>= 1;
	}
	return s;
}
int main(){
	int i, j, n, m, T, l, r, ji, ans;
	inv[1] = 1;
	for(i = 2; i <= maxn; i++){
		inv[i] = z * (mod - mod / i) * inv[mod % i] % mod;
	}
	mu[1] = 1;
	for(i = 2; i <= maxn; i++){
		if(!x[i]) x[i] = i,p[++cnt] = i, mu[i] = -1;
		for(j = 1; j <= cnt; j++){
			t = z * i * p[j];
			if(t > maxn) break;
			x[t] = p[j];
			if(i % p[j] == 0) break;
			mu[t] = -mu[i];
		}
	}
	f[1] = 1; 
	for(i = 2; i <= maxn; i++) f[i] = (f[i - 1] + f[i - 2]) % mod;
	for(i = 0; i <= maxn; i++) g[i] = 1;
	for(i = 1; i <= maxn; i++){
		for(j = 1; z * i * j <= maxn; j++){
			if(mu[j] == -1) g[i * j] = z * g[i * j] * ksm(f[i], mod - 2, mod) % mod;
			if(mu[j] == 1) g[i * j] = z * g[i * j] *  f[i] % mod;
		}
	}
	s[0] = 1;
	for(i = 1; i <= maxn; i++) s[i] = z * s[i - 1] * g[i] % mod;
	scanf("%d", &T);
	while(T--){
		ans = 1;
		scanf("%d%d", &n, &m);
		if(n > m) swap(n, m);
		for(l = 1; l <= n; l = r + 1){
			r = min(n / (n / l), m / (m / l));
			ji = z * s[r] * ksm(s[l - 1], mod - 2, mod) % mod;
			ji =ksm(ji, n / l, mod);
			ji = ksm(ji, m / l, mod);
			ans = z * ans * ji % mod;
		}
		printf("%d\n", ans);
	}
	return 0;
}
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值