【BZOJ】【P4407】【于神之怒加强版】【题解】【数论】

传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4407

这两天刚好在给学弟讲数论，bzoj上就来了一道裸题……

答案就是\sum_D F(D)*n/d*m/d 

F(D=\sum{d|D} d^K mu(D/d)

线性筛，分块求

F(p)=p^K-1

F(x*p)=F(x)*F(p) (x,p)=1

F(x*p)=F(x)*p (x,p)/=1

Code:

#include<bits/stdc++.h>
using namespace std;
const int maxn=5e6+10;
typedef long long LL;
const LL MOD=1e9+7;
int T,K;
bool notp[maxn];
LL p[maxn];
LL f[maxn];
LL sum[maxn];
LL pw(LL x,LL k,LL p){
LL ans=1;
for(;k;k>>=1){
if(k&1)ans=ans*x%p;
x=x*x%p;
}return ans;
}
void pre(){
f[1]=1;
for(int i=2;i<maxn;i++){
if(!notp[i]){
f[i]=pw(i,K,MOD)-1;
p[++p[0]]=i;
}
for(int j=1;j<=p[0];j++){
if(i*p[j]>=maxn)break;
notp[i*p[j]]=1;
if(i%p[j]==0){
f[i*p[j]]=f[i]*pw(p[j],K,MOD)%MOD;
break;
}else{
f[i*p[j]]=f[i]*f[p[j]]%MOD;
}
}
}
for(int i=1;i<maxn;i++)
sum[i]=(sum[i-1]+f[i])%MOD;
}
int solve(int n,int m){
LL ans=0;
for(int i=1,j=0;i<=n;i=j+1){
j=min(n/(n/i),m/(m/i));
ans+=(sum[j]-sum[i-1]+MOD)%MOD*(n/i)%MOD*(m/i)%MOD;
ans%=MOD;
}return ans%MOD;
}
int main(){
scanf("%d%d",&T,&K);
pre();
// for(int i=1;i<=20;i++)
// cout<<i<<" "<<f[i]<<endl;
while(T--){
int n,m;scanf("%d%d",&n,&m);
if(n>m)swap(n,m);
printf("%d\n",int(solve(n,m)%MOD));
}
return 0;
}