数学计数——基础筛法

[参考资料：筛法（OI-wiki）](https://oi-wiki.org/math/number-theory/sieve/) ### Chapter 1. 知识点 注：以下整数均指 $>1$ 的整数。 首先声明之后代码里的主要变量含义：$\texttt{cnt}$ 表示素数个数；$\texttt{p}$ 为记录所有素数的数组；$\texttt{vis}$ 表示每个数是否是素数。 1. 普通筛法 - 思想：整数的整数倍数是合数； - 时间复杂度：$O(n\log n)$。 ~~~cpp for(int i=2;i<=n;i++) { if(!vis[i]) p[++cnt]=i; for(int j=2*i;j<=n;j+=i) //整数的整数倍数是合数 vis[j]=1; } ~~~ 2. 埃氏筛 - 素数的整数倍数是合数； - 时间复杂度：$O(n\log \log n)$。 ~~~cpp for(int i=2;i<=n;i++) if(!vis[i]) { p[++cnt]=i; for(int j=2*i;j<=n;j+=i) //素数的整数倍数是合数 vis[j]=1; } ~~~ 3. 线性筛 - 思想：整数的素数倍数是合数； - 时间复杂度：$O(n)$。 - 时间复杂度的证明：枚举整数的素数倍数，当枚举到的素数是当前整数的约数时停止。设 $n$ 的最小质因子是 $x_n$，$n$ 会且仅会在 $i=\frac{n}{x_n}$ 的 $x_n$ 倍时倍筛到。所以每个数仅仅会被筛到一次，时间复杂度是线性的。 ~~~cpp for(int i=2;i<=n;i++) { if(!vis[i]) p[++cnt]=i; for(int j=1;j<=cnt&&p[j]*i<=n;j++) //整数的素数倍数是合数 { vis[p[j]*i]=1; if(i%p[j]==0) break; } } ~~~ 4. 区间筛 - 此算法为埃氏筛思想的一种延申运用。 - 问题描述：求区间 $[a,b]$ 内所有的素数。$a\le b\le 10^{12},b-a\le 10^6$。 - 算法流程：1. 线性筛出 $1\sim \sqrt b$ 内的所有素数；2. 运用埃氏筛的思想筛掉 $[a,b]$ 中的所有合数。 ~~~cpp init(ceil(sqrt(r))); //线性筛 //筛去 [a,b]中的合数（埃氏筛） for(int i=1;i<=cnt;i++) //枚举每一个在 [1,sqrt(b)]内的素数 for(int j=(a-1)/pri[i]+1;1ll*j*pri[i]<=b;j++) ss[1ll*j*pri[i]-a]=1; //ss[i]=1/0: a+i是否为合数 int ans=0; for(int i=0;i<=r-l;i++) //注意循环上下界 if(!ss[i]) ans++; ~~~ 5. 线性筛求积性函数* - **后续会继续补充该方面的知识。** - 积性函数的定义：考虑一个定义在正整数上的函数 $f$，若 $\forall a,b\in \mathbb{N}^+,\gcd(a,b)=1$，都有 $f(a)\times f(b)=f(ab)$，则称 $f$ 为积性函数。 - 常见的两个积性函数： - 欧拉函数 $\varphi$：$\varphi(n)=\sum_{k=1}^n [\gcd(i,n)=1]$。 - 莫比乌斯函数 $\mu$：当 $n$ 有平方因子时，$\mu(n)=0$；否则 $\mu(n)=(-1)^k$，$k$ 为 $n$ 的素因子个数。 - 可以用线性筛计算欧拉函数的原因；对于每个合数 $n$，在线性筛的过程中都会拆成 $\frac{n}{x_n}$ 和 $x_n$。设 $\gcd(x_n^∞,n)=x_n^k$，则 $\gcd(\frac{n}{x_n^k},x_n^k)=1$。 - 算法流程：1. 直接计算出 $f(x_n^k)$；2. 对于其他合数 $n$ ，$f(n)=f(\frac{n}{x_n^k})\times f(x_n^k)$。 - 下面给出以计算欧拉函数为例的核心代码（$\texttt{phi}$ 表示欧拉函数值）： ~~~cpp for(int i=2;i<=n;i++) { if(!vis[i]) { p[++cnt]=i; phi[i]=i-1; //由于欧拉函数的特性，只需特殊计算 phi(p) 即可 } for(int j=1;j<=t&&p[j]*i<=n;j++) //整数的素数倍数是合数 { vis[p[j]*i]=1; if(i%p[j]) phi[p[j]*i]=phi[i]*(p[j]-1); else { phi[p[j]*i]=phi[i]*p[j]; break;} } } ~~~ ### Chapter 2. 题目 1. [P3383 【模板】线性筛素数](https://www.luogu.com.cn/problem/P3383) 2. [P3912 素数个数](https://www.luogu.com.cn/problem/P3912) - 两道线性筛模板题。 3. [P1835 素数密度](https://www.luogu.com.cn/problem/P1835) - 区间筛模板题。 4. [P1865 A % B Problem](https://www.luogu.com.cn/problem/P1865) - 线性筛模板 $+$ 前缀和。 第 $1\sim 4$ 题都是比较简单的模板题，之后的题将会有些难度，很多是涉及到欧拉函数的题。 5. [P2568 GCD](https://www.luogu.com.cn/problem/P2568) - 若 $\gcd(x,y)=p$（ $p$ 为素数，$x\le y$），则 $\gcd(\frac{x}{p},\frac{y}{p})=1$。 - 我们对每个素数 $p$ 考虑满足 $1\le x'\le y'\le \frac{n}{p},\gcd(x',y')=1$ 的正整数对 $(x',y')$ 的个数。 - 那么先考虑对 $1\le y\le p$ 有多少个 $1\le x\le y$ 满足 $\gcd(x,y)=1$。可以发现这就是欧拉函数的定义。 - 那么我们只需要通过线性筛预处理欧拉函数，然后求出欧拉函数的前缀和，最后再对每个素数算出结果并相加得到答案即可。 ~~~cpp init(); //线性筛 + 预处理欧拉函数 for(int i=1;i<=n;i++) sum[i]=sum[i-1]+phi[i]; //欧拉函数前缀和 for(int i=1;i<=cnt&&pri[i]<=n;i++) ans+=2*sum[n/pri[i]]-1; //求和 ~~~ 6. [P2303 [SDOI2012] Longge 的问题](https://www.luogu.com.cn/problem/P2303) - $\sum\limits_{i=1}^n \gcd(i, n)=\sum\limits_{k=1}^n[n\bmod k=0]k\times\varphi(\frac{n}{k})$。 - 时间复杂度：$O(\text{因子个数}\times\sqrt n)$，直接计算可过。 ~~~cpp long long phi(long long x) { long long res=x; for(long long i=2;i*i<=x;i++) if(x%i==0) //i为 x的素因数 { res=res/i*(i-1); while(x%i==0) x/=i; } if(x>1) res=res/x*(x-1); return res; } long long getans(long long n) { long long ans=0,i=1; while(i*i<n) { if(n%i==0) ans+=(i+(n/i))*getphi(i); //同时计算 i和 x/i的贡献 i++; } if(i*i==x) ans+=i*getphi(i); //sqrt(x)和 x/sqrt(x) return ans; } ~~~ 本题还有更好的做法：[link](https://www.luogu.com.cn/blog/PinkRabbit/solution-p2303)。 7. [P2158 [SDOI2008] 仪仗队](https://www.luogu.com.cn/problem/P2158) - 设左下角坐标为 $(0,0)$，则点能被看到的点的坐标 $(x,y)$ 满足 $\gcd(x,y)=1$。 - 所以答案就是 $2(\varphi(1)+\cdots+\varphi(n))+1$。 8. [P2398 GCD SUM](https://www.luogu.com.cn/problem/P2398) - 考虑对每个 $d=1,2,\cdots,n$ 计算满足 $\gcd(i,j)=d$ 的 $(i,j)$ 有多少个，即 $\gcd(x,y)=1,x,y\le \frac{n}{d}$，由第七题就可以知道答案。 9. [P1390 公约数的和](https://www.luogu.com.cn/problem/P1390) - 与第八题一样的。 10. [P1447 [NOI2010] 能量采集](https://www.luogu.com.cn/problem/P1447) - 就是要求 $\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)$ 的值。 ```cpp const int MAXN=1e5+10; int n,m; long long num[MAXN],ans; int main() { scanf("%d%d",&n,&m); if(n>m) swap(n,m); for(int i=n;i;i--) //注意循环的顺序 { num[i]=1ll*(n/i)*(m/i); for(int j=2*i;j<=n;j+=i) num[i]-=num[j]; ans+=(2*i-1)*num[i]; //计算贡献 } printf("%lld",ans); return 0; } ```