P4449 于神之怒加强版 题解

CSDN同步

原题链接

简要题意：

给定一个固定的 \(k\)，\(T\) 组询问求：

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

其中 \(n,m,k \leq 5 \times 10^6\)，\(T \leq 2 \times 10^3\).

没什么好说的，没有部分分。根据 莫比乌斯反演性质 开始推式子！

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

\[= \sum_{i=1}^{\min(n,m)} d^k \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} [\gcd(i,j)==1] \]

\[= \sum_{i=1}^{\min(n,m)} d^k \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{k|\gcd(i,j)} \mu_k \]

\[= \sum_{i=1}^{\min(n,m)} d^k \sum_{k=1}^{\lfloor \frac{\min(n,m)}{d} \rfloor} \lfloor \frac{n}{kd} \rfloor \lfloor \frac{m}{kd} \rfloor \]

\[= \sum_{i=1}^T \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{d|T} d^k \times \mu_{\frac{T}{d}} \]

显然 \(O(T \sqrt{n})\) 的时间足以通过，前面的部分整除分块足够。但是对于后面这一块较难解决，则设：

\[g_x= \sum_{d|x} d^k \times \mu_{\frac{x}{d}} \]

通过观察可知这是个 积性函数。那么显然，对于：

\[x = \prod_{i=1}^q p_i^{a_i} \]

存在：

\[g_x = \prod_{i=1}^q g(p_i^{a_i}) \]

\[= \prod_{i=1}^q (p_i^{k \times (a_i-1)} \times \mu_{p_i} + p_i^{k \times a_i} \times \mu_1) \]

\[= \prod_{i=1}^q p_i^{k \times (a_i-1)} \times (p_i^k-1) \]

首先，上面对于 \(g\) 的几步需要解释。那个关于 \(x\) 的式子实际上是 \(x\) 质因数分解以后的结果，\(p_i \in \text{prime}\).

然后第一步是线性筛积性函数的套路。

第二步是一个变形，考虑展开 \(g(p_i^{a_i})\) 的过程.

第三步就是个乘法分配律，算出 \(\mu\) 而已。

好了，看看这个式子，每个部分都可以筛。结束。

时间复杂度：\(O(n + T \sqrt{n} + p \log k)\). \(p\) 的值之后解释。

实际得分：\(100pts\).

#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N=5e6+1;
const ll MOD=1e9+7;

inline int read(){char ch=getchar(); int f=1; while(ch<'0' || ch>'9') {if(ch=='-') f=-f; ch=getchar();}
	int x=0; while(ch>='0' && ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar(); return x*f;}

inline void write(ll x) {
	if(x<0) {putchar('-');write(-x);return;}
	if(x<10) {putchar(char(x%10+'0'));return;}
	write(x/10);putchar(char(x%10+'0'));
}

int prime[N],mu[N];
ll g[N],f[N]; bool h[N];
int T,k,cnt=0; ll ans=0;

inline ll pw(ll x,ll y) {
	ll ans=1; while(y) {
		if(y&1) ans=ans*x%MOD;
		x=x*x%MOD; y>>=1;
	} return ans;
}

inline void Euler(int n) {
	f[1]=1; for(int i=2;i<=n;i++) {
		if(!h[i]) prime[++cnt]=i,g[cnt]=pw(i,k),f[i]=(g[cnt]-1+MOD)%MOD;
			//f 表示 分解质因数后后面的一堆东西
		for(int j=1;j<=cnt && i*prime[j]<=n;j++) {
			h[i*prime[j]]=1;
			if(i%prime[j]==0) {f[i*prime[j]]=1ll*f[i]*g[j]%MOD;break;}
			f[i*prime[j]]=1ll*f[i]*f[prime[j]]%MOD; //线性筛模板
		}
	}
	for(int i=1;i<=n;i++) f[i]=(f[i]+f[i-1])%MOD; //前缀和
}

int main() {
	T=read(); k=read(); Euler(N-1);
	while(T--) {
		int n=read(),m=read();
		ll ans=0; int r;
		for(int i=1;i<=min(n,m);i=r+1) { //整除分块
			r=min(n/(n/i),m/(m/i));
			ans=(ans+1ll*(f[r]-f[i-1]+MOD)*(n/i)%MOD*(m/i)%MOD)%MOD;
		} write(ans); putchar('\n');
	}
	return 0;
}

显然，你发现对于素数我们需要暴力计算 \(p^k\). 那么，首先 \(5 \times 10^6\) 以内的素数有多少个呢？代码测试：

Link 代码

得到结果约 \(3.4 \times 10^5\)，而 \(\log k = \log (5 \times 10^6)\)，大概是 \(22\)，不会出问题，显然 \(3.4 \times 10^5 \times 22\) 是 \(10^7\) 级的，这也是程序较慢的原因（不过不影响通过，照样 \(2.55s\)）