bzoj4407: 于神之怒加强版

题目链接

bzoj4407: 于神之怒加强版

题解

求这个东西

[sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k ]

然后是套路

[egin{aligned}sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k &=sum_{d=1}^{min(n,m)}d^ksum_{i=1}^nsum_{j=1}^mleft[(i,j)=d ight]\ &=sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{min(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor)}mu(i)lfloorfrac{n}{id} floorlfloorfrac{m}{id} floor\ &=sum_{T=1}^{min(n,m)}lfloorfrac{n}{T} floorlfloorfrac{m}{T} floorsum_{dmid T}d^kmu(frac{T}{d})end{aligned} ]

显然,(sum_{dmid T}d^kmu(frac{T}{d}))它是个积性函数
尝试线筛
(F(T) = sum_{dmid T}d^kmu(frac{T}{d}))
考虑一个只有一个质因数的情况(F(p_i^{k_i}) = (p_i^{ki})^k - (p_i^{k_i - 1})^k)
所以(F(p_i^{k_i} imes p_i)=F(p_i^{k_i}) imes 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
",(solve(n,m) % mod + mod) % mod);  
	} 
	return 0; 
} 
原文地址:https://www.cnblogs.com/sssy/p/9525715.html