BZOJ4816 SDOI2017 数字表格 莫比乌斯反演

传送门

---

做莫比乌斯反演题显著提高了我的(LaTeX)水平

推式子(默认(N leq M)，分数下取整，会省略大部分过程)

(egin{align*} prodlimits_{i=1}^N prodlimits_{j=1}^M f[gcd(i,j)] & = prodlimits_{d=1}^N f[d]^{sumlimits_{i=1}^frac{N}{d} sumlimits_{j=1}^frac{M}{d}[gcd(i,j)==1]} \ & = prodlimits_{d = 1}^N f[d]^{sumlimits_{p=1}^frac{N}{d} mu(p) frac{N}{dp} frac{M}{dp}} end{align*})

推到这里可以(O(TN))地通过两个数论分块得出答案，可以获得70pts

当然这还不够，按照老套路枚举(dp)继续推式子

(egin{align*} prodlimits_{d = 1}^N f[d]^{sumlimits_{p=1}^frac{N}{d} mu(p) frac{N}{dp} frac{M}{dp}} & = prodlimits_{T=1} ^ N prodlimits_{d | T} f[d] ^ {mu (frac{T}{d}) frac{N}{T} frac{M}{T}} \ & = prodlimits_{T=1} ^ N prodlimits_{d | T} (f[d] ^ {mu (frac{T}{d})} ) ^ { frac{N}{T} frac{M}{T}} end{align*})

(frac{N}{T} frac{M}{T})可以数论分块，所以要求(f[d]^{mu(frac{T}{d})})的前缀积

当你以为要线性筛的时候……(N leq 10^6)直接枚举倍数乘进去就行了……

总复杂度(O(NlogN + Tsqrt{N}))

#include<bits/stdc++.h>
//This code is written by Itst
using namespace std;

inline int read(){
    int a = 0;
    char c = getchar();
    bool f = 0;
    while(!isdigit(c)){
        if(c == '-')
            f = 1;
        c = getchar();
    }
    while(isdigit(c)){
        a = (a << 3) + (a << 1) + (c ^ '0');
        c = getchar();
    }
    return f ? -a : a;
}

const int MOD = 1e9 + 7 , MAXN = 1e6 + 7;
bool nprime[MAXN];
int fib[MAXN] , inv[MAXN] , mu[MAXN] , times[MAXN] , prime[MAXN];
int N , M , cnt;

inline int poww(long long a , int b){
	int times = 1;
	while(b){
		if(b & 1)
			times = times * a % MOD;
		a = a * a % MOD;
		b >>= 1;
	}
	return times;
}

void init(){
    fib[1] = inv[1] = times[1] = times[0] = mu[1] = 1;
    for(int i = 2 ; i <= N ; ++i){
    	times[i] = 1;
    	fib[i] = (fib[i - 1] + fib[i - 2]) % MOD;
    	inv[i] = poww(fib[i] , MOD - 2);
	}
	for(int i = 2 ; i <= N ; ++i){
		if(!nprime[i]){
			prime[++cnt] = i;
			mu[i] = -1;
		}
		for(int j = 1 ; j <= cnt && i * prime[j] <= N ; ++j){
			nprime[i * prime[j]] = 1;
			if(i % prime[j] == 0)
				break;
			mu[i * prime[j]] = -1 * mu[i];
		}
	}
	for(int i = 1 ; i <= N ; ++i)
		for(int j = 1 ; j * i <= N ; ++j)
			if(mu[j])
				times[j * i] = 1ll * times[j * i] * (mu[j] > 0 ? fib[i] : inv[i]) % MOD;
	for(int i = 2 ; i <= N ; ++i)
		times[i] = 1ll * times[i - 1] * times[i] % MOD;
	for(int i = 0 ; i <= N ; ++i)
		inv[i] = poww(times[i] , MOD - 2);
}

int main(){
	N = 1e6;
    init();
    for(int T = read() ; T ; --T){
		N = read();
    	M = read();
    	if(N > M)
    	    swap(N , M);
    	int ans = 1;
    	for(int i = 1 , pi ; i <= N ; i = pi + 1){
    		pi = min(N / (N / i) , M / (M / i));
    		ans = 1ll * ans * poww(1ll * times[pi] * inv[i - 1] % MOD , 1ll * (N / i) * (M / i) % (MOD - 1)) % MOD;
		}
		printf("%d
" , ans);
	}
    return 0;
}