51nod 1847 奇怪的数学题 莫比乌斯反演+min_25筛+杜教筛

题目分析

莫比乌斯反演

所谓的$sgcd(i,j)$，就是$gcd(i,j)$除以其最小的一个质因子。我们记$g(x)=(\frac{x}{minpri(x)}{)}^K$，答案就是求${\sum }_{i=1}^n{\sum }_{j=1}^{n}g(gcd(i,j))$

既然涉及到gcd，就用莫比乌斯反演搞一搞：

${\sum }_{d=1}^{n}g(d){\sum }_{i=1}^n{\sum }_{j=1}^n[gcd(i,j)=d]$

${\sum }_{d=1}^{n}g(d){\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{t|i,t|j}\mu (t)$

${\sum }_{d=1}^{n}g(d){\sum }_{t=1}^{\lfloor \frac{n}{d}\rfloor }\mu (t)(\lfloor \frac{n}{dt}\rfloor {)}^2$

${\sum }_{T=1}^n(\lfloor \frac{n}{T}\rfloor {)}^2{\sum }_{d|T}g(d)\mu (\frac{T}{d})$

前面的$T$部分可以数论分块，所以关键是求后面那一块的前缀和。

杜教筛

求$(g\ast \mu )(i)$的前缀和，想到杜教筛。杜教筛中使用的辅助函数，就用$\mu$常用的辅助函数$I$即可（这个函数是每一位都为1）

所以设$S(n)$表示那个卷积函数前$n$项的前缀和，有：

$I(1)S(n)={\sum }_{i=1}^n(g\ast \mu \ast I)(i)-{\sum }_{i=2}^{n}I(i)S(\lfloor \frac{n}{i}\rfloor )$

也就是：

$S(n)={\sum }_{i=1}^{n}g(i)-{\sum }_{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )$

又可以数论分块，关键问题又变成了求$g(i)$的前缀和

min_25筛

嘛，其实我还是喜欢叫这个东西扩展埃氏筛（因为那个下划线在用markdown的时候有点不方便，雾）。

设$s(i)$表示前$i$个数中的质数个数。

$f(i)$表示前$i$个数中质数的$k$次方和。

$g(i)$表示前$i$个数中合数都去掉自己的最小质因子后留下的数字的$k$次方和（和前面的$g$意义不同了）。

前两个很简单，转移式子分别是对于每个质因子$p$去掉以其为最小质因子的合数的贡献：$s(i)-=s(\lfloor \frac{i}{p}\rfloor )-s(p-1)$和$f(i)-=p^k(f(\lfloor \frac{i}{p}\rfloor )-f(p-1))$

发现每次$f(i)$丢掉的$f(\lfloor \frac{i}{p}\rfloor )-f(p-1)$是一些应该加进$g(i)$中的贡献，所以每次把$g(i)$加上这一串。

那么原式子中求的$g$，就是现在求出来的这些东西中的$g(i)+s(i)$（$s$是算质数的贡献）

现在又双叒叕有个问题，扩埃筛执行时，需要处理前$i$个数的函数贡献和。也就是我们要求前$i$个数的$K$次方和。

第二类斯特林数

这个问题我知道！用伯努利数就可以了！ヾ(≧∇≦*)ゝ

(╯°Д°)╯……个鬼啦！模数没逆元！

那就用第二类斯特林数吧，${\sum }_{i=1}^n{i}^K={\sum }_{i=1}^K{S}_2(K,i)\frac{(n+1{)}^{\underline{i+1}}}{i+1}$

（如果你以前没听说过这个式子，请别把下降幂看成幂了）

代码

#include<bits/stdc++.h>
using namespace std;
typedef unsigned int uint;
const int N=31630;
int K;uint n,ans,sqn;int vis1[N],vis2[N];
uint s1[N],s2[N],f1[N],f2[N],g1[N],g2[N],S2[55][55],mp1[N],mp2[N];

uint ksm(uint x,uint y) {
	uint re=1;
	for(;y;y>>=1,x=x*x) if(y&1) re=re*x;
	return re;
}
void prework() {
	S2[0][0]=1;
	for(uint i=1;i<=K;++i)
		for(uint j=1;j<=i;++j)
			S2[i][j]=S2[i-1][j-1]+j*S2[i-1][j];
}
uint power_sum(uint x) {
	uint re=0;
	for(uint i=0;i<=K;++i) {
		uint kl=S2[K][i];
		for(uint j=0;j<=i;++j)
			if((x+1-j)%(i+1)) kl*=x+1-j;
			else kl*=(x+1-j)/(i+1);
		re+=kl;
	}
	return re;
}
void min25() {
	if(n<=1) return;
	sqn=sqrt(n);
	for(uint i=1;i<=sqn;++i) {
		s1[i]=i-1,s2[i]=n/i-1;
		f1[i]=power_sum(i)-1,f2[i]=power_sum(n/i)-1;
	}
	for(uint p=2;p<=sqn;++p) {
		if(s1[p]==s1[p-1]) continue;
		uint tmps=s1[p-1],tmpf=f1[p-1],pk=ksm(p,K);
		for(uint i=1;i<=sqn/p;++i) {
			s2[i]-=s2[i*p]-tmps;
			f2[i]-=pk*(f2[i*p]-tmpf);
			g2[i]+=f2[i*p]-tmpf;
		}
		for(uint i=sqn/p+1;i<=n/p/p&&i<=sqn;++i) {
			s2[i]-=s1[n/i/p]-tmps;
			f2[i]-=pk*(f1[n/i/p]-tmpf);
			g2[i]+=f1[n/i/p]-tmpf;
		}
		for(uint i=sqn;i>=p*p;--i) {
			s1[i]-=s1[i/p]-tmps;
			f1[i]-=pk*(f1[i/p]-tmpf);
			g1[i]+=f1[i/p]-tmpf;
		}
	}
}
uint du_sieve(uint x) {
	if(x<=sqn?vis1[x]:vis2[n/x]) return x<=sqn?mp1[x]:mp2[n/x];
	uint re=(x<=sqn?g1[x]+s1[x]:g2[n/x]+s2[n/x]);
	for(uint i=2,j;i<=x;i=j+1) j=x/(x/i),re-=(j-i+1)*du_sieve(x/i);
	(x<=sqn?vis1[x]:vis2[n/x])=1,(x<=sqn?mp1[x]:mp2[n/x])=re;
	return re;
}

int main()
{
	scanf("%u%d",&n,&K);
	prework(),min25();
	for(uint i=1,j,las=0,now;i<=n;i=j+1,las=now)
		j=n/(n/i),now=du_sieve(j),ans+=(n/i)*(n/i)*(now-las);
	printf("%u\n",ans);
	return 0;
}