P4449 于神之怒加强版

求：\(\sum_{i=1}^n \sum_{j=1}^m gcd(i,j)^k\)

\(\sum_{d=1}^n d^k\sum_{i=1}^{n/d} \sum_{j=1}^{m/d} \sum_{k|gcd(i,j)} \mu(k)\)

\(\sum_{d=1}^n d^k \sum_{k=1}^{n/d} \mu(x) \lfloor \frac{n}{dk} \rfloor \lfloor \frac{m}{dk} \rfloor\)

设 \(T=dk\)

\(\sum_{T=1}^n \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{d|T}d^k \mu(\frac{T}{d})\)

em，然后就发现前面可以除法分块，后面这一坨 \(f(T)=\sum_{d|T}d^k \mu(\frac{T}{d})\) 是积性函数，可以线性筛前缀和（雾

然后 \(p\in prime,f(p^c)=p^{c\times (k-1)} \times (p^k+\mu(p))\)

#include<bits/stdc++.h>
#define R register int
using namespace std;
namespace Luitaryi {
inline int g() { R x=0,f=1;
  register char s; while(!isdigit(s=getchar())) f=s=='-'?-1:f;
  do x=x*10+(s^48); while(isdigit(s=getchar())); return x*f;
} const int N=5000000,M=1000000007;
int n,m,K,T,p[N>>1],cnt,mx,f[N+10],h[N+10],nn[2010],mm[2010]; bool vis[N+10];
inline int qpow(int a,int b) { R ret=1;
  for(;b;b>>=1,a=1ll*a*a%M) if(b&1) ret=1ll*ret*a%M; return ret;  
}
inline void PRE(int N) { f[1]=1;
  for(R i=2;i<=N;++i) {
    if(!vis[i]) p[++cnt]=i,h[i]=qpow(i,K),f[i]=(h[i]-1+M)%M;
    for(R j=1,k;j<=cnt&&i*p[j]<=N;++j) {
      k=i*p[j]; vis[k]=true;
      if(i%p[j]==0) {
        f[k]=1ll*f[i]*h[p[j]]%M;
        break;
      } f[k]=1ll*f[i]*f[p[j]]%M;
    }
  } for(R i=1;i<=N;++i) f[i]=(f[i]+f[i-1])%M;
}
inline void main() {
  T=g(),K=g();
  for(R i=1;i<=T;++i) mx=max(mx,min(nn[i]=g(),mm[i]=g()));
  PRE(mx); 
  for(R i=1;i<=T;++i) { R ans=0;
    n=nn[i],m=mm[i]; if(n>m) swap(n,m);
    for(R l=1,r;l<=n;l=r+1) {
      r=min(n/(n/l),m/(m/l));
      ans=(ans+1ll*(f[r]-f[l-1])*(n/l)%M*(m/l))%M;
    } printf("%d\n",(ans+M)%M);
  }
}
} signed main() {Luitaryi::main(); return 0;}

---

2019.12.16