[BZOJ4407]于神之怒加强版
Description
给下N,M,K.求
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m \gcd(i,j)^k \%(10^9+7)
\]
Input
输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。
Output
如题
Sample Input
1 2
3 3
Sample Output
20
HINT
1<=N,M,K<=5000000,1<=T<=2000
我们令n<m，然后继续推柿子
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m \gcd(i,j)^k
\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d=1}^n d^k [\gcd(i,j)=1]
\]
\[\sum\limits_{d=1}^n d^k\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum\limits_{x|i,x|j}\mu(x)
\]
\[\sum\limits_{d=1}^n d^k\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\lfloor\dfrac{n}{dx}\rfloor\lfloor\dfrac{n}{dx}\rfloor
\]
然后我们令\(T=dx\)，柿子变为
\[\sum\limits_{T=1}^n\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor\sum\limits_{d|T}d^k\mu(\dfrac{T}{d})
\]
令\(g(T)=\sum\limits_{d|T}d^k\mu(\dfrac{T}{d})\)，因此我们只要预处理出\(g(T)\)及其前缀和，我们便又可以开心地分块了
由于\(d^k\)是积性函数，\(g(T)\)也是积性函数，所以有
\[g(T)=\prod\limits_{i=1}^t g(P_i^{x_i})
\]
\[=\prod\limits_{i=1}^t(P_i^{k\times (x_i-1)}\times \mu(P_i)+P_i^{k\times x_i}\times \mu(1))
\]
\[=\prod\limits_{i=1}^t P_i^{k\times (x_i-1)}\times(P_i^k-1)
\]
然后就可以线筛了\(\sim\sim\sim\)
/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
int x=0,f=1;char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
inline void print(int x){
if (x>=10) print(x/10);
putchar(x%10+'0');
}
const int N=5e6,p=1e9+7;
int prime[N+10],miu[N+10],g[N+10],f[N+10];
bool inprime[N+10];
int Data,k,tot;
int mlt(int a,int b){
int res=1;
for (;b;b>>=1,a=1ll*a*a%p) if (b&1) res=1ll*res*a%p;
return res;
}
void prepare(){
f[1]=1;
for (int i=2;i<=N;i++){
if (!inprime[i]){
prime[++tot]=i;
g[tot]=mlt(i,k);
f[i]=(g[tot]-1+p)%p;
}
for (int j=1;j<=tot&&i*prime[j]<=N;j++){
inprime[i*prime[j]]=1;
if (i%prime[j]==0){
f[i*prime[j]]=1ll*f[i]*g[j]%p;
break;
}
f[i*prime[j]]=1ll*f[i]*f[prime[j]]%p;
}
}
for (int i=1;i<=N;i++) f[i]=(f[i]+f[i-1])%p;
}
int main(){
Data=read(),k=read();
prepare();
while (Data--){
int n=read(),m=read(),T=min(n,m),pos,Ans=0;
for (int i=1;i<=T;i=pos+1){
pos=min(n/(n/i),m/(m/i));
Ans=(Ans+1ll*(f[pos]-f[i-1]+p)*(n/i)%p*(m/i)%p)%p;
}
printf("%d\n",Ans);
}
return 0;
}
