【模板】线性筛 Siever

posted on 2021-11-19 15:41:01 | under 模板 | source
posted on 2021-06-12 09:53:53 | under 学术 | source
posted on 2022-07-30 14:50:21 | under 模板 | source

problem

\(\forall x\leq n\)，求出 \(isprime(x),\varphi(x),\mu(x)\)，要求线性。

\[\varphi(x)=x\cdot\prod_i \frac{p_i-1}{p_i},\text{where }x=\prod_i p_i^{c_i}. \]

\[\mu(x)=\begin{cases} 0,&(\exists c_i>1),\\ (-1)^k,&(k=\sum_i [c_i=1]). \end{cases}\quad\text{where }x=\prod_i p_i^{c_i}.\]

solution-bruteforce

\(isprime(x)\)：验证 \([2,\sqrt{x}]\) 内的数是否整除 \(x\)。

点击查看代码

bool isprime(LL x){
	if(x<=2) return x==2;
	for(LL i=2;i*i<=x;i++) if(x%i==0) return 0;
	return 1;
}

\(\varphi(x)\)：大力质因数分解。

点击查看代码

LL phi(LL x){
	LL res=x;
	for(LL i=2;i*i<=x;i++){
		if(x%i==0){
			res=res/i*(i-1);
			while(x%i==0) x/=i;
		}
	}
	if(x>1) res=res/x*(x-1);
	return res;
}

\(\mu(x)\)：分解质因数。

点击查看代码

LL mu(LL x){
	LL res=1;
	for(LL i=2;i*i<=x;i++){
		if(x%i==0){
			res=-res,x/=i;
			if(x%i==0) return 0;
		}
	}
	if(x>1) res=-res;
	return res;
}

还有埃氏筛，但是一般不会用埃氏筛单纯地筛质数。可以做区间筛（筛根号里的质数再筛区间）。

点击查看代码

bool isp[100010];
void sieve(int n){
	memset(isp,1,sizeof isp),isp[0]=isp[1]=0;
	for(int i=1;i<=n;i++){
		if(!isp[i]) continue;
		for(int j=i+i;j<=n;j+=i) isp[j]=0;
	}
}

solution-primes

线性筛，我们用 \(i\cdot prime_j\) 来产生一个合数，显然这个 \(prime_j\) 就是 \(i\cdot prime_j\) 的最小质因数，证明：

先看一次代码：

for(int i=1;i<=N;i++){
	if(isp[i]) pri[++cnt]=i;
	for(int j=1;j<=cnt&&i*pri[j]<=N;j++){
		isp[i*pri[j]]=0;
		if(i%pri[j]==0) break;
	}
}

先解释 if(i%pri[j]==0) break：
当 \(prime_j|i\) 时，说明 \(i\) 已经被筛过了（是个合数），而且是被 \(prime_j\) 筛掉的（\(prime\) 递增），下一个 \(i\cdot prime_{j+1}\) 肯定也被筛了（被 \(prime_j\) 筛的），所以马上 break，保证每个合数只筛一次（这里保证了复杂度的线性）。

我们保证了每个合数只被筛一次，那么我们担心的无非是这次筛的时候 \(prime_j\) 不是最小质因数，那就假设一个合数 \(a\) 是被 \(i\cdot prime_j\) 筛掉，而真正的最小质因数是 \(prime_x\)，显然 \(i\mid prime_x\)（要不然这个因数在哪里），又有 \(prime_x<prime_j\)，\(i\) 会被 \(prime_x\) break 掉，就不存在 \(a\) 被 \(i\cdot prime_j\) 筛掉这件事。

我们已经证明「每个合数只筛一次」和「筛的 \(i\cdot prime_j\) 不可能不是最小质因数」，结合，则原命题成立。

这是一个有用的小技巧，如果质因数分解的数不大，可以大力跳最小质因子，复杂度变成 \(O(\log n)\)。习题：P3538 [POI2012]OKR-A Horrible Poem。

solution-phi

\(\varphi\) 是积性函数：\(\varphi(ab)=\varphi(a)\varphi(b)\)（其中 \(a\perp b\)）。（注：\(a\perp b\) 表示 \(a,b\) 互质，\(\gcd(a,b)=1\)）

用 \(i\times pri_j\) 筛时，若 \(i\perp pri_j\) 则可以直接转移：\(\varphi(i\times pri_j)=\varphi(i)\varphi(pri_j)\)。

