min_25筛

Posted 2021-02-18 sktt1faker

算法简介

min_25筛是min_25发明的能在(O(frac{n^{frac 3 4}}{log_2n}))复杂度内解决一类积性函数前缀和的算法，这类函数在自变量为质数的取值必须为多个完全积性函数的线性组合。

理论

约定(minp(i))代表(i)的最小质因子，(p_j)为从小到大第(i)个质数，(f^*(i))为一个在自变量取质数值时与(f(i))的值完全相同的完全积性函数，(P)为全体质数集合，(|P|)为集合(P)的大小。
首先得求出来(sumlimits_{i=1}^n{[i∈P]f(i)})。
先定义(g(n,j)=sumlimits_{i=1}^n{[i∈P or minp(i)>p_j]f^*(i)})，这相当于把所有(p_j)的倍数筛掉了。现在我们需要求出这个东西。
我们考虑用(g(n,j-1))转移(g(n,j))，这相当于在(g(n,j))基础上再筛掉(p_j)为最小质因子的数，那么分情况：

1. (p_j^2>n)，以(p_j)为最小质因子的最小数都大于(n)了，你筛什么筛？于是(g(n,j)=g(n,j-1))；
2. (p_j^2<=n)，那么我们要筛去所有最小质因子为(p_j)的数。根据完全积性函数的性质，我们先把(p_j)提出来，那么后面的那一部分就是(g(frac n {p_j},j-1)-g(p_j-1,j-1))，这代表从所有有(p_j)这个质因子的数中减去最小质因子小于(p_j)的数，于是(g(n,j)=g(n,j-1)-f^*(p_j)(g(frac n {p_j},j-1)-g(p_j-1,j-1)))。
   终于，我们推出了这个东西的转移：

[g(n,j)= egin{cases}g(n,j-1),p_j^2>ng(n,j-1)-f^*(p_j)(g(frac n {p_j},j-1)-g(p_j-1,j-1)),p_j^2<=nend{cases} ]

那么(sumlimits_{i=1}^n{[i∈P]f(i)}=g(n,|P|))。

接下来我们要用这个东西来推导所求答案。
令(S(n,j)=sumlimits_{i=1}^n{[minp(i)ge p_j]f(i)})，即所有最小质因子大于等于(p_j)的数的(f)的和。
我们把(S(n,j))分为(2)个部分分别计算：(i)为质数、(i)为合数，当然(1)你要单独算。

1. 对于所有的质数，当然就是从(g(n,j))中剔除那些小于(p_j)的质数，也就是(g(n,j)-sumlimits_{i=1}^{j-1}{p_j})；
2. 对于所有的合数，我们枚举最小质因子及其次数并把它们提出来，计算剩下的部分，最后加上我们处理不了的(p_k)整数次幂的情况，就是(sumlimits_{kge j}sumlimits_{e=1}^{p_k^{e+1}< n}(f(p_k^{e})S(frac n {p_k^e},k+1)+f(p_k^{e+1})))。
   因此，我们推出

[S(n,j)=g(n,|P|)-sumlimits_{i=1}^{j-1}{p_i}+ sumlimits_{kge j}^{p_k^2le n}sumlimits_{e=1}^{p_k^{e+1}le n}(f(p_k^{e})S(frac n {p_k^e},k+1)+f(p_k^{e+1})) ]

而我们要求的答案，不就是(S(n,1)+f(1))吗？

实现

这个东西实现还是有点小难。
我们在计算(g)的时候，可以把第二维滚掉。
对于(g)和(S)第一维，我们直接存存不下，但我们发现这个东西只有(lfloor frac n i floor)会被用到，于是可以来一手整除分块把(1)~(n)强行映射到(1)~(sqrt n)上：开一个(w[])存所有整除分块的值(lfloor frac n i floor)，而对于每个(lfloor frac n i floor)，如果它(<sqrt n)，那么在第一个数组中的(lfloor frac n i floor)位置存(w)的下标，否则在另一个数组中的(lfloor frac n {lfloor frac n i floor} floor)位置存(w)的下标。这样很抽象，我也没有理解得很清楚，但OI不就是一门不求甚解（？）的学科吗？背下来就完事了。
求(S)时直接递归，不需要记忆化，我也不知道为什么。
洛谷模板

#include<bits/stdc++.h>
#define int long long
typedef long long ll;
const int mod=1e9+7;
int n,sqr,tot,pri[100010],np[100010],sp1[100010],sp2[100010];
int sw,g1[200010],g2[200010],id1[100010],id2[100010];
ll w[200010];
inline ll Read()
{
	ll x=0,f=1;char ch=getchar();
	while(!isdigit(ch)&&ch!=‘-‘) ch=getchar();
	ch==‘-‘?f=-1:x=ch-‘0‘;
	while(isdigit(ch=getchar())) x=x*10+ch-‘0‘;
	return x*f;
}
inline void sieve()
{
	np[1]=1;
	for(register int i=2;i<=sqr;++i)
	{
		if(!np[i]) 
		{
			pri[++tot]=i;
			sp1[tot]=(sp1[tot-1]+i)%mod,sp2[tot]=(sp2[tot-1]+i*i)%mod;
		}
		for(register int j=1;j<=tot&&i*pri[j]<=sqr;++j)
		{
			np[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}
ll S(ll x,int y)
{
	if(pri[y]>=x) return 0;
	int k=(x<=sqr?id1[x]:id2[n/x]);
	ll res=(((ll)g2[k]-g1[k]-sp2[y]+sp1[y])%mod+mod)%mod;
	for(register int i=y+1;i<=tot&&(ll)pri[i]*pri[i]<=x;++i)
	{
		ll pe=pri[i];
		for(register int e=1;pe<=x;++e,pe*=pri[i])
		{
			int xx=pe%mod;
			res=(res+(ll)xx*(xx-1)%mod*(S(x/pe,i)+(e!=1)))%mod;
		}
	}
	return res;
}
signed main()
{
	n=Read();
	sqr=sqrt(n);
	sieve();
	for(register ll l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		w[++sw]=n/r;
		g1[sw]=w[sw]%mod;
		g2[sw]=((ll)g1[sw]*(g1[sw]+1)%mod*(2*g1[sw]+1)%mod*166666668%mod-1)%mod;
		g1[sw]=((ll)g1[sw]*(g1[sw]+1)%mod*500000004%mod-1)%mod;
		if(w[sw]<=sqr) id1[w[sw]]=sw;
		else id2[r]=sw;
	}
	for(register int i=1;i<=tot;++i)
	{
		ll p2=(ll)pri[i]*pri[i];
		for(register int j=1;j<=sw&&p2<=w[j];++j)
		{
			ll k=w[j]/pri[i];
			k=(k<=sqr?id1[k]:id2[n/k]);
			g1[j]=(g1[j]-(ll)pri[i]*(g1[k]-sp1[i-1]+mod)%mod+mod)%mod;
			g2[j]=(g2[j]-p2%mod*(g2[k]-sp2[i-1]+mod)%mod+mod)%mod;
		}
	}
	printf("%lld",(S(n,0)+1)%mod);
	return 0;
}