P4449 于神之怒加强版

Posted 2021-03-22 jackpei

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

(sum_{d=1}^n d^ksum_{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} floor lfloor frac{m}{dk} floor)

设 (T=dk)

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

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

然后 (pin prime,f(p^c)=p^{c imes (k-1)} imes (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
",(ans+M)%M);
  }
}
} signed main() {Luitaryi::main(); return 0;}

---

2019.12.16