【51NOD 1847】奇怪的数学题（莫比乌斯反演，杜教筛，min_25筛，第二类斯特林数）

【51NOD 1847】奇怪的数学题（莫比乌斯反演，杜教筛，min_25筛，第二类斯特林数）

题面

51NOD

[sum_{i=1}^nsum_{j=1}^nsgcd(i,j)^k ]

其中(sgcd)表示次大公约数。

题解

明摆着(sgcd)就是在(gcd)的基础上除掉(gcd)的最小因数。
所以直接枚举(gcd)。

[egin{aligned} ans&=sum_{i=1}^nsum_{j=1}^n sgcd(i,j)^k\ &=sum_{i=1}^nsum_{j=1}^n (frac{gcd(i,j)}{min_p(gcd(i,j))})^k\ &=sum_{d=1}^n (frac{d}{min_p(d)})^ksum_{i=1}^nsum_{j=1}^n [gcd(i,j)=1]\ &=sum_{d=1}^n (frac{d}{min_p(d)})^ksum_{i=1}^{n/d}mu(i)[frac{n}{id}]^2\ &=sum_{T=1}^n [frac{n}{T}]^2sum_{d|T}(frac{d}{min_p(d)})^kmu(frac{T}{d}) end{aligned}]

好啦好啦，又有一个看起来要求前缀和的狄利克雷卷积啦。
令(displaystyle f(d)=(frac{d}{min_p(d)})^k)
那么后面那个东西就是((f*mu)(T))
要求它的前缀和，emmmm，杜教筛。
看到(mu)，emmm，构造(g(x)=1)
那么令(displaystyle S(n)=sum_{i=1}^n (f*mu)(i))
那么套杜教筛的式子：

[g(1)S(n)=sum_{i=1}^n (f*mu*g)(i)-sum_{i=2}^n g(i)S([frac{n}{i}]) ]

把(g(x)=1)带进去就得到了：

[S(n)=sum_{i=1}^n f(i)-sum_{i=2}^nS([frac{n}{i}]) ]

行，来筛(f)前缀和，涉及到了最小质因子？？？(min\_25)筛。
考虑一下(f)的前缀和是什么东西，首先所有质数的贡献都是(1)，这里需要计算质数个数。
然后就是枚举最小质因子求剩下部分的(k)次幂的和就行了。
而很有趣的一点就是(f)所要求的所有合数的贡献，恰好就是计算所有质数(k)次方贡献时减去的部分。所以直接加上就好啦。
然后这里涉及到了自然数幂和的问题。推导戳这里

[sum_{i=1}^n i^k=sum_{j=0}^kegin{Bmatrix}k\jend{Bmatrix}frac{(n+1)^{underline {j+1}}}{j+1} ]

这样一来就可以直接算啦。

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define uint unsigned int
#define MAX 100000
uint n,blk;int k;
bool zs[MAX];
int pri[MAX],tot;
uint prik[MAX],sprik[MAX];
uint w[MAX],g[MAX],h[MAX],f[MAX];
int id1[MAX],id2[MAX],m;
uint S[60][60];
int getid(uint x){return (x<=blk)?id1[x]:id2[n/x];}
uint fpow(uint a,int b)
{
	uint s=1;
	while(b){if(b&1)s*=a;a*=a;b>>=1;}
	return s;
}
void pre(int n)
{
	for(int i=2;i<=n;++i)
	{
		if(!zs[i])pri[++tot]=i,prik[tot]=fpow(i,k);
		for(int j=1;j<=tot&&i*pri[j]<=n;++j)
		{
			zs[i*pri[j]]=true;
			if(i%pri[j]==0)break;
		}
	}
	for(int i=1;i<=tot;++i)sprik[i]=sprik[i-1]+prik[i];
}
uint Sumk(int n)
{
	uint ret=0;
	for(int i=1;i<=k;++i)
	{
		uint s=1;
		for(int j=0;j<=i;++j)
			if((n+1-j)%(i+1))s*=n+1-j;
			else s*=(n+1-j)/(i+1);
		ret+=S[k][i]*s;
	}
	return ret;
}
bool vis[MAX];
uint M[MAX];
uint Solve(uint n)
{
	if(vis[getid(n)])return M[getid(n)];
	uint ret=f[getid(n)];
	for(uint i=2,j;i<=n;i=j+1)
		j=n/(n/i),ret-=(j-i+1)*Solve(n/i);
	vis[getid(n)]=true;return M[getid(n)]=ret;
}
int main()
{
	scanf("%u%d",&n,&k);pre(blk=sqrt(n));S[0][0]=1;
	for(int i=1;i<=k;++i)
		for(int j=1;j<=i;++j)
			S[i][j]=S[i-1][j-1]+j*S[i-1][j];
	for(uint i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);w[++m]=n/i;
		g[m]=w[m]-1;h[m]=Sumk(w[m])-1;
		if(w[m]<=blk)id1[w[m]]=m;
		else id2[n/w[m]]=m;
	}
	for(int j=1;j<=tot;++j)
		for(int i=1;i<=m&&1u*pri[j]*pri[j]<=w[i];++i)
		{
			int k=getid(w[i]/pri[j]);
			g[i]-=g[k]-(j-1);
			h[i]-=prik[j]*(h[k]-sprik[j-1]);
			f[i]+=h[k]-sprik[j-1];
		}
	for(int i=1;i<=m;++i)f[i]+=g[i];
	uint ans=0;
	for(uint i=1,j,lt=0,nw;i<=n;i=j+1)
	{
		j=n/(n/i);uint s=n/i;s*=s;
		nw=Solve(j);ans+=s*(nw-lt);lt=nw;
	}
	printf("%u
",ans);
	return 0;
}