【51Nod1847】奇怪的数学题

​ 　　
　　
　　记\(f(x)=\)\(x\)的次大因数，那么\(sgcd(i,j)=f(gcd(i,j))\)。
　　
　　下面来推式子：

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^nsgcd(i,j)^k&=\sum_{i=1}^n\sum_{j=1}^nf(gcd(i,j))^k&记d=gcd(i,j)\\ &=\sum_{d=1}^nf(d)^k\sum_{\frac i d=1}^{\lfloor \frac n d\rfloor}\sum_{\frac j d=1}^{\lfloor \frac n d\rfloor}[gcd(\frac i d,\frac j d)=1]&f(1)=0\\ &=\sum_{d=2}^nf(d)^k\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac n d\rfloor}[gcd(i,j)=1]\\ &=\sum_{d=2}^nf(d)^k(2\sum_{i=1}^{\lfloor \frac n d\rfloor}\varphi(i)-1) \end{aligned} \]

　　最后一步的括号是用欧拉函数的定义直接替换出来的。
　　
​　　我们发现，可以按\(\lfloor \frac n d \rfloor\)的取值数论分段，因为括号显然只受\(\lfloor \frac n d\rfloor\)的取值影响，关键是如何求\(f(x)^k\)的前缀和？
　　
​　　记\(S_x=\sum_{i=2}^xf(x)^k\)，
　　
​　　令min_25中的\(g\)数组以\(s(x)=x^k\)的定义计算\(g\)，考虑由\(g_{i,j-1}\)推到\(g_{i,j}\)的时候，减去的是什么？是\(s(p_j)*(g_{\lfloor\frac n {p_j}\rfloor,j-1}-g_{p_j-1,j-1})\)，后面括号的部分是什么呢？恰好是那些最小质因子为\(p_j\)且除去\(p_j\)后剩余部分最小质因数不小于\(p_j\)的数的\(k\)次方之和，我们发现这些数的\(f\)值之和就是后面的括号。又因为每一个数至多被如此枚举到1次，所以对于每一个\(i\)，我们把括号的值累加起来，这就是那些\([2,i]\)中非质数的\(f\)值之和；再加上\([2,i]\)中的质数的\(f\)值之和，也就是质数个数，我们就求得了\(S_i\)。
　　
　　min_25真牛逼。
　　
　　所以就做完了，\(g\)的初始化需要用到自然数幂求和，考虑到\(k\)比较小，可以用第二类斯特林数求解。

\[\begin{aligned} Sum_k(n)&=\sum_{i=1}^ni^k\\ &=\sum_{j=1}^k\begin{Bmatrix}k\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1} \end{aligned} \]

　　
　　

Code

　　

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <map>
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef map<ll,ui> mlu;
const int SQRTN=32000,LEN=1000001;
bool vis[LEN];
ui p[LEN],pcnt,pphi[LEN],s[51][51],pk[LEN];
ll n,k,sqrtn,m;
ll a[SQRTN*2],cnt,pos1[SQRTN],pos2[SQRTN];
ui g[SQRTN*2],sum[SQRTN*2],g1[SQRTN*2];
mlu record;
ui ksm(ui x,int y){
	ui res=1;
	for(;y;x=x*x,y>>=1)
		if(y&1) res=res*x;
	return res;
}
void prework(){
	pphi[1]=1;
	for(int i=2;i<LEN;i++){
		if(!vis[i]){
			p[++pcnt]=i;	
			pphi[i]=i-1;
		}
		for(int j=1;j<=pcnt&&i*p[j]<LEN;j++){
			int x=i*p[j];
			vis[x]=true;
			if(i%p[j]==0){
				pphi[x]=pphi[i]*p[j];
				break;
			}
			pphi[x]=pphi[i]*(p[j]-1);
		}
	}
	for(int i=2;i<LEN;i++) pphi[i]+=pphi[i-1];
	s[0][0]=s[0][1]=1;
	for(int i=1;i<=50;i++){
		s[i][1]=1;
		for(int j=2;j<=i;j++)
			s[i][j]=s[i-1][j]*(ui)j+s[i-1][j-1];
	}
}
inline int gp(ll x){
	return x<=sqrtn?pos1[x]:pos2[n/x];
}
void Diz(){
	for(ll i=1,j;i<=n;i=j+1){
		j=n/(n/i);
		a[++cnt]=n/i;
	}
	reverse(a+1,a+1+cnt);
	for(int i=1;i<=cnt;i++)
		if(a[i]<=sqrtn) pos1[a[i]]=i;
		else pos2[n/a[i]]=i;
}
ui sumup(ll l,ll r){
	ll a=l+r,b=r-l+1;
	if(a&1) b/=2; else a/=2;
	return (ui)a*(ui)b;
}
ui getPhi(ll n){
	if(n<LEN) return pphi[n];
	mlu::iterator pos=record.find(n);
	if(pos!=record.end()) return pos->second;
	ui res=sumup(1,n);
	for(ll i=2,j;i<=n;i=j+1){
		j=n/(n/i);
		res-=getPhi(n/i)*(j-i+1);
	}
	record[n]=res;
	return res;
}
ui getSumk(ll n){
	ui res=0,t;
	for(ll j=1;j<=k;j++){
		t=1;	
		int md=(n-j+1)%(j+1);
		for(ll i=n-j+1;i<=n+1;i++,md++)
			if(!md||md==(j+1)) t*=(ui)i/(j+1);
			else t*=(ui)i;
		res+=t*s[k][j];
	}
	return res;
}
void calc_g(){
	for(int i=2;i<=cnt;i++) g[i]=getSumk(a[i])-1,g1[i]=a[i]-1,sum[i]=0;
	for(int j=1;j<=m;j++)
		for(int i=cnt;i>=2&&a[i]>=p[j]*p[j];i--){
			g1[i]-=g1[gp(a[i]/p[j])]-g1[gp(p[j]-1)];
			ui delta=(g[gp(a[i]/p[j])]-g[gp(p[j]-1)]);
			g[i]-=pk[j]*delta;
			sum[i]+=delta;
		}
}
ui calc(ll n){
	return sum[gp(n)]+g1[gp(n)];
}
int main
	prework();
	scanf("%lld%lld",&n,&k);
	sqrtn=(ll)sqrt(n);
	m=upper_bound(p+1,p+1+pcnt,sqrtn)-p-1;
	for(int i=1;i<=m;i++) pk[i]=ksm(p[i],k);
	Diz();
	calc_g();
	ui ans=0,last=0,tmp;
	for(ll i=2,j;i<=n;i=j+1){
		j=n/(n/i);
		tmp=calc(j);
		ans+=(getPhi(n/i)*2-1)*(tmp-last);
		last=tmp;
	}
	printf("%lu\n",ans);
	return 0;
}