Min_25 筛学习笔记

My Blogs

杜教筛是一种能在 \(\mathcal O(n^{\frac 2 3})\) 的时间复杂度内求积性函数前缀和的筛法。虽然复杂度比较优秀，但是被筛的积性函数需要满足特殊性质。

Min_25 筛由 Min_25 发明，相对更通用，其时间复杂度为 \(\mathcal O(\frac{n^{\frac 3 4}}{\log n})\)。

首先构造一个完全积性函数 \(f'(x)\) 满足其在质数处和 \(f(x)\) 取值相同，其余位置任意。下设 \(\mathbb{P}\) 为质数集合，\(p_j\) 为第 \(j\) 个质数。

设 \(g(n,j)\) 表示 \(\sum_{i\leq n}f'(i)[i\in \mathbb{P}\lor\min_{x|i,x\in\mathbb{P}}>p_j]\)。一开始 \(g(n,0)\) 就是 \(\sum_{i\leq n}f'(i)\)，考虑 \(j\) 变大了之后会减少那些数的贡献：

若 \(p_j^2>n\)，则显然 \(g(n,j)=g(n,j-1)\)。

否则，有贡献的是最小质因子恰为 \(p_j\) 的数。可以考虑整体除掉一个 \(p_j\)，因为满足完全积性，所以所有应当删去的树贡献就全部少乘了 \(f'(p_j)\)。

这样就变成了 \(g(\lfloor\frac n {p_j}\rfloor,j-1)\)。但是这样回多算 \([1,\lfloor\frac n {p_j}\rfloor]\) 内质数的贡献，需要再减去。

转移方程为：

\[g(n,j)=\begin{cases} g(n,j-1)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad p_j^2>n\\ g(n,j-1)-f'(p_j)(g(\lfloor\frac n {p_j}\rfloor,j-1)-\sum_{j\leq \lfloor\frac n {p_j}\rfloor,j\in\mathbb{P}}f'(j))\;\;\;\;\;\;p_j^2\leq n \end{cases} \]

这样最后得到的 \(g(n,\lvert\mathbb P\rvert)\) 就是 \(n\) 以内所有质数处的函数值。经过上述操作可以把错误的函数值舍去，最后剩下的是对的函数值。

预处理 \(g\) 之后设 \(S(n,j)\) 表示 \(\sum_{i\leq n}f(i)(\min_{p|i,p\in \mathbb P}p>p_j)\)。

首先对于 \(i\) 是质数的部分，答案就是 \(g(n,\lvert\mathbb P\rvert)\)，这样会多算 \(<j\) 的质数，再减去 \(\sum_{i\in\mathbb P,i<j} f(i)\)。

否则枚举其最小的因子 \(p\) 和次幂 \(e\)，则有：

\[S(n,j)=g(n,\lvert\mathbb P\rvert)-\sum_{i\in\mathbb P,i<j} f(i)+\sum_{i>j}\sum_{e\geq 1}^{p_i^e\leq n}f(p_i^e)(S(\lfloor\frac n {p_i^e}\rfloor,i)+[e>1]) \]

加入 \([e>1]\) 是因为质数处的点值已经被计算过，但是 \(p^k(k>1)\) 的点值还没有被计算。

边界条件是 \(S(n,j)=0(p_j>n)\)。但是这样预处理 \(g\) 的复杂度达到了 \(\mathcal O(n\log n)\)。

注意到 \(g(i,j)\) 有用当且仅当存在 \(k\) 使得 \(\lfloor\frac n k\rfloor=i\)，所以有用的点只有 \(\mathcal O(\sqrt n)\) 个。总复杂度是 \(\mathcal O(\frac{n^{\frac 3 4}}{\log n})\) 不会证。

	ll n,w[300010],g1[300010],g2[300010],sum1[300010],sum2[300010],pr[300010];
	int N,m,cnt,id1[300010],id2[300010],iv6;
	bool v[300010];
	inline void sieve()
	{
		for(int i=2;i<=N;++i)
		{
			if(!v[i])pr[++cnt]=i;
			for(int j=1;j<=cnt&&pr[j]*i<=N;++j)
			{v[pr[j]*i]=1;if(i%pr[j]==0)break;}
		}
		for(int i=1;i<=cnt;++i)
		{
			sum1[i]=Cadd(sum1[i-1],pr[i]);
			sum2[i]=(sum2[i-1]+1ll*pr[i]*pr[i])%MOD;
		}
	}
	inline int f1(ll x){x%=MOD;return x*(x+1)/2%MOD;}
	inline int f2(ll x){x%=MOD;return x*(x+1)%MOD*(2*x+1)%MOD*iv6%MOD;}
	inline int getid(ll x){return x<=N?id1[x]:id2[n/x];}
	ll S(ll x,int j)
	{
		if(pr[j]>x)return 0;
		ll ans=Cdel(Cdel(g2[getid(x)],g1[getid(x)]),Cdel(sum2[j],sum1[j]));
		for(int i=j+1;i<=cnt&&pr[i]*pr[j]<=x;++i)
		{
			for(ll e=1,sp=pr[i];sp<=x;sp*=pr[i],++e)
			ans=(ans+sp%MOD*(sp%MOD-1)%MOD*(S(x/sp,i)+(e>1)))%MOD;
		}
		return ans;
	}
	inline void mian()
	{
		read(n),N=sqrt(n),sieve(),iv6=power(6,MOD-2);
		for(ll l=1,r;l<=n;l=r+1)
		{
			r=n/(n/l),w[++m]=n/l;
			g1[m]=f1(w[m])-1,g2[m]=f2(w[m])-1;
			if(w[m]<=N)id1[w[m]]=m;else id2[n/w[m]]=m;
		}
		for(int i=1;i<=cnt;++i)
		{
			for(int j=1;j<=m&&pr[i]*pr[i]<=w[j];++j)
			{
				g1[j]=((g1[j]-1ll*pr[i]*(g1[getid(w[j]/pr[i])]-sum1[i-1]))%MOD+MOD)%MOD;
				g2[j]=((g2[j]-1ll*pr[i]*pr[i]%MOD*(g2[getid(w[j]/pr[i])]-sum2[i-1]))%MOD+MOD)%MOD;
			}
		}
		cout<<S(n,0)+1;
	}