description
51nod
求$$\sum_{i=1}{n}\sum_{j=1}{n}sgcd(i,j)^k$$其中\(sgcd(i,j)\)表示\(i,j\)的次大公约数,如果\(gcd(i,j)=1\)那么\(sgcd(i,j)=0\)。
solution
记答案为\(Ans\)。
首先考虑直接枚举\(sgcd(i,j)\)。
\[Ans=\sum_{d=1}^{n}\xi^k(d)\sum_{i=1}^{n}\sum_{j=1}^n[gcd(i,j)==d]
\]
\]
其中当\(n\not=1\)时\(\xi^k(n)=\frac{n}{p_{min}(n)}\)。
这个时候如果你像菜鸡fdf一样瞎反演出
\(\sum_{i=1}^{n}\sum_{j=1}^n[gcd(i,j)==d]=\sum_{t|d}\lfloor\frac{n}{t}\rfloor^2\mu(\frac{t}{d})\)就会变成这样
\[\begin{aligned}
Ans=&\sum_{d=1}^{n}\xi^k(d)\sum_{t|d}\lfloor\frac{n}{t}\rfloor^2\mu(\frac{t}{d}) \\
=&\sum_{t=1}^{n}\lfloor\frac{n}{t}\rfloor^2\sum_{t|d}\xi^k(d)\mu(\frac{t}{d}) \\
\end{aligned}
\]
Ans=&\sum_{d=1}^{n}\xi^k(d)\sum_{t|d}\lfloor\frac{n}{t}\rfloor^2\mu(\frac{t}{d}) \\
=&\sum_{t=1}^{n}\lfloor\frac{n}{t}\rfloor^2\sum_{t|d}\xi^k(d)\mu(\frac{t}{d}) \\
\end{aligned}
\]
假设我们对于前面数论分块,那么现在要求
\[\begin{aligned}
\sum_{i=1}^{n}\sum_{n|d}\xi^k(d)\mu(\frac{d}{n})=&\sum_{i=1}^{n}\xi^k(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)\\
\end{aligned}
\]
\sum_{i=1}^{n}\sum_{n|d}\xi^k(d)\mu(\frac{d}{n})=&\sum_{i=1}^{n}\xi^k(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)\\
\end{aligned}
\]
你会发现这玩意又要数论分块一次。
然后菜鸡fdf就华丽地得到了一个\(O(n)\)的做法。O(n)可能也可以做
我们注意\(\sum_{i=1}^{n}\sum_{j=1}^n[gcd(i,j)==d]\)中\(i,j\)的上界都是\(n\)。
转化一下变成求
\[\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(i,j)==1]=2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\phi(i)-1
\]
\]
这玩意不就是枚举其中一个数+去重么。
这应该是一个经典等式,然而我因为太菜没有做过。
于是原式变为
\[Ans=\sum_{d=1}^{n}\xi^k(d)(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\phi(i)-1)
\]
\]
这样只要一次数论分块即可。
Code
我不管我不管我不会杜教筛只会暴力Min_25
#include<bits/stdc++.h>
#define FL "a"
using namespace std;
typedef unsigned int ll;
typedef long double dd;
const int N=1e5+10;
const int mod=1e9+7;
inline ll read(){
ll data=0,w=1;char ch=getchar();
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch<='9'&&ch>='0')data=data*10+ch-48,ch=getchar();
return data*w;
}
inline void file(){
freopen(FL".in","r",stdin);
freopen(FL".out","w",stdout);
}
inline ll poww(ll a,ll b){
ll res=1;
for(;b;b>>=1,a*=a)
if(b&1)res*=a;
return res;
}
int cnt;ll n,k,pri[N],pw[N][52],sum[N],sumk[N],s[52][52];bool vis[N];
inline void sieve(){
register int i,j;
for(vis[1]=1,i=2;i<N;i++){
if(!vis[i])pri[++cnt]=i;
for(j=1;j<=cnt&&i*pri[j]<N;j++){
vis[i*pri[j]]=1;if(i%pri[j]==0)break;
}
}
for(i=1;i<=cnt;i++)
for(j=pw[i][0]=1;j<=50;j++)pw[i][j]=pw[i][j-1]*pri[i];
for(i=1;i<=cnt;i++){
sum[i]=sum[i-1]+pri[i];
sumk[i]=sumk[i-1]+poww(pri[i],k);
}
s[0][0]=1;
for(i=0;i<=50;i++)
for(j=1;j<=i;j++)
s[i][j]=s[i-1][j-1]+j*s[i-1][j];
}
inline ll getsum(ll sn,ll sk){
register ll res=0,tmp,i,j;
for(i=0;i<=sk;i++){
tmp=1;
for(j=0;j<=i;j++)
if((sn-j+1)%(i+1)==0)tmp*=(sn-j+1)/(i+1);
else tmp*=(sn-j+1);
res+=s[sk][i]*tmp;
}
return res;
}
ll w[N],g[N],xi[N],ps[N],phi[N];int sqr,tot,id1[N],id2[N];
#define ID(x) (x<=sqr?id1[x]:id2[n/(x)])
inline void Min25(){
register ll i,j,a,b,e,ans;
for(sqr=sqrt(n),tot=0,i=1;i<=n;i=j+1){
j=n/(n/i);w[++tot]=n/i;
w[tot]<=sqr?id1[w[tot]]=tot:id2[n/w[tot]]=tot;
}
for(i=1;i<=tot;i++)g[i]=getsum(w[i],k)-1;
for(i=1;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++)
for(j=1;1ll*pri[i]*pri[i]<=w[j];j++){
g[j]-=pw[i][k]*(g[ID(w[j]/pri[i])]-sumk[i-1]);
xi[j]+=g[ID(w[j]/pri[i])]-sumk[i-1];
}
for(i=1;i<=tot;i++)g[i]=w[i]-1;
for(i=1;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++)
for(j=1;1ll*pri[i]*pri[i]<=w[j];j++)
g[j]-=g[ID(w[j]/pri[i])]-i+1;
for(i=1;i<=tot;i++)ps[i]-=g[i],xi[i]+=g[i];
for(i=1;i<=tot;i++)g[i]=1ll*w[i]*(w[i]+1)/2-1;
for(i=1;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++)
for(j=1;1ll*pri[i]*pri[i]<=w[j];j++)
g[j]-=pri[i]*(g[ID(w[j]/pri[i])]-sum[i-1]);
for(i=1;i<=tot;i++)ps[i]+=g[i];
for(i=cnt;i;i--)
if(1ll*pri[i]*pri[i]<=n)
for(j=1;1ll*pri[i]*pri[i]<=w[j];j++)
for(e=1,a=pri[i];1ll*a*pri[i]<=w[j];e++,a*=pri[i])
phi[j]+=(pw[i][e]-pw[i][e-1])*(phi[ID(w[j]/a)]+ps[ID(w[j]/a)]-(sum[i]-i))+(pw[i][e+1]-pw[i][e]);
for(i=1;i<=tot;i++)phi[i]+=ps[i]+1;
ans=0;
for(i=1;i<=n;i=j+1)
j=n/(n/i),ans+=(xi[ID(j)]-xi[ID(i-1)])*(2*phi[ID(n/i)]-1);
printf("%u\n",ans);
}
int main()
{
n=read();k=read();sieve();Min25();
return 0;
}