P4449 于神之怒加强版 （题解）

题目链接P4449 于神之怒加强版

题目大意：求

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

\(数据范围n,m\leq5e6\)
\(二话不说，先开导式子(假定n<m)：\)
\begin{aligned}
ans&=\sum _{i=1}^{n}\sum _{j=1}^{m}gcd(i,j) ^k\\
\end{aligned}

\[先套路地枚举d=gcd(i,j) \]

\begin{aligned}
ans&=\sum _{i=1}^{n}\sum _{j=1}^{m}\sum _{d=1}^{n}d ^k(d=gcd(i,j))
\end{aligned}

\[把i,j同时除一个d得到 \]

\begin{aligned}
ans&=\sum _{d=1}^{n}\sum _{i=1}^{\lfloor{\frac{n}{d}}\rfloor}\sum _{j=1}^{\lfloor{\frac{m}{d}}\rfloor}d ^k(gcd(i,j)=1)\\
ans&=\sum _{d=1}^{n}\sum _{i=1}^{\lfloor{\frac{n}{d}}\rfloor}\sum _{j=1}^{\lfloor{\frac{m}{d}}\rfloor}d ^k\sum _{e|gcd(i,j)}u(e)\\
ans&=\sum _{d=1}^{n}\sum _{e=1}^{\lfloor{\frac{n}{d}}\rfloor}{\lfloor{\frac{n}{de}}\rfloor}{\lfloor{\frac{m}{de}}\rfloor}d ^ku(e)\\
\end{aligned}

\[式子化到这里你就已经成功了十分之一了，因为这个式子的复杂度是n\sqrt n的，距离5e6的标准还差很多 \]

\[所以我们此处再次套路的使T=de 并改为枚举T,则有 \]

\begin{aligned}
ans&=\sum _{T=1}^{n}\sum _{d|T}\lfloor{\frac{n}{T}}\rfloor \lfloor{\frac{m}{T}}\rfloor d ^ku(\frac{T}{d})
\end{aligned}

\[换一下枚举顺序，把d的枚举放到后面 \]

\begin{aligned}
ans&=\sum _{T=1}^{n}\lfloor{\frac{n}{T}}\rfloor \lfloor{\frac{m}{T}}\rfloor\sum _{d|T} d ^ku(\frac{T}{d})
\end{aligned}

\[这时候的式子就顺眼多了，前面那个显然可以用除法分块在\sqrt n的时间内算出，这时考虑后面那一坨的性质 \]

\[首先显然的是f(d)=d^k是一个积性函数，u(d)更毋庸置疑，两大积性函数相乘，结果也显然是积性函数 \]

\[假定g(T)=\sum _{d|T} d ^ku(\frac{T}{d})那现在问题来了，怎么去利用线性筛去筛这个积性函数 \]

\[先单独考虑积性函数的性质，易知g(a*b)=g(a)*g(b)(gcd(a,b)=1) \]

\[再单独去考虑这个函数的一些性质： \]

\[首先在x为质数时g(x)=\sum _{d|x} d ^ku(\frac{x}{d})中的d既然可以整除x，所以d的值只能取1或x \]

\[首先d=1时，g(x)=1^k*u(x)=u(x)=-1 \]

\[其次d=x时，g(x)=x^k*u(1)=x^k*1=x^k \]

\[求和即为g(x)=x^k-1(x\in prime) \]

\[先导出一个结论：g(p^{x+1})=?(p \in prime),开导！ \]

\begin{aligned}
设op=p^{x+1},则有\\
g(op)=\sum _{d|op} d ^ku(\frac{op}{d})
\end{aligned}

\[因为op必然含有平方因子，所以当且仅当d=p^{x+1}和d=p^x时原式有值，分别为 \]

\[g(op)=p^{xk+k}(d=p^{x+1}),g(op)=-p^{xk}(d=p^x) \]

\[求和即为g(p^{x+1})=p^{xk+k}-p^{xk}(p \in prime) \]

\[随后即可导入最后一步g(x*prime)=? (prime|x) \]

\[因为要用g(x)的值去转移g(x*prime)，所以我们设g(x*prime)=g(x)*A，再次开导！！ \]

\[g(x*prime)=g(x) * A \]

\[感性理解一下，因为prime|x，所以说x*prime和x的质因数种类相同，只是质因数prime的数量多了一 \]

\[我们设a1,a2,a3...an为x不同的质数，c1,c2,c3...cn为其对应的数量，则有他们之间两两互质 \]

\[故g(x)=g(a1^{c1})*g(a2^{c2})*...*g(an^{cn})，带入到等式两边即为 \]

\[g(a1^{c1})*g(a2^{c2})*...*g(prime^{totprime+1})*...*g(an^{cn})=g(a1^{c1})*g(a2^{c2})*...*g(prime^{totprime})*...*g(an^{cn})*A \]

\[约分化简得g(prime^{totprime+1})=g(prime^{totprime})*A \]

\[这时候我们刚推得式子就有用了，带入得： \]

\[g(prime^{totprime+1})=g(prime^{totprime})*A \]

\[(prime^{totprime*k+k}-prime^{totprime*k})=(prime^{totprime*k}-prime^{totprime*k-k})*A \]

\[即可显然解得A=prime^k \]

\[至此，所有线性筛中的转移已经全部完成，再搭配上前缀和优化，完美的在O(n+T \sqrt n)的复杂度内完成了此题 \]

最后就是代码了：

#include<bits/stdc++.h>
using namespace std;
#define int long long 
const int N=5e6+5,p=1e9+7;
int n,m,k,t,x,y,z;
int s[N],f[N];
int prime[N],tot;
bool not_prime[N];
int qpow(int a,int b){
	int ans=1;
	while(b){
		if(b&1)ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}return ans;
}
void iint(){
	f[1]=1;
	for(int i=2;i<=5000000;i++){
		if(!not_prime[i])prime[++tot]=i,f[i]=(qpow(i,k)-1+p)%p;
		for(int j=1;j<=tot&&i*prime[j]<=5000000;++j){
			not_prime[i*prime[j]]=1;
			if(i%prime[j]==0){
				f[prime[j]*i]=(f[i]*((f[prime[j]]+1)%p))%p;
				break;
			}
			f[i*prime[j]]=f[i]*f[prime[j]]%p;
		}
	}
	for(int i=1;i<=5000000;i++)s[i]=(s[i-1]+f[i])%p;
}
int get(int n,int m){
	int ans=0;
	if(n>m)swap(n,m);
	for(int l=1,r=1;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans=(ans+(s[r]-s[l-1]+p)%p*(n/l)%p*(m/l)%p)%p;
	}return ans;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	// freopen("in.in","r",stdin);
	// freopen("out.out","w",stdout);
	cin>>t>>k;
	iint();
	while(t--){
		cin>>n>>m;
		cout<<get(n,m)<<endl;
	}
	return 0;
}

润喽~~~~~