题解 P4449 于神之怒加强版
这道题算是我完完整整推的第一道题，写篇题解纪念一下。
废话不多说，直接开始推式子（给新手准备，过程较详细，大佬可自行跳步），以下过程中均假设 \((n\le m)\)，\([d=1]\) 类似于代码中的 (d==1) 。
\[\sum_{i=1}^{n}\sum_{j=1}^mgcd(i,j)^k
\]
直接按套路提取
\[\sum_{d=1}^{n}d^k\sum_{i=1}^{n}\sum_{j=1}^m[gcd(i,j)=d]
\]
等价于
\[\sum_{d=1}^{n}d^k\sum_{i=1}^{n}\sum_{j=1}^m[gcd(\frac{i}{d},\frac{j}{d})=1]
\]
我们改变一下枚举项，枚举 \(\frac{i}{d},\frac{j}{d}\)
\[\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]
\]
接下来就要用到莫比乌斯函数的性质
\[\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{x|i}\sum_{x|j}\mu(x)
\]
我们改成枚举 \(\frac{i}{x},\frac{j}{x}\)
\[\sum_{d=1}^{n}d^k\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor\lfloor\frac{m}{dx}\rfloor\mu(x)
\]
这个 \(dx\) 让人很不爽，我们把它换个元，令 \(T=dx\)
\[\sum_{T=1}^n\sum_{d|T}d^k\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\mu(\frac{T}{d})
\]
把无关项往前提一提
\[\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}d^k\mu(\frac{T}{d})
\]
整理好后我们会发现这个式子太友善了，后面一部分
\[\sum_{d|T}d^k\mu(\frac{T}{d})
\]
完全就是 \(id_k*\mu\) ( \(*\) 为 \(Dirichlet\) 卷积)，所以我们可以知道 \(f(T)=id_k*\mu\) 是一个积性函数。
所以我们只需要研究 \(T|p^x\kern 0.4emx\in N_+ \kern 0.4emp\in prime\)
而由于 \(\mu\) 的性质，(当 \(x=1\) 时)
\[d^k\mu(\frac{T}{d})=\left\{
\begin{aligned}
-1\kern 1.0em(d=1)\\
T^k\kern 1.0em(d=T)\\
\end{aligned}
\right.
\]
所以 \(f(T)=T^k-1\)
当 \(x>1\) 时，有
\[d^k\mu(\frac{T^x}{d})=\left\{
\begin{array}{lcl}
-(T^{x-1})^k\kern 1.0em(d=T^{x-1})\\
(T^x)^k\kern 2.7em(d=T^x)\\
\end{array}
\right.
\]
所以我们线性筛时
\[f_{i×prime_j}=\left\{
\begin{array}{lcl}
f_i×f_{prime_j}\kern 1.0em gcd(prime_j,i)=1\\
f_i×prime_j^k\kern 1.0em (prime_j|i)
\end{array}
\right.
\]
这样我们就可以愉快的线性筛了，而由于那个快被我们遗忘的前半部分用数论分块，所以总体复杂度为 \(O(n+(\text{素数个数})×logk+T×\sqrt{n})\)
\(AC\kern 0.4em CODE:\)
#include<bits/stdc++.h>
#define ri register int
#define p(i) ++i
using namespace std;
const int MOD=1e9+7,N=5e6+7;
inline int read() {
int x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
return x*f;
}
int prime[N],vis[N],f[N],cm[N],k;
inline int fpow(int x,int y) {
int res=1;
while(y) {
if (y&1) res=1ll*res*x%MOD;
x=1ll*x*x%MOD;y>>=1;
}
return res;
}
void getmju(int n) {
int tot=0;f[1]=1;
for (ri i(2);i<=n;p(i)) {
if (!vis[i]) prime[p(tot)]=i,cm[tot]=fpow(i,k),f[i]=1ll*(cm[tot]-1+MOD)%MOD;
for (ri j(1);j<=tot&&prime[j]*i<=n;p(j)) {
vis[prime[j]*i]=1;
if (!(i%prime[j])) {
f[i*prime[j]]=1ll*cm[j]*f[i]%MOD;
break;
}
f[i*prime[j]]=1ll*f[i]*f[prime[j]]%MOD;
}
}
for (ri i(2);i<=n;p(i)) f[i]=(f[i]+f[i-1])%MOD;
}
int n[2020],m[2020],mxn;
int main() {
//freopen("rr.out","w",stdout);
int T=read();k=read();
for (ri i(1);i<=T;p(i)) n[i]=read(),m[i]=read(),mxn=max(mxn,min(n[i],m[i]));//小小的优化一下
getmju(mxn);
for (ri i(1);i<=T;p(i)) {
int nn=n[i],mm=m[i];
if (nn>mm) swap(nn,mm);
int ans=0;
for (ri l(1),r;l<=nn;l=r+1) {
r=min(nn/(nn/l),mm/(mm/l));
ans=(ans+1ll*((f[r]-f[l-1]+MOD)%MOD)*(nn/l)%MOD*(mm/l)%MOD)%MOD;
}
printf("%d\n",ans);
}
return 0;
}
