51nod 1847 奇怪的数学题

Link
设\(f(n)\)表示\(n\)的次大因子，那么\(sgcd(i,j)=f(\gcd(i,j))\)。

\[\begin{aligned} ans&=\sum\limits_{i=1}^n\sum\limits_{j=1}^nf(\gcd(i,j))^k\\ &=\sum\limits_{i=1}^nf(i)^k\sum\limits_{u=1}^n\sum\limits_{v=1}^n[\gcd(i,j)=i]\\ &=\sum\limits_{i=1}^nf(i)^k\sum\limits_{u=1}^{\lfloor\frac ni\rfloor}\sum\limits_{v=1}^{\lfloor\frac ni\rfloor}\sum\limits_{d|\gcd(u,v)}\mu(d)\\ &=\sum\limits_{i=1}^nf(i)^k\sum\limits_{d=1}^{\lfloor\frac ni\rfloor}\mu(d)\lfloor\frac n{id}\rfloor^2 \end{aligned} \]

设\(g=f^k\times\mu\)（幂在点乘意义下），那么

\[\begin{aligned} ans&=\sum\limits_{i=1}^nf(i)^k\sum\limits_{d=1}^{\lfloor\frac nk\rfloor}\mu(d)\lfloor\frac n{id}\rfloor^2\\ &=\sum\limits_{i=1}^ng(i)\lfloor\frac ni\rfloor^2 \end{aligned} \]

对\(\lfloor\frac ni\rfloor^2\)整除分块，那么我们要计算的就是\(g\)在所有\(\lfloor\frac ni\rfloor\)的前缀和。
注意到\(g\times I=f\)，因此如果我们能够计算\(f^k\)在所有\(\lfloor\frac ni\rfloor\)的前缀和的前缀和，就能够用杜教筛计算\(\sum g\)了。

\[\sum\limits_{i=1}^nf(i)^k=\pi(n)+\sum\limits_{p\le\sqrt n}\sum\limits_{i=1}^{\lfloor\frac np\rfloor}[\min(i)\ge p]i^k \]

其中\(\min(i)\)表示\(i\)的最小非平凡质因子，若\(i=1\vee i\in\mathbb P\)则\(\min(i)=-\infty\)。
后面和Min_25筛计算\(\sum p^k\)时求\(h\)的过程中被筛掉的部分比较类似，模拟即可。
同时可以顺便求出需要用到的\(\pi(n)\)。
求\(h\)的过程中需要求自然数幂和，考虑利用Stirling数处理。

\[\begin{aligned} &\sum\limits_{i=1}^ni^k\\ =&\sum\limits_{i=1}^n\sum\limits_{j=0}^k\left\{k\atop j\right\}i^{\underline j}\\ =&\sum\limits_{j=0}^k\left\{k\atop j\right\}\frac{(n+1)^{\underline{j+1}}}{j+1} \end{aligned} \]

\(j+1\)不一定有逆元，注意到显然有\(j+1|(n+1)^{\underline{j+1}}\)，因此暴力计算下降幂时除掉即可。

#include<cstdio>
using u32=unsigned;
const int N=100007,K=57;
int n,k,tot,lim,cnt,pr[N],pos[N];u32 S[K][K],pw[N],sum_pw[N],sum[N];
int read(){int x;scanf("%d",&x);return x;}
int getid(int x){return x<=lim? x:cnt+1-n/x;}
u32 pow(u32 p,int e){u32 q=1;for(;e;e>>=1,p*=p)if(e&1)q*=p;return q;}
u32 power_sum(int n)
{
    u32 ans=0;
    for(int i=1;i<=k&&i<=n;++i)
    {
	u32 now=1;int f=0;
	for(int j=n+1;j>n-i;--j) if(!(j%(i+1))&&!f) now*=j/(i+1),f=1; else now*=j;
	ans+=now*S[k][i];
    }
    return ans;
}
void init()
{
    static int low[N];S[0][0]=1;
    for(lim=1;(lim+1)*(lim+1)<=n;++lim);
    for(int i=1;i<=k;++i) for(int j=1;j<=i;++j) S[i][j]=S[i-1][j-1]+S[i-1][j]*j;
    for(int i=2;i<=lim;++i)
    {
	if(!low[i]) pr[++tot]=low[i]=i,sum_pw[tot]=sum_pw[tot-1]+(pw[tot]=pow(i,k));
	for(int j=1;j<=tot&&i*pr[j]<=lim&&pr[j]<=low[i];++j) low[i*pr[j]]=pr[j];
    }
}
void Min_25_Sieve()
{
    static u32 pi[N],t[N];
    for(int l=1,r;l<=n;l=r+1) r=n/(n/l),pos[++cnt]=r,pi[cnt]=r-1,t[cnt]=power_sum(r)-1;
    for(int i=1;i<=tot&&pr[i]<=lim;++i) for(int j=cnt,p;j&&pr[i]*pr[i]<=pos[j];--j) p=getid(pos[j]/pr[i]),pi[j]-=pi[p]-(i-1),sum[j]+=t[p]-sum_pw[i-1],t[j]-=pw[i]*(t[p]-sum_pw[i-1]);
    for(int i=1;i<=cnt;++i) sum[i]+=pi[i];
}
void Xudyh_Sieve()
{
    for(int i=1;i<=cnt;++i) for(int l=2,r;l<=pos[i];l=r+1) r=pos[i]/(pos[i]/l),sum[i]-=(r-l+1)*sum[getid(pos[i]/l)];
}
void calc()
{
    u32 ans=0;
    for(int l=1,r,p;l<=n;l=r+1) p=getid(r=n/(n/l)),ans+=(sum[p]-sum[p-1])*(n/l)*(n/l);
    printf("%u",ans);
}
int main()
{
    n=read(),k=read();
    init(),Min_25_Sieve(),Xudyh_Sieve(),calc();
}