51nod1847 奇怪的数学题

题目大意： 求

\[\sum _i^N\sum_j^Nsgcd(i,j)^k \]

\(N \le 1e9 ,k \le 50\)

先推式子：

设\(x\)的最大不等于x的约数是\(f(x)​\)

\[\begin{aligned} &\sum _i^N\sum_j^Nsgcd(i,j)^k \\ =&\sum _d^N \sum _i^{\lfloor\frac nd\rfloor}\sum_j^{\lfloor\frac nd\rfloor} [gcd(i,j)=1]f(d)^k \\ =&\sum _d^N (2\sum _i^{\lfloor\frac nd\rfloor}\varphi(i) -1) * f(d)^k \end{aligned} \]

前面可以杜教筛，现在关注如何求\(\sum f(d)^k\)

考虑min25筛，min25筛容斥的一步正好求的是最小质因子是\(prime_i\)的一堆值。将这一步减去的全部加到答案上，就行了。具体如下

\[g(i,0)=\sum _{i=2}^N i^k \\ g(i,j)=g(i,j-1)-(g(\lfloor\frac i {prime_j}\rfloor,j-1)-\sum _{p=1}^{j-1}prime_p^k)*(prime_j^k) \\ h(i)+=(g(\lfloor\frac i {prime_j}\rfloor,j-1)-\sum _{p=1}^{j-1}prime_p^k) \]

需要注意的是\(\sum _{i=2}^N i^k\)怎么求：

\[\begin{aligned} &\sum _i^n i^k \\ =&\sum _i^n \sum _{j=0}^k \begin{Bmatrix} k \\ j \end{Bmatrix}i^\underline{j} \\ =&\sum _{j=0}^k \begin{Bmatrix} k \\ j \end{Bmatrix} \sum _i^n i^\underline{j} \\ =&\sum _{j=0}^k \begin{Bmatrix} k \\ j \end{Bmatrix} \frac {(n+1)^\underline{j+1}} {j+1} \end{aligned} \]

代码

#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/hash_policy.hpp>
using namespace __gnu_pbds;
using namespace std;
typedef long long ll;
typedef unsigned int uint;
const int N=1e6+5;
int n,prime[N],pcnt,v[N],phi[N],k;
uint pre[N],preprime[N];

inline void sieve(int n=1000000){
	phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!v[i])prime[++pcnt]=i,phi[i]=i-1;
		for(int j=1;j<=pcnt&&1ll*prime[j]*i<=n;j++){
			int nxt=i*prime[j];v[nxt]=1;
			if(i%prime[j]) phi[nxt]=phi[i]*(prime[j]-1);
			else {
				phi[nxt]=phi[i]*prime[j];
				break;
			}
		}
	}
	for(int i=1;i<=n;i++)pre[i]=pre[i-1]+phi[i];
}
gp_hash_table<int,uint> Table;
uint calc(int x){
	if(x&1)return (uint)x*((x+1)>>1);
	else return (uint)(x>>1)*(x+1);
}
uint PHI(int n){
	if(n<=1000000)return pre[n];
	if(Table.find(n)!=Table.end())return Table[n];
	uint ans=calc(n);
	for(int u=2,v;u<=n;u=v+1){
		v=n/(n/u);
		ans-=(uint)(v-u+1)*PHI(n/u);
	}
	return Table[n]=ans;
}
uint S[100][100];
inline void stir_init(int n=60){
	S[0][0]=1;
	for(int i=1;i<=n;i++)for(int j=1;j<=i;j++){
		S[i][j]=S[i-1][j-1]+(uint)j*S[i-1][j];
	}
}

inline uint qs(int n,int m=k){
	uint ans=0;
	for(int j=0;j<=m;j++){
		uint cur=1;bool flag=0;
		for(int k=0;k<=j;k++){
			if(!flag && (n+1-k)%(j+1)==0)flag=1,cur*=(uint)((n+1-k)/(j+1));
			else cur*=(uint)(n+1-k);
		}
		ans+=cur*S[m][j];
	}
	return ans;
}

int Sqr;
int id1[N],id2[N],w[N],m;
uint f[N],g[N],Count[N];

uint qpow(uint a,int b){uint ret=1;for(;b;b>>=1,a=a*a)if(b&1)ret=ret*a;return ret;}

inline void init(){
	Sqr=sqrt(n);
	for(int u=1,v;u<=n;u=v+1){
		v=n/(n/u);
		w[++m]=n/u;
		if(w[m]<=Sqr)id1[w[m]]=m;
		else id2[n/w[m]]=m;
		f[m]=qs(n/u)-1;
		Count[m]=n/u-1;
	}
	for(int j=1;1ll*prime[j]*prime[j]<=n;j++){
		uint Pw=qpow(prime[j],k);preprime[j]=preprime[j-1]+Pw;
		for(int i=1;i<=m && 1ll*prime[j]*prime[j] <= w[i];i++){
			int idd=(w[i]/prime[j])<=Sqr?id1[w[i]/prime[j]]:id2[n/(w[i]/prime[j])];
			f[i]-=Pw*(f[idd]-preprime[j-1]);
			g[i]+=f[idd]-preprime[j-1];
			Count[i]-=(Count[idd]-j+1);
		}
	}
}
int getid(int x){
	if(x==0)return 0;
	return x<=Sqr?id1[x]:id2[n/x];
}

int main()
{
	sieve();stir_init();
	cin >> n >> k;
	init();
	uint ans=0;
	for(int u=1,v;u<=n;u=v+1){
		v=n/(n/u);
		
		int id=getid(v),pst=getid(u-1);
		ans+=((g[id]+Count[id]) - (g[pst]+Count[pst]))*((uint)2*PHI(n/u)-1);
	}
	cout << ans << endl;
}