【51NOD 1847】奇怪的數學題（莫比烏斯反演，杜教篩，min_25篩，第二類斯特林數）
題面
\[\sum_{i=1}^n\sum_{j=1}^nsgcd(i,j)^k \]
其中\(sgcd\)表示次大公約數。
題解
明擺着\(sgcd\)就是在\(gcd\)的基礎上除掉\(gcd\)的最小因數。
所以直接枚舉\(gcd\)。
\[\begin{aligned} ans&=\sum_{i=1}^n\sum_{j=1}^n sgcd(i,j)^k\\ &=\sum_{i=1}^n\sum_{j=1}^n (\frac{gcd(i,j)}{min_p(gcd(i,j))})^k\\ &=\sum_{d=1}^n (\frac{d}{min_p(d)})^k\sum_{i=1}^n\sum_{j=1}^n [gcd(i,j)=1]\\ &=\sum_{d=1}^n (\frac{d}{min_p(d)})^k\sum_{i=1}^{n/d}\mu(i)[\frac{n}{id}]^2\\ &=\sum_{T=1}^n [\frac{n}{T}]^2\sum_{d|T}(\frac{d}{min_p(d)})^k\mu(\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}^k\begin{Bmatrix}k\\j\end{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\n",ans);
return 0;
}