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
ot=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}lfloorfrac{n}{t}
floor^2mu(frac{t}{d}))就会变成这样
[egin{aligned}
Ans=&sum_{d=1}^{n}xi^k(d)sum_{t|d}lfloorfrac{n}{t}
floor^2mu(frac{t}{d}) \
=&sum_{t=1}^{n}lfloorfrac{n}{t}
floor^2sum_{t|d}xi^k(d)mu(frac{t}{d}) \
end{aligned}
]
假设我们对于前面数论分块,那么现在要求
[egin{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}^{lfloorfrac{n}{i}
floor}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}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{n}{d}
floor}[gcd(i,j)==1]=2sum_{i=1}^{lfloorfrac{n}{d}
floor}phi(i)-1
]
这玩意不就是枚举其中一个数+去重么。
这应该是一个经典等式,然而我因为太菜没有做过。
于是原式变为
[Ans=sum_{d=1}^{n}xi^k(d)(2sum_{i=1}^{lfloorfrac{n}{d}
floor}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
",ans);
}
int main()
{
n=read();k=read();sieve();Min25();
return 0;
}