https://www.luogu.org/problemnew/show/P4449
(F(n)=sumlimits_{i=1}^{n}sumlimits_{i=1}^{m} gcd(i,j)^k)
首先加方括号,枚举g,提g:((min)表示(min(n,m)))
(sumlimits_{g=1}^{min} g^k sumlimits_{i=1}^{n} sumlimits_{i=1}^{m} [gcd(i,j)==g])
(sumlimits_{g=1}^{min} g^k sumlimits_{i=1}^{lfloorfrac{n}{g} floor} sumlimits_{i=1}^{lfloorfrac{m}{g} floor} [gcd(i,j)==1])
后面莫比乌斯反演:(k被你用了真恶心)
(sumlimits_{g=1}^{min} g^k sumlimits_{x=1}^{min} mu(x){lfloorfrac{n}{gx}
floor} {lfloorfrac{m}{gx}
floor})
众所周知,这种情况要枚举(T=gx):
(sumlimits_{T=1}^{min}sumlimits_{g|T} g^k mu(frac{T}{g}) {lfloorfrac{n}{T}
floor} {lfloorfrac{m}{T}
floor})
提T:
(sumlimits_{T=1}^{min} {lfloorfrac{n}{T}
floor} {lfloorfrac{m}{T}
floor} sumlimits_{g|T} g^k mu(frac{T}{g}))
所以:
(F(n)=sumlimits_{i=1}^{n}sumlimits_{i=1}^{m} gcd(i,j)^k = sumlimits_{T=1}^{min} {lfloorfrac{n}{T}
floor} {lfloorfrac{m}{T}
floor} sumlimits_{g|T} g^k mu(frac{T}{g}))
(n,m)这么小那我给你搞个线性筛吧。记 (G(n)=sumlimits_{g|n} g^k mu(frac{n}{g})) ,这个可以 (O(nlogn)) 筛出来。但是我偏偏要线性筛。
每次回答一个分块了事了。
小心取模,靠。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read() {
int x=0;
char c=getchar();
while(c<'0'||c>'9')
c=getchar();
do {
x=(x<<3)+(x<<1)+c-'0';
c=getchar();
} while(c>='0'&&c<='9');
return x;
}
inline void write(ll x) {
if(x>9) {
write(x/10);
}
putchar(x%10+'0');
return;
}
const int MAXN=5e6;
const int MOD=1e9+7;
int pri[MAXN+1];
int &pritop=pri[0];
ll G[MAXN+1];
int pk[MAXN+1];
int k;
inline ll qpow(ll x,int n) {
ll res=1;
while(n) {
if(n&1) {
res*=x;
if(res>=MOD)
res%=MOD;
}
x*=x;
if(x>=MOD)
x%=MOD;
n>>=1;
}
return res;
}
void sieve(int n=MAXN) {
pk[1]=1;
G[1]=1;
for(int i=2; i<=n; i++) {
if(!pri[i]) {
pri[++pritop]=i;
pk[i]=i;
G[i]=qpow(i,k)-1ll;
if(G[i]<0)
G[i]+=MOD;
}
for(int j=1; j<=pritop; j++) {
int &p=pri[j];
int t=i*p;
if(t>n)
break;
pri[t]=1;
if(i%p) {
pk[t]=pk[p];
//积性函数
G[t]=G[i]*G[p];
if(G[t]>=MOD)
G[t]%=MOD;
} else {
pk[t]=pk[i]*p;
if(pk[t]==t) {
//t是质数的幂次
G[t]=qpow(p,k)*G[i];
if(G[t]>=MOD)
G[t]%=MOD;
} else {
//积性函数
G[t]=G[pk[t]]*G[t/pk[t]];
if(G[t]>=MOD)
G[t]%=MOD;
}
break;
}
}
}
for(int i=1;i<=n;i++){
G[i]+=G[i-1];
if(G[i]>=MOD)
G[i]-=MOD;
}
}
inline ll ans(int n,int m) {
ll res=0;
int N=min(n,m);
for(int l=1,r;l<=N;l=r+1){
r=min(n/(n/l),m/(m/l));
ll tmp1=1ll*(n/l)*(m/l);
if(tmp1>=MOD)
tmp1%=MOD;
ll tmp2=G[r]-G[l-1];
if(tmp2<0)
tmp2+=MOD;
tmp1*=tmp2;
if(tmp1>=MOD)
tmp1%=MOD;
res+=tmp1;
if(res>=MOD)
res%=MOD;
}
return res;
}
inline void solve() {
int t=read();
k=read();
sieve();
while(t--) {
int n=read(),m=read();
write(ans(n,m));
putchar('
');
}
}
int main() {
#ifdef Yinku
freopen("Yinku.in","r",stdin);
#endif // Yinku
solve();
return 0;
}