否则，由于 \(pri_j\) 是质数，且 \(pri_j<i\)，剩下一种情况为 \(pri_j|i\)，这时 \(i\) 中包含了 \(pri_j\) 这一质因子，在外面乘上。\(\varphi(i\times pri_j)=\varphi(i)pri_j\)。

Recall: $\varphi$ 的定义

\[\varphi(x)=\boxed{x}\cdot\prod_i \frac{p_i-1}{p_i},\text{where }x=\prod_i p_i^{c_i}. \]

边界：质数 \(p\)，则 \(\varphi(p)=p-1\) 显然。

solution-mu

\(\mu\) 是积性函数：\(\mu(ab)=\mu(a)\mu(b)\)（其中 \(a\perp b\)）。

用 \(i\times pri_j\) 筛时，若 \(i\perp pri_j\) 则可以直接转移：\(\mu(i\times pri_j)=\mu(i)\mu(pri_j)\)。

否则，由于 \(pri_j\) 是质数，且 \(pri_j<i\)，剩下一种情况为 \(pri_j|i\)，这时 \(i\) 中包含了 \(pri_j\) 这一质因子，乘起来相当于有两个 \(pri_j\)，那么 \(\mu(i\times pri_j)=0\)。

边界：质数 \(p\) 的 \(\mu(p)=-1\)。

其它的积性函数应该都可以的。

code

注意到在标记非质数时，被标记的数一定 \(>i\)，又注意到当前质数个数 \(\leq i\)，所以可以将 \(isprime\) 和 \(prime\) 混用。

phi[n] 指 \(\varphi(n)\)，mu[n] 指 \(\mu(n)\)。

primes (version 1)

template<int N> class LineSiever{
    private:
        bool isp[N+10];
        int n,p[N+10];
    public:
        LineSiever():n(0){
            memset(isp,1,sizeof isp);
            isp[0]=isp[1]=0;
            for(int i=2;i<=N;i++){
                if(isp[i]) p[++n]=i;
                for(int j=1;j<=n&&1LL*i*p[j]<=N;j++){
                    if(isp[i*p[j]]=0,i%p[j]==0) break;
                }
            }
        }
        operator int(){return n;}
        bool operator()(int x){return isp[x];}
        int operator[](int x){return p[x];}
};
LineSiever<100> p;
cout<<p<<endl;//100 内有多少素数
cout<<(p(2)?"Yes":"No")<<endl;//2 是不是素数
cout<<p[3]<<endl;//第 3 个素数（5）
cout<<LineSiever<(int)1e8>()<<endl;//简写求出 1e8 内素数个数

primes (version 2)

template<int N> struct siever{
	int pri[N+10],cnt;
	siever(){
		memset(pri,1,sizeof pri),pri[0]=pri[1]=0;
		for(int i=1;i<=N;i++){
			if(pri[i]) pri[++cnt]=i;
			for(int j=1;j<=cnt&&i*pri[j]<=N;j++){
				pri[i*pri[j]]=0;
				if(i%pri[j]==0) break;
			}
		} 
	}
    int operator[](int i){return pri[i];}
};

primes,$\varphi$

template<int N> struct siever{
	int pri[N+10],phi[N+10],cnt;
	bool vis[N+10];
	siever(){
		memset(vis,1,sizeof vis);
		vis[phi[1]=1]=0;
		for(int i=1;i<=N;i++){
			if(!vis[i]) phi[pri[++cnt]=i]=i-1;
			for(int j=1;j<=cnt&&i*pri[j]<=N;j++){
				if(vis[i*pri[j]]=1,i%pri[j]){
					phi[i*pri[j]]=phi[pri[j]]*phi[i];
				}else{
					phi[i*pri[j]]=pri[j]*phi[i];
					break;
				}
			}
		}
	}
};

primes,$\varphi,\mu$

template<int N> struct siever{
	int cnt,p[N+10],phi[N+10],mu[N+10];
	sieve():cnt(0){
		memset(p,phi[1]=mu[1]=1,sizeof p),p[1]=0;
		for(int i=1;i<=N;i++){
			if(p[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1;
			for(int j=1;j<=cnt&&p[j]*i<=N;j++){
				p[p[j]*i]=0;
				if(i%p[j]==0){phi[p[j]*i]=p[j]*phi[i],mu[p[j]*i]=0;break;}
				else phi[p[j]*i]=phi[p[j]]*phi[i],mu[p[j]*i]=-mu[i];
			}
		}
	}
};