Description
给出 N,K ,请计算下面这个式子:
其中,sgcd(i, j)表示(i, j)的所有公约数中第二大的,特殊地,如果gcd(i, j) = 1, 那么sgcd(i, j) = 0。
考虑到答案太大,请输出答案对2^32取模的结果.
1≤N≤109,1≤K≤50 1 ≤ N ≤ 10 9 , 1 ≤ K ≤ 50
样例解释:
因为gcd(i, j)=1时sgcd(i,j)=0对答案没有贡献,所以我们只考虑gcd(i,j)>1的情况.
当i是2时,j是2时,sgcd(i,j)=1,它的K次方是1
当i是2时,j是4时,sgcd(i,j)=1,它的K次方是1
当i是3时,j是3时,sgcd(i,j)=1,它的K次方是1
当i是4时,j是2时,sgcd(i,j)=1,它的K次方是1
当i是4时,j是4时,sgcd(i,j)=2,它的K次方是8
当i是5时,j是5时,sgcd(i,j)=1,它的K次方是1
Solution
本题要用Min_25筛,
显然,sgcd就是gcd再乘上它最小的质因子,
先上反演标准式子,假设只求
∑ni=1∑nj=1gcd(i,j)
∑
i
=
1
n
∑
j
=
1
n
gcd
(
i
,
j
)
如果你像我这样暴力,直接设 S(x,j) S ( x , j )
转移十分暴力,总复杂度很大(不会算),压线过,
如果你聪明一点,设
S(x,j)=∑xi=2(i/(i的最小质因子))k
S
(
x
,
j
)
=
∑
i
=
2
x
(
i
/
(
i
的
最
小
质
因
子
)
)
k
转移就简单多了,计算答案也很简单,复杂度也有保证,
Code
#include <cstdio>
#include <algorithm>
#include <cmath>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define efo(i,q) for(int i=A[q];i;i=B[i][0])
#define min(q,w) ((q)>(w)?(w):(q))
#define max(q,w) ((q)<(w)?(w):(q))
#define abs(q) ((q)<0?-(q):(q))
using namespace std;
typedef long long LL;
typedef unsigned int Uin;
const int N=2000500;
const LL mo=4294967296LL;
int read(int &n)
{
char ch=' ';int q=0,w=1;
for(;(ch!='-')&&((ch<'0')||(ch>'9'));ch=getchar());
if(ch=='-')w=-1,ch=getchar();
for(;ch>='0' && ch<='9';ch=getchar())q=q*10+ch-48;n=q*w;return n;
}
LL m,n;
LL ans;
int K;
int M;
bool prz[N];
int pr[N];
Uin prs[N],prk[N];
Uin phi[N],phi1[N];
Uin Ss[55][55],ds[55];
Uin s1[N],s2[N];
bool s1z[N],s2z[N];
LL d1[N];
Uin ksm(Uin q,int w)
{
Uin ans=1;
for(;w;w>>=1,q=q*q)if(w&1)ans=ans*q;
return ans;
}
void Pre(int n)
{
phi[1]=s1[1]=1;
fo(i,2,n)
{
if(!prz[i])pr[++pr[0]]=i,phi[i]=i-1,prs[pr[0]]=prs[pr[0]-1]+(s1[i]=ksm(i,K));
fo(j,1,pr[0])
{
int t=i*pr[j];
if(t>n)break;
prz[t]=1;
phi[t]=phi[i]*pr[j];
s1[t]=s1[i]*s1[pr[j]];
if(i%pr[j]==0)break;
phi[t]=phi[i]*(pr[j]-1);
}
}
fo(i,2,n)phi[i]+=phi[i-1],s1[i]+=s1[i-1];
}
Uin Gsum(Uin n)
{
if(!n)return 0;
if((int)n<=M)return s1[n];
else if(s2z[(::n)/n])return s2[(::n)/n];
if(n&1)ds[0]=n;
else ds[0]=n;
fo(i,1,K)
{
ds[i]=1;
for(Uin j=n+1;(int)n-(int)j<i;--j)if(j%(i+1))ds[i]*=j;
else ds[i]*=(Uin)j/(i+1);
fo(j,0,i-1)if((i+j)&1)ds[i]+=ds[j]*Ss[i][j];else ds[i]-=ds[j]*Ss[i][j];
}
s2[(::n)/n]=ds[K],s2z[(::n)/n]=1;
return ds[K];
}
Uin Gphi(Uin n)
{
if(n<=M)return phi[n];
if(phi1[(::n)/n])return phi1[(::n)/n];
Uin ans;
if(n&1)ans=((Uin)(n+1)/2)*n;
else ans=n/2*(Uin)(n+1);
for(Uin i=2,nx;i<=n;i=nx+1)
{
nx=n/(n/i);
ans=(ans-Gphi(n/i)*(nx-i+1))%mo;
}
return phi1[(::n)/n]=ans;
}
Uin f[N],f1[N];
Uin d[N],id[N],id1[2000];
Uin Gt;
Uin Gdoit(Uin n)
{
Uin ans=0;Gt=0;
for(Uin i=2,nx;i<=n;i=nx+1)
{
nx=n/(n/i);
Gt+=Gphi(n/i)*(Uin)(nx-i+1);
ans+=Gphi(n/i)*(Gsum(nx)-Gsum(i-1));
}
ans=ans*2;
Gt=Gt*2-n+1;
ans-=Gsum(n);
return ans+1;
}
int gcd(int x,int y){return y?(gcd(y,x%y)):x;}
int gv[N],gv1[2000];
Uin g[N],g1[2000];
Uin gg[N],gg1[2000];
Uin Gjian(int n,int m)
{
int S=1;
Uin ans=0;Gt=0;
if(n<=M)S=max(S,gv[n]+1),ans=g[n],Gt=gg[n];
else S=max(S,gv1[(::n)/n]+1),ans=g1[(::n)/n],Gt=gg1[(::n)/n];
for(int i=S,nx;i<=m;i+=nx)
{
for(nx=1;i+(nx<<1)<m&&n/pr[i+(nx<<1)-1]==n/pr[i];nx<<=1);
ans+=Gphi(n/pr[i])*(prk[nx+i-1]-prk[i-1]);
Gt+=Gphi(n/pr[i])*(Uin)(nx);
}
if(n<=M)gv[n]=m,g[n]=ans,gg[n]=Gt;
else gv1[(::n)/n]=m,g1[(::n)/n]=ans,gg1[(::n)/n]=Gt;
Gt=Gt*2-(Uin)m;
return ans*2-prs[m];
}
Uin Gs(Uin n)
{
for(Uin i=1,nx;i<=n;i=nx+1)
{
nx=n/(n/i);
int t=n/i;
d[++d[0]]=t;
if(t<=M)id[t]=d[0];
else id1[n/t]=d[0];
f[d[0]]=Gdoit(t);
f1[d[0]]=Gt;
}
Uin ans=0;
fo(j,1,pr[0]-1)
{
if(pr[j]*pr[j]>n)break;
int q=n/pr[j];
Uin t=0;
if(n/pr[j]<=M)t=f[id[q]];
else t=f[id1[n/(q)]];
t-=Gjian(q,j-1);
ans+=t;
Uin pk=prk[j]=ksm(pr[j],K);
prk[j]+=prk[j-1];
for(int i=1;i<=d[0]&&d[i]>=(LL)pr[j]*pr[j];++i)
{
int t=d[i]/pr[j];
if(t<=M)t=id[t];else t=id1[n/t];
f[i]-=(f[t]-Gjian(d[i]/pr[j],j-1))*pk;
f1[i]-=(f1[t]-Gt);
}
}
return ans+f1[1];
}
int main()
{
scanf("%lld%d",&n,&K);
M=2e6;
Pre(M);
Ss[1][1]=d1[1]=1;
fo(i,2,K)
{
fod(j,i,1)d1[j]=(d1[j-1]+(i-1LL)*d1[j])%mo,Ss[i][j]=abs(d1[j]);
}
printf("%u\n",Gs(n));
return 0;
}