首页 > 代码库 > BZOJ 4407 于神之怒加强版
BZOJ 4407 于神之怒加强版
题目链接：于神之怒加强版
这个式子还是很妙的，只是我已经思维僵化了
\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k \\
=&\sum_{g=1}^ng^k\sum_{i=1}^{\lfloor \frac{n}{g} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{g} \rfloor}\sum_{d|i,d|j}\mu(d) \\
=&\sum_{g=1}^ng^k\sum_{d=1}^{\lfloor \frac{n}{g} \rfloor}\mu(d)\lfloor \frac{n}{dg} \rfloor\lfloor \frac{m}{dg} \rfloor \\
=&\sum_{g=1}^ng^k\sum_{g|x}^{n}\mu(\frac{x}{g}) \lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor \\
=&\sum_{x=1}^n \lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor \sum_{g|x} g^k \mu(\frac{x}{g})
\end{aligned}
然后我们可以令\(f(x)=\sum_{g|x} g^k \mu(\frac{x}{g})\)，显然\(f(x)\)是\(id^k\)与\(\mu\)的狄利克雷卷积，所以\(f(x)\)也是一个积性函数
由于\(n\)的范围有\(5\times 10^6\)，所以我们来考虑一下如何筛这个\(f\)函数，也就是考虑在\(i\)是\(p\)的倍数的情况下计算\(f(ip)\)(\(p\)为质数)
令\(i=p^kq(q\perp p)\)，然后分情况讨论一下：
当\(q\ne 0\)时，由于\(i\)的最小质因子为\(p\)，那么\(p^{k+1}<i\)，所以\(f(\frac{i}{p^k})\times f(p^{k+1})\)就是所求；
当\(q=0\)时，由于\(f(p^x)=p^{kx}-p^{kx-k}\)，所以\(f(i)\times p^k\)即为所求，\(p^k\)则可以通过预处理所有质数的\(k\)次方解决。
最后再分块计算就好了。
下面贴代码：
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
#define maxn 5000010
#define mod 1000000007
using namespace std;
typedef long long llg;
int T,k,n,m,pr[maxn],lp;
int g[maxn],c[maxn];
bool vis[maxn];
llg f[maxn];
void gi(llg &x){if(x>=mod) x%=mod;}
llg mi(llg a,int b){
llg s=1;
while(b){
if(b&1) s=s*a,gi(s);
a=a*a,gi(a); b>>=1;
}
return s;
}
int main(){
File("a");
scanf("%d %d",&T,&k); f[1]=c[1]=1;
for(int i=2;i<maxn;i++){
if(!vis[i]) pr[++lp]=i,g[i]=i,c[i]=mi(i,k),f[i]=c[i]-1;
for(int j=1;i*pr[j]<maxn;j++){
if(i*pr[j]==31823){
int aa;
aa++;
}
vis[i*pr[j]]=1;
if(i%pr[j]) g[i*pr[j]]=pr[j],f[i*pr[j]]=f[i]*f[pr[j]]%mod;
else{
f[i*pr[j]]=g[i]!=i?f[i/g[i]]*f[g[i]*pr[j]]:f[i]*c[pr[j]];
gi(f[i*pr[j]]); g[i*pr[j]]=g[i]*pr[j]; break;
}
}
}
for(int i=2;i<maxn;i++) f[i]+=f[i-1],gi(f[i]);
while(T--){
scanf("%d %d",&n,&m);
if(n>m) n^=m^=n^=m;
llg ans=0;
for(int i=1,nt;i<=n;i=nt+1){
nt=min(n/(n/i),m/(m/i));
ans+=(f[nt]-f[i-1]+mod)*(n/i)%mod*(m/i)%mod;
if(ans>=mod) ans%=mod;
}
printf("%lld\n",ans);
}
return 0;
}
BZOJ 4407 于神之怒加强版