51nod 1847 奇怪的数学题(min25)

51nod 1847 奇怪的数学题(min25)

http://www.51nod.com/Challenge/Problem.html#problemId=1847

设\(f(n)\)为\(n>1\)的次小约数(\(f(1)=0\))，很显然\(f(n)={n\over {\rm minp} (n)}\)，其中\({\rm minp}(n)\)表示n的最小质因数。

原式就变成

\[\sum_i^n\sum_j^n f(\gcd(i,j))^k \]

先枚举一下\(\gcd\)就变成了

\[\sum_{d=2} f(d)^k \sum_i^{[n/d]}\sum_j^{[n/d]} [i\perp j]=\sum_{d=2} f(d)^k \big(-1+2\sum_i^{[n/d]}\varphi(i)\big) \]

前后两个部分，后面那个部分是板子，主要看前面\(f(d)^k\)的求法

考虑下min25筛的tao核lu心式子

\[g(n,j)=g(n,j-1)-p^k\big(g(\lfloor {n\over p}\rfloor,j-1)-{\rm sump} (j-1)\big) \]

总数-\({\rm minp} (n)=p_j\)的所有合数的\(k\)次和

因此我们在转移的时候记录好\(\big(g(\lfloor {n\over p}\rfloor,j-1)-{\rm sump} (j-1)\big)\)的和就行了

还有个小问题是\(f(p)=1\)，我们没有在上面统计出质数的贡献，顺便处理下质数个数和就行。

//@winlere
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<cmath>

using namespace std;  typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int uint;
inline int qr(){
	int ret=0,f=0,c=getchar();
	while(!isdigit(c)) f|=c==45,c=getchar();
	while( isdigit(c)) ret=ret*10+c-48,c=getchar();
	return f?-ret:ret;
}
const int maxn=1e6+5;
int n,k,N,pr[maxn],val[maxn],cnt,id;
uint phi[maxn],sphi[maxn];
uint sumpk[maxn],sump0[maxn];
uint s0[maxn],sk[maxn],str[51][51],tot[maxn];

uint ksm(uint base,int p){
	uint ret=1;
	while(p){
		if(p&1) ret*=base;
		base*=base; p>>=1;
	}
	return ret;
}

void seive(int n){
	static bool usd[maxn]; phi[1]=1; sphi[1]=1;
	for(int t=2;t<=n;++t){
		if(!usd[t]) pr[++cnt]=t,sumpk[cnt]=sumpk[cnt-1]+ksm(t,k),sump0[cnt]=sump0[cnt-1]+1u,phi[t]=t-1u;
		for(int i=1;i<=cnt;++i){
			int v=pr[i]*t;
			if(v>n) break;
			usd[v]=1; 
			if(t%pr[i]) phi[v]=phi[t]*phi[pr[i]];
			else {phi[v]=phi[t]*pr[i];break;}
		}
		sphi[t]=sphi[t-1]+phi[t];
	}
}

int&getId(int val){
	static int id[maxn],id2[maxn];
	return val<=N?id[val]:id2[n/val];
}

uint sum0(ull n){return n;}
uint sumk(ull n){
	uint ret=0;
	for(uint t=0,ed=min<uint>(k,n);t<=ed;++t){
		uint mi=1;
		for(uint i=0,f=0;i<=t;++i)
			if(f||(n+1u-i)%(t+1)) mi*=n+1u-i;
			else mi*=(n+1u-i)/(t+1),f=1;
		ret+=str[k][t]*mi;
	}
	return ret;
}

void init(){
	str[0][0]=1;
	for(int t=1;t<=50;++t)
		for(int i=1;i<=t;++i)
			str[t][i]=str[t-1][i-1]+str[t-1][i]*i;
	for(int l=1,r=1;r<=n;l=++r){
		r=n/(n/l); getId(n/l)=++id;
		val[id]=n/l; s0[id]=sum0(n/l)-1; sk[id]=sumk(n/l)-1;
		//cerr<<"insert::"<<n/l<<endl;
	}
	for(int t=1;t<=cnt&&1ll*pr[t]*pr[t]<=val[1];++t){
		uint w=ksm(pr[t],k);
		for(int i=1;i<=id&&pr[t]*pr[t]<=val[i];++i){
			int k=getId(val[i]/pr[t]);
			sk[i]-=w*(sk[k]-sumpk[t-1]);
			s0[i]-=  (s0[k]-sump0[t-1]);
			tot[i]+=sk[k]-sumpk[t-1];
		}
	}
	for(int t=1;t<=id;++t) tot[t]+=s0[t];//,cerr<<val[t]<<" :=: "<<s0[t]<<endl;
}

uint que(uint x){
	return tot[getId(x)];
}

map<int,uint> rem;
uint S(int n){
	if(n<=1000000) return sphi[n];
	if(rem.find(n)!=rem.end()) return rem[n];
	uint ret=1ull*n*(n+1u)>>1;
	for(int l=2,r=2;r<=n;l=++r)
		r=n/(n/l),ret-=(r-l+1u)*S(n/l);
	return rem[n]=ret;
} 

int main(){
	n=qr(); k=qr(); N=sqrt(n);
	seive(1e6);
	init();
	uint ans=0;
	for(int l=1,r=1;r<=n;l=++r)
		r=n/(n/l),ans+=(que(r)-que(l-1))*(2u*S(n/l)-1);
	printf("%u\n",ans);
	return 0;
}