简要题意：
给定一个固定的 (k)，(T) 组询问求：
[sum_{i=1}^n sum_{j=1}^m gcd(i,j)^k
]
其中 (n,m,k leq 5 imes 10^6)，(T leq 2 imes 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}
floor} sum_{j=1}^{lfloor frac{n}{d}
floor} [gcd(i,j)==1]
]
[= sum_{i=1}^{min(n,m)} d^k sum_{i=1}^{lfloor frac{n}{d}
floor} sum_{j=1}^{lfloor frac{n}{d}
floor} sum_{k|gcd(i,j)} mu_k
]
[= sum_{i=1}^{min(n,m)} d^k sum_{k=1}^{lfloor frac{min(n,m)}{d}
floor} lfloor frac{n}{kd}
floor lfloor frac{m}{kd}
floor
]
[= sum_{i=1}^T lfloor frac{n}{T}
floor lfloor frac{m}{T}
floor sum_{d|T} d^k imes mu_{frac{T}{d}}
]
显然 (O(T sqrt{n})) 的时间足以通过，前面的部分整除分块足够。但是对于后面这一块较难解决，则设：
[g_x= sum_{d|x} d^k imes 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 imes (a_i-1)} imes mu_{p_i} + p_i^{k imes a_i} imes mu_1)
]
[= prod_{i=1}^q p_i^{k imes (a_i-1)} imes (p_i^k-1)
]
首先，上面对于 (g) 的几步需要解释。那个关于 (x) 的式子实际上是 (x) 质因数分解以后的结果，(p_i in ext{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('
');
}
return 0;
}
显然，你发现对于素数我们需要暴力计算 (p^k). 那么，首先 (5 imes 10^6) 以内的素数有多少个呢？代码测试：
得到结果约 (3.4 imes 10^5)，而 (log k = log (5 imes 10^6))，大概是 (22)，不会出问题，显然 (3.4 imes 10^5 imes 22) 是 (10^7) 级的，这也是程序较慢的原因（不过不影响通过，照样 (2.55s)）