bzoj4407: 于神之怒加强版

题目链接

bzoj4407: 于神之怒加强版

题解

求这个东西

\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k \]

然后是套路

\[\begin{aligned}\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k &=\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^n\sum_{j=1}^m\left[(i,j)=d\right]\\ &=\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor\\ &=\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d\mid T}d^k\mu(\frac{T}{d})\end{aligned} \]

显然，\(\sum_{d\mid T}d^k\mu(\frac{T}{d})\)它是个积性函数
尝试线筛
另\(F(T) = \sum_{d\mid T}d^k\mu(\frac{T}{d})\)
考虑一个只有一个质因数的情况\(F(p_i^{k_i}) = (p_i^{ki})^k - (p_i^{k_i - 1})^k\)
所以\(F(p_i^{k_i}\times p_i)=F(p_i^{k_i})\times p_i^{k}\)
对于不同质因数,因为是积性函数,直接乘起来

代码

#include<cstdio> 
#include<algorithm> 
inline int read() { 
	 int x = 0,f = 1 ; 
	 char c = getchar(); 
	 while(c < '0' || c > '9')c = getchar(); 
	 while(c <= '9' && c >= '0') x = x * 10 + c - '0',c = getchar(); 
	 return x * f; 
} 
const int mod = 1e9 + 7; 
int t,k; 
const int maxn = 5000007; 
int num,p[maxn],mu[maxn],F[maxn],tmp[maxn]; 
bool np[maxn]; 
int fstpow(int x,int k) {
	int ret = 1; 
	for(;k;k >>= 1,x = 1ll * x * x % mod) 
		if(k & 1) ret = 1ll * ret * x % mod; 
	return ret; 
} 
const int M = 5000000; 
void pre(int k) { 
	F[1] = mu[1] = 1; 
	for(int i = 2;i <= M;++ i) { 
		if(!np[i]) p[++ num] = i, mu[i] = -1, tmp[num] = fstpow(i,k),F[i] = (tmp[num] - 1) % mod; 
		for(int t,j = 1;j <= num && i * p[j] <= M;++ j) { 
			t = i * p[j]; 
			np[t] = 1; 
			if(i % p[j]) mu[t] = -mu[i],F[t] = 1ll * F[i] * F[p[j]] % mod; 
			else {mu[t] = 0,F[t] = 1ll * F[i] * tmp[j] % mod; } 
		} 
	} 
	for(int i = 2;i <= M;++ i) F[i] += F[i - 1],F[i] >= mod && (F[i] -= mod); 
} 
long long solve(int n,int m) { 
	long long ret = 0; 	
	for(int i = 1,nxt;i <= std::min(n,m);i = nxt + 1) { 
		nxt = std::min(n / (n / i),m / (m / i)); 
		ret += 1ll * (F[nxt] - F[i - 1] + mod) % mod * (n / i) % mod * (m / i) % mod; 
	} 
	return ret; 
} 
int main() { 
	t = read(),k = read(); 
	pre(k); 
	for(int n,m,i = 1;i <= t;++ i) { 
		n = read(),m = read(); 
		printf("%lld\n",(solve(n,m) % mod + mod) % mod);  
	} 
	return 0; 
}