前言
万年不写公开博客了，这次填个坑
题目相关
题目大意
求(sum_{i=1}^Nsum_{j=1}^Nsgcd(i,j)^kpmod {2^{32}})
sgcd(i,j)为i,j的第二大公约数
特殊的，若gcd(i,j)=1,则sgcd(i,j)=0
数据范围
(nle10^9,kle50)
题解
设m(x)为x的次大因数，特殊的，m(1)=0
[egin{aligned}
sum_{i=1}^Nsum_{j=1}^Nsgcd(i,j)^k&=sum_{i=1}^Nsum_{j=1}^Nm(gcd(i,j))^k&根据定义\
&=sum_{g=1}^Nm(g)^ksum_{i=1}^{leftlfloorfrac{N}{g}
ight
floor}sum_{j=1}^{leftlfloorfrac{N}{g}
ight
floor}[gcd(i,j)==1]&提出gcd\
&=sum_{g=1}^Nm(g)^kleft(2left(sum_{i=1}^{leftlfloorfrac{N}{g}
ight
floor}varphi(i)
ight)-1
ight)&大小互换也可以所以要乘以2，并且(1,1)只能算一次所以减一\
&=sum_{t=1}^nsum_{leftlfloorfrac{N}{g}
ight
floor=t}m(g)^kleft(2left(sum_{i=1}^{t}varphi(i)
ight)-1
ight)\
end{aligned}]
我们发现，满足条件的t只有(mathcal O(sqrt n ))个，而最后一项可以通过杜教筛来做
容易发现，对于同一个t，其满足条件的g是连续的一串，那么只要对于所有满足条件的t求(sum_{g=1}^{maxg}m(g)^k)，这个用min_25筛算
因为本质上min25筛会枚举最小因数，所以把最小因数排除掉即可
在预处理上述运算的时候要求(sum_{g=1}^{n}g^k)
我们引入斯特林数
[x^n=sum_{i=0}^xegin{Bmatrix}n\iend{Bmatrix}x^{idownarrow}
]
证明可以点进那篇博客看
这样的话，求sigma就是
[egin{aligned}
sum_{x=1}^nx^k&=sum_{x=1}^nsum_{i=0}^xegin{Bmatrix}k\iend{Bmatrix}x^{idownarrow}\
&=sum_{i=0}^negin{Bmatrix}k\iend{Bmatrix}sum_{x=i}^nx^{idownarrow}\
&=sum_{i=0}^negin{Bmatrix}k\iend{Bmatrix}sum_{x=i}^ninom{x}{i}i!\
&=sum_{i=0}^ni!egin{Bmatrix}k\iend{Bmatrix}sum_{x=i}^ninom{x}{i}\
&=sum_{i=0}^ni!egin{Bmatrix}k\iend{Bmatrix}inom{n+1}{i+1}\
&=sum_{i=0}^ni!egin{Bmatrix}k\iend{Bmatrix}frac{(n+1)!}{(i+1)!(n-i)!}\
&=sum_{i=0}^negin{Bmatrix}k\iend{Bmatrix}frac{(n+1)!}{(i+1)(n-i)!}\
&=sum_{i=0}^negin{Bmatrix}k\iend{Bmatrix}frac{(n+1)^{i+1downarrow}}{i+1}\
end{aligned}]
此处复杂度(mathcal O(k^2))
总复杂度(mathcal O(frac{n^{frac34}}{log n }+sqrt nk^2+n^{frac23}))
#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
",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
",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
",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
",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('
');
return 0;
}
后记
本题可以用min_25而不用杜教筛
但我太菜了没有学过非递归min_25
QAQ