前言
万年不写公开博客了，这次填个坑
题目相关
题目大意
求
∑
i
=
1
N
∑
j
=
1
N
s
g
c
d
(
i
,
j
)
k
(
m
o
d
2
32
)
\sum_{i=1}^N\sum_{j=1}^Nsgcd(i,j)^k\pmod {2^{32}}
i=1∑Nj=1∑Nsgcd(i,j)k(mod232)
sgcd(i,j)为i,j的第二大公约数
特殊的，若gcd(i,j)=1,则sgcd(i,j)=0
数据范围
n ≤ 1 0 9 , k ≤ 50 n\le10^9,k\le50 n≤109,k≤50
题解
设m(x)为x的次大因数，特殊的，m(1)=0
∑
i
=
1
N
∑
j
=
1
N
s
g
c
d
(
i
,
j
)
k
=
∑
i
=
1
N
∑
j
=
1
N
m
(
g
c
d
(
i
,
j
)
)
k
根
据
定
义
=
∑
g
=
1
N
m
(
g
)
k
∑
i
=
1
⌊
N
g
⌋
∑
j
=
1
⌊
N
g
⌋
[
g
c
d
(
i
,
j
)
=
=
1
]
提
出
g
c
d
=
∑
g
=
1
N
m
(
g
)
k
(
2
(
∑
i
=
1
⌊
N
g
⌋
φ
(
i
)
)
−
1
)
大
小
互
换
也
可
以
所
以
要
乘
以
2
，
并
且
(
1
,
1
)
只
能
算
一
次
所
以
减
一
=
∑
t
=
1
n
∑
⌊
N
g
⌋
=
t
m
(
g
)
k
(
2
(
∑
i
=
1
t
φ
(
i
)
)
−
1
)
\begin{aligned} \sum_{i=1}^N\sum_{j=1}^Nsgcd(i,j)^k&=\sum_{i=1}^N\sum_{j=1}^Nm(gcd(i,j))^k&根据定义\\ &=\sum_{g=1}^Nm(g)^k\sum_{i=1}^{\left\lfloor\frac{N}{g}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{N}{g}\right\rfloor}[gcd(i,j)==1]&提出gcd\\ &=\sum_{g=1}^Nm(g)^k\left(2\left(\sum_{i=1}^{\left\lfloor\frac{N}{g}\right\rfloor}\varphi(i)\right)-1\right)&大小互换也可以所以要乘以2，并且(1,1)只能算一次所以减一\\ &=\sum_{t=1}^n\sum_{\left\lfloor\frac{N}{g}\right\rfloor=t}m(g)^k\left(2\left(\sum_{i=1}^{t}\varphi(i)\right)-1\right)\\ \end{aligned}
i=1∑Nj=1∑Nsgcd(i,j)k=i=1∑Nj=1∑Nm(gcd(i,j))k=g=1∑Nm(g)ki=1∑⌊gN⌋j=1∑⌊gN⌋[gcd(i,j)==1]=g=1∑Nm(g)k⎝⎜⎛2⎝⎜⎛i=1∑⌊gN⌋φ(i)⎠⎟⎞−1⎠⎟⎞=t=1∑n⌊gN⌋=t∑m(g)k(2(i=1∑tφ(i))−1)根据定义提出gcd大小互换也可以所以要乘以2，并且(1,1)只能算一次所以减一
我们发现，满足条件的t只有
O
(
n
)
\mathcal O(\sqrt n )
O(n)个，而最后一项可以通过杜教筛来做
容易发现，对于同一个t，其满足条件的g是连续的一串，那么只要对于所有满足条件的t求
∑
g
=
1
m
a
x
g
m
(
g
)
k
\sum_{g=1}^{maxg}m(g)^k
∑g=1maxgm(g)k，这个用min_25筛算
因为本质上min25筛会枚举最小因数，所以把最小因数排除掉即可
在预处理上述运算的时候要求
∑
g
=
1
n
g
k
\sum_{g=1}^{n}g^k
∑g=1ngk
我们引入斯特林数
x
n
=
∑
i
=
0
x
{
n
i
}
x
i
↓
x^n=\sum_{i=0}^x\begin{Bmatrix}n\\i\end{Bmatrix}x^{i\downarrow}
xn=i=0∑x{ni}xi↓
证明可以点进那篇博客看
这样的话，求sigma就是
∑
x
=
1
n
x
k
=
∑
x
=
1
n
∑
i
=
0
x
{
k
i
}
x
i
↓
=
∑
i
=
0
n
{
k
i
}
∑
x
=
i
n
x
i
↓
=
∑
i
=
0
n
{
k
i
}
∑
x
=
i
n
(
x
i
)
i
!
=
∑
i
=
0
n
i
!
{
k
i
}
∑
x
=
i
n
(
x
i
)
=
∑
i
=
0
n
i
!
{
k
i
}
(
n
+
1
i
+
1
)
=
∑
i
=
0
n
i
!
{
k
i
}
(
n
+
1
)
!
(
i
+
1
)
!
(
n
−
i
)
!
=
∑
i
=
0
n
{
k
i
}
(
n
+
1
)
!
(
i
+
1
)
(
n
−
i
)
!
=
∑
i
=
0
n
{
k
i
}
(
n
+
1
)
i
+
1
↓
i
+
1
\begin{aligned} \sum_{x=1}^nx^k&=\sum_{x=1}^n\sum_{i=0}^x\begin{Bmatrix}k\\i\end{Bmatrix}x^{i\downarrow}\\ &=\sum_{i=0}^n\begin{Bmatrix}k\\i\end{Bmatrix}\sum_{x=i}^nx^{i\downarrow}\\ &=\sum_{i=0}^n\begin{Bmatrix}k\\i\end{Bmatrix}\sum_{x=i}^n\binom{x}{i}i!\\ &=\sum_{i=0}^ni!\begin{Bmatrix}k\\i\end{Bmatrix}\sum_{x=i}^n\binom{x}{i}\\ &=\sum_{i=0}^ni!\begin{Bmatrix}k\\i\end{Bmatrix}\binom{n+1}{i+1}\\ &=\sum_{i=0}^ni!\begin{Bmatrix}k\\i\end{Bmatrix}\frac{(n+1)!}{(i+1)!(n-i)!}\\ &=\sum_{i=0}^n\begin{Bmatrix}k\\i\end{Bmatrix}\frac{(n+1)!}{(i+1)(n-i)!}\\ &=\sum_{i=0}^n\begin{Bmatrix}k\\i\end{Bmatrix}\frac{(n+1)^{i+1\downarrow}}{i+1}\\ \end{aligned}
x=1∑nxk=x=1∑ni=0∑x{ki}xi↓=i=0∑n{ki}x=i∑nxi↓=i=0∑n{ki}x=i∑n(ix)i!=i=0∑ni!{ki}x=i∑n(ix)=i=0∑ni!{ki}(i+1n+1)=i=0∑ni!{ki}(i+1)!(n−i)!(n+1)!=i=0∑n{ki}(i+1)(n−i)!(n+1)!=i=0∑n{ki}i+1(n+1)i+1↓
此处复杂度
O
(
k
2
)
\mathcal O(k^2)
O(k2)
总复杂度
O
(
n
3
4
l
o
g
n
+
n
k
2
+
n
2
3
)
\mathcal O(\frac{n^{\frac34}}{log\ n }+\sqrt nk^2+n^{\frac23})
O(log nn43+nk2+n32)
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<map>
#define rg register
typedef long long LL;
typedef unsigned int ULL;
template <typename T> inline T max(const T a,const T b){return a>b?a:b;}
template <typename T> inline T min(const T a,const T b){return a<b?a:b;}
template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;}
template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;}
template <typename T> inline T abs(const T a){return a>0?a:-a;}
template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);}
template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;}
template <typename T> inline T square(const T x){return x*x;};
template <typename T> inline void read(T&x)
{
char cu=getchar();x=0;bool fla=0;
while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
if(fla)x=-x;
}
template <typename T> inline void printe(const T x)
{
if(x>=10)printe(x/10);
putchar(x%10+'0');
}
template <typename T> inline void print(const T x)
{
if(x<0)putchar('-'),printe(-x);
else printe(x);
}
namespace du
{
const int MAX=2000000;
bool isprime[MAX+1];
int n,prime[MAX+1],primesize;LL phi[MAX+1],phi1[MAX+1];
inline void getlist(const LL listsize)
{
memset(isprime,1,sizeof(isprime));
isprime[1]=0,phi[1]=1;;
for(rg int i=2;i<=listsize;i++)
{
if(isprime[i])prime[++primesize]=i,phi[i]=i-1;
for(rg int j=1;j<=primesize&&(LL)i*prime[j]<=listsize;j++)
{
isprime[i*prime[j]]=0;
if(i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
}
inline LL Val(const int x)
{
if(x<=MAX)return phi[x];
return phi1[n/x];
}
void pre(int ttt)
{
n=ttt;
getlist(MAX);
for(rg int i=1;i<=MAX;i++)phi[i]+=phi[i-1];
for(rg int i=min((int)sqrt(n),n/(MAX+1));i>=1;i--)
{
const int u=n/i;
LL res=(LL)u*(u+1)/2;
int limit=sqrt(u);
for(rg int j=1;j<=limit;j++)res-=Val(u/j);
limit++;
for(rg int j=1;limit*j<=u;j++)res-=phi[j]*((u/j)-(u/(j+1)));
phi1[i]=res;
}
}
}
int n,k,part,tot,qz[700001],prime[700001],sq[700001],primesize,sum0[700001],id[700001];
inline int ck(const int x){return x<=part?x:tot-n/x+1;}
ULL S[51][51],sum2[700001],sum3[700001];
LL sum1[700001];
void pre()
{
S[0][0]=1;
for(rg int i=1;i<=k;i++)
for(rg int j=1;j<=i;j++)
S[i][j]=S[i-1][j-1]+j*S[i-1][j];
}
inline ULL qzktime(int x)
{
ULL ans=0,dq=0;
for(rg int i=0;i<=k;i++)
for(rg int j=x-i+1;j<=x+1;j++)
if(j%(i+1)==0)
{
dq=S[k][i];
for(rg int l=x-i+1;l<=x+1;l++)
if(j==l)dq*=l/(i+1);
else dq*=l;
ans+=dq;
break;
}
return ans;
}
int val[700001];
ULL R[700001],ans,QZ[700001];
ULL pow(ULL x,int y)
{
ULL res=1;
for(;y;y>>=1,x=x*x)if(y&1)res*=x;
return res;
}
int main()
{
read(n),read(k);
du::pre(n),pre();
part=sqrt(n);
tot=primesize=0;
for(rg int i=1;i<=part;i++)id[++tot]=i,sum0[tot]=i-1,sum1[tot]=((((LL)i+1)*i)>>1)-1,sum2[tot]=qzktime(i)-1,sum3[tot]=i-1;
id[++tot]=n/part;
if(id[tot]!=part)sum0[tot]=id[tot]-1,sum1[tot]=((LL)(((LL)id[tot]+1)*id[tot])>>1)-1,sum2[tot]=qzktime(id[tot])-1,sum3[tot]=id[tot]-1;
else tot--;
for(rg int i=part-1;i>=1;i--)id[++tot]=n/i,sum0[tot]=id[tot]-1,sum1[tot]=((((LL)id[tot]+1)*id[tot])>>1)-1,sum2[tot]=qzktime(id[tot])-1,sum3[tot]=id[tot]-1;
for(rg int i=2;i<=part;i++)
if(sum0[i]!=sum0[i-1])
{
const int limit=i*i;
ULL I=pow(i,k);
for(rg int j=tot;id[j]>=limit;j--)
{
const int t=ck(id[j]/i);
sum3[j]+=sum2[t]-QZ[primesize]-sum0[t]+primesize;
sum2[j]-=(sum2[t]-QZ[primesize])*I;
sum1[j]-=(sum1[t]-qz[primesize])*i;
sum0[j]-=sum0[t]-primesize;
}
prime[++primesize]=i,sq[primesize]=i*i,qz[primesize]=qz[primesize-1]+i,QZ[primesize]=QZ[primesize-1]+I;
}
for(rg int i=1;i<=tot;i++)sum1[i]-=sum0[i];//printf("%lld\n",sum1[i]);//sum1为phi的前缀和
for(rg int i=1;i<=part;i++)val[i]=n/i,R[i]=sum3[ck(val[i])];//printf("%d %u %lld\n",val[i],R[i],sum1[i]+1);
for(rg int i=part+1;i<=tot;i++)val[i]=tot-i+1,R[i]=sum3[ck(val[i])];//printf("%d %u %lld\n",val[i],R[i],sum1[i]+1);
for(rg int i=1;i<=part;i++)ans+=(R[i]-R[i+1])*(2*(ULL)du::Val(i)-1);//printf("%lld\n",sum1[i]);
for(rg int i=part+1;i<=tot;i++)ans+=(R[i]-R[i+1])*(2*(ULL)du::Val(n/(tot-i+1))-1);
print(ans),putchar('\n');
return 0;
}
后记
本题可以用min_25而不用杜教筛
但我太菜了没有学过非递归min_25
QAQ
2283
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
