奇怪的数学题
Description
Solution
考虑枚举
gcd
g
c
d
。
minp(d)
m
i
n
p
(
d
)
表示
d
d
的最小质因子。
一个很直观的思想就是分块,对于这一部分就是经典的杜教筛问题。
关于前面的部分,可以用
Min_25
M
i
n
_
25
筛处理,但注意到
f(d)=(dminp(d))
f
(
d
)
=
(
d
m
i
n
p
(
d
)
)
并非是个积性函数,其实只需维护两个函数即可,即维护
S(n,i)=∑ni=2f(d)k[i的最小质因子不小于pi或i是质数]
S
(
n
,
i
)
=
∑
i
=
2
n
f
(
d
)
k
[
i
的
最
小
质
因
子
不
小
于
p
i
或
i
是
质
数
]
S′(n,i)=∑ni=2dk[i的最小质因子不小于pi或i是质数]
S
′
(
n
,
i
)
=
∑
i
=
2
n
d
k
[
i
的
最
小
质
因
子
不
小
于
p
i
或
i
是
质
数
]
g,g′
g
,
g
′
也是同理。
于是我们便能得到转移式
还有一个细节,由于没有逆元求自然数幂和要用第二类斯特林数,取模时最好用自然溢出。
Code
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)
using namespace std;
typedef unsigned long long ll;
const ll N=12e5,M=1e6,mo=((ll)1<<32)-1,mm=mo+1,U=100;
ll n;
int p[N],k,oo,phi[N];
ll g[N],h[N],zh[N],w[N],f1[N],f2[N];
bool bz[N];
int id1[M],id2[M];
ll xs[U];
ll f[U][U];
inline ll ksm(ll o,ll t)
{
ll y=1;
for(;t;t>>=1,o=o*o)
if(t&1)y=y*o;
return y;
}
inline void prepare()
{
f[1][1]=1; f[1][0]=0;
fo(i,2,k){
f[i][0]=0;
fo(l,1,i)f[i][l]=(f[i-1][l-1]+f[i-1][l]*l);
}
fo(i,0,k)xs[i]=f[k][i];
}
inline ll qq(ll o)
{
ll ans=0; ++o;
fo(i,1,k){
int ok=0; ll lj=1;
fo(l,0,i)if((o-l)%(i+1)==0&&ok==0)lj=(lj*((o-l)/(i+1))),ok=1;
else lj=lj*(o-l);
ans=(ans+lj*xs[i]);
}
return ans;
}
ll ask(ll o)
{
if(o<=M)return f1[o];
ll k=n/o;
if(f2[k])return f2[k];
f2[k]=(o*(o+1)>>1);
ll l=2,r;
while(l<=o){
r=o/(o/l);
f2[k]=f2[k]-(r-l+1)*ask(o/l);
l=r+1;
}
return f2[k];
}
int main()
{
cin>>n>>k;
prepare();
fo(i,2,M){
if(!bz[i])p[++oo]=i,phi[i]=i-1,zh[oo]=(zh[oo-1]+ksm(i,k));
fo(j,1,oo)if(i*p[j]<=M){
bz[i*p[j]]=true;
phi[i*p[j]]=phi[i]*(p[j]-1);
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
}else break;
}
phi[1]=1;
fo(i,1,M)f1[i]=(f1[i-1]+phi[i]);
ll P=sqrt(n); int u=0,m=0;
fo(i,1,oo)if(p[i]<=P)u=i;else break;
ll l=1,r;
while(l<=n){
ll len=n/l; r=n/len;
if(len<=P)id1[len]=++m;else id2[r]=++m;
g[m]=len-1; h[m]=qq(len)-1; w[m]=len;
l=r+1;
}
fo(i,1,u){
register ll ww=ksm(p[i],k),U=(ll)p[i]*p[i]; register int op;
fo(j,1,m)if(U<=w[j])
{
op=(w[j]/p[i]<=P)?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
g[j]=g[j]-(g[op]-i+1);
h[j]=h[j]-(h[op]-zh[i-1])*ww;
}else break;
}
fd(i,u,1){
register ll ww=ksm(p[i],k),U=(ll)p[i]*p[i];
fo(j,1,m)if(U<=w[j]){
register ll uu=w[j]/p[i];
register int ok;
for(register ll o=p[i],o3=ww,u=w[j]/o;o<=uu;o=o*p[i],o3=o3*ww,u=u/p[i])
{
ok=(u<=P)?id1[u]:id2[n/u];
h[j]=h[j]+(h[ok]-zh[i])*o3;
h[j]=h[j]+o3*ww;
}
for(register ll o1=1,o2=p[i],u=w[j]/p[i];o2<=uu;o1=o1*ww,o2=o2*p[i],u=u/p[i])
{
ok=(u<=P)?id1[u]:id2[n/u];
g[j]=g[j]+(h[ok]-zh[i])*o1;
g[j]=g[j]+o1*ww;
}
}else break;
}
l=1; ll ans=0;
while(l<=n){
r=n/(n/l);
int o1,o2;
o2=(r<=P)?id1[r]:id2[n/r];
o1=l-1<=P?id1[l-1]:id2[n/(l-1)];
ll cz=g[o2]-g[o1];
ll ok=2*ask(n/r)-1;
ans=ans+cz*ok;
l=r+1;
}
ans=ans&mo;
cout<<ans;
}