【51nod 1874】 奇怪的数学题

题目

求

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

首先这个次大公约数显然就是\(gcd\)除一下最小质因子了

于是

\[\sum_{i=1}^n\sum_{j=1}^n(\frac{(i,j)}{minp((i,j))})^k \]

显然可以枚举\(gcd\)，之后利用欧拉函数来求和

\[\sum_{d=1}^n(\frac{d}{minp(d)})^k(2\sum_{i=1}^{\frac{n}{d}}\varphi(i)-1) \]

后面那个东西显然可以直接杜教筛了，但是我们想要整除分块的话就需要前面这个函数的前缀和

考虑到最小质因子是一个跟埃氏筛非常有关系的东西，所以我们得知这题的正解应该是min_25

考虑我们用min_25去筛一个函数叫\(f(i)=i^k\)，我们只需要在min_25的过程中把每一个质数作为最小质因子的那些数的\(k\)次方和筛出来就好了

显然我们这样并没有计算\(d\)为质数的情况，于是我们顺便筛一个质数个数就好了

现在的问题变成了如何求自然数得幂次方和这样一个经典的问题，显然我们随便插个值就好了

但是就会发现模数不是质数，需要求逆元的方法就都凉了

有这样一个非常牛逼的做法是不需要逆元的，就是斯特林数

\[\sum_{i=1}^ni^k=\sum_{i=1}^kS_2(k,i)\frac{(n+1)^{\underline {(j+1)}}}{j+1} \]

不想推了，就抄yyb的吧

至于上面那个柿子为明明有分母却不需要逆元是因为\(n+1\)
的\(j+1\)次下降幂一定会有一个数是\(j+1\)的倍数，只要我们乘这个数得时候乘上这个数除以\(j+1\)的商就好了

代码

#include<cmath>
#include<cstdio>
#define re register
#define uint unsigned int
const int maxn=1e5+5;
int n,k,U,Sqr,o,m;
uint phi[1500005],g[maxn],S[2][55],s[maxn],pre[maxn],sum[maxn],h[maxn];
int f[1500005],p[500005],w[maxn],id1[maxn],id2[maxn];
inline int chk(int x) {return x<=Sqr?id1[x]:id2[n/x];}
inline uint ksm(uint a,int b) {
	uint S=1;
	while(b) {if(b&1) S*=a;b>>=1;a*=a;}
	return S;
}
uint Sphi(int x) {
	if(x<=U) return phi[x];
	if(sum[n/x]) return sum[n/x];
	uint now=1ll*x*(x+1)/2ll;
	for(re int l=2,r;l<=x;l=r+1) {
		r=x/(x/l);
		now-=Sphi(x/l)*(r-l+1);
	}
	return sum[n/x]=now;
}
inline uint calc(int n) {
	uint cnt=0;
	for(re int i=1;i<=k;i++) {
		uint ans=1;
		int flag=0;
		for(re int t=i+1,j=n+1;t;--j,--t) 
			if(!flag&&j%(i+1)==0) flag=1,ans*=(j/(i+1));else ans*=j;  
		cnt+=S[o][i]*ans;
	}
	return cnt;
} 
int main() {
	scanf("%d%d",&n,&k);U=std::pow(n,0.666)+1;Sqr=std::sqrt(n)+1;
	f[1]=1,phi[1]=1;o=0;S[0][1]=1;
	for(re int i=2;i<=U;i++) {
		if(!f[i]) p[++p[0]]=i,phi[i]=i-1,pre[p[0]]=pre[p[0]-1]+ksm(i,k);
		for(re int j=1;j<=p[0]&&p[j]*i<=U;j++) {
			f[p[j]*i]=1;
			if(i%p[j]==0) {phi[p[j]*i]=p[j]*phi[i];break;}
			phi[p[j]*i]=(p[j]-1)*phi[i];
		}
	}
	for(re int i=2;i<=k;i++,o^=1)
		for(re int j=1;j<=i;j++)
		 	S[o^1][j]=S[o][j-1]+S[o][j]*(uint)j;
	for(re int i=1;i<=U;i++) phi[i]+=phi[i-1];
	for(re int l=1,r;l<=n;l=r+1) {
		r=n/(n/l);w[++m]=n/l;
		if(w[m]<=Sqr) id1[w[m]]=m;else id2[n/w[m]]=m;
		g[m]=calc(w[m])-1;h[m]=w[m]-1;
	}
	for(re int j=1;j<=p[0]&&1ll*p[j]*p[j]<=n;j++) {
		uint tmp=ksm(p[j],k);
		for(re int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++) {
			int t=chk(w[i]/p[j]);
			uint tot=tmp*(g[t]-pre[j-1]);
			s[i]+=g[t]-pre[j-1];g[i]-=tot;
			h[i]-=(h[t]-j+1);
		}
	}
	for(re int i=1;i<=m;i++) s[i]+=h[i];
	uint ans=0;
	for(re int l=1,r;l<=n;l=r+1) {
		r=n/(n/l);
		ans+=(2ll*Sphi(n/l)-1)*(s[chk(r)]-s[chk(l-1)]);
	}
	printf("%u\n",ans);
	return 0;
}