BZOJ4407 :于神之怒加强版

题面

戳我

Sol

(ans=sum_{d=1}^{N}d^k*sum_{i=1}^{lfloorfrac{N}{d} floor}mu(i)*lfloorfrac{N}{d*i} floor*lfloorfrac{M}{d*i} floor)
(将d*i换成S)
(原式=sum_{S=1}^{N}(lfloorfrac{N}{S} floor)*(lfloorfrac{M}{S} floor)*sum_{i|S}(frac{S}{i})^k*mu(i))
(设f(n)=sum_{i|n}(frac{n}{i})^k*mu(i)),它是个积性函数,可以线性筛

# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(5e6 + 1), MOD(1e9 + 7);

IL ll Read(){
    char c = '%'; ll x = 0, z = 1;
    for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;
    for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';
    return x * z;
}

int prime[_], num, mu[_], f[_], k, po[_], s[_];
bool isprime[_];

IL int Pow(RG ll x, RG ll y){
	RG ll ret = 1;
	for(; y; y >>= 1, x = x * x % MOD) if(y & 1) ret = ret * x % MOD;
	return ret;
}

IL void Prepare(){
	isprime[1] = 1; s[1] = f[1] = 1;
	for(RG int i = 2; i < _; ++i){
		if(!isprime[i]) prime[++num] = i, po[i] = Pow(i, k), f[i] = (po[i] - 1 + MOD) % MOD;
		for(RG int j = 1; j <= num && i * prime[j] < _; ++j){
			isprime[i * prime[j]] = 1;
			if(i % prime[j]) f[i * prime[j]] = 1LL * f[i] * f[prime[j]] % MOD;
			else{  f[i * prime[j]] = 1LL * f[i] * po[prime[j]] % MOD; break;  }
		}
		s[i] = (f[i] + s[i - 1]) % MOD;
	}
}

int main(RG int argc, RG char *argv[]){
	RG ll T = Read(), n, m; k = Read(); Prepare();
	for(; T; --T){
		RG ll ans = 0; n = Read(); m = Read();
		if(n > m) swap(n, m);
		for(RG ll i = 1, j; i <= n; i = j + 1){
			j = min(n / (n / i), m / (m / i));
			(ans += 1LL * (s[j] - s[i - 1] + MOD) % MOD * (n / i) % MOD * (m / i) % MOD) %= MOD;
		}
		printf("%lld
", (ans + MOD) % MOD);
	}
	return 0;
}

原文地址:https://www.cnblogs.com/cjoieryl/p/8268628.html