质数筛法
朴素筛法
枚举素数即可,判断一个数 x x x是不是素数,从 2 2 2开始到 x \sqrt{x} x枚举他的约数即可。
因为一个合数 x x x可以写成是两个约数的乘积 m n = x mn=x mn=x,假设 m > x m \gt \sqrt{x} m>x,那么 n < x n \lt \sqrt{x} n<x,因此,对于任意的合数 x x x,必定存在一个约数 k ≤ x k \leq \sqrt{x} k≤x。
class Solution {
public:
bool isPrime(int x) {
for (int i = 2; i * i <= x; ++i) {
if (x % i == 0) {
return false;
}
}
return true;
}
int countPrimes(int n) {
int ans = 0;
for (int i = 2; i < n; ++i) {
ans += isPrime(i);
}
return ans;
}
};
时间复杂度为 O ( n n ) O(n\sqrt{n}) O(nn),空间复杂度为 O ( 1 ) O(1) O(1)。
埃氏筛
我们考虑这样一个事实:如果 x x x是质数,那么大于 x x x的 x x x的倍数 2 x , 3 x , … 2x,3x,\ldots 2x,3x,… 一定不是质数,因此我们可以从这里入手。又因为唯一素约数分解定理,一个合数必定可以写成一个小于他的素数和一个整数的乘积,因此我们只需要碰见一个素数,就关掉它的倍数即可。
class Solution {
public:
int countPrimes(int n) {
vector<int> isPrime(n, 1);
int ans = 0;
for (int i = 2; i < n; ++i) {
if (isPrime[i]) {
ans += 1;
if ((long long)i * i < n) {
for (int j = i * i; j < n; j += i) {
isPrime[j] = 0;
}
}
}
}
return ans;
}
};
假设 x x x是质数,应该从 x ∗ x x * x x∗x开始关掉,因为 2 x , 3 x , 4 x , 5 x , … , x ( x − 1 ) 2x,3x,4x,5x,\ldots,x(x-1) 2x,3x,4x,5x,…,x(x−1),已经被 2 , 3 , 4 , 5 , … , x − 1 2,3,4,5,\ldots,x-1 2,3,4,5,…,x−1关掉,不需要再考虑了。
开根优化
同时可以加上开根优化,对于枚举小于等于 n n n的素数,我们只需要筛到 n \sqrt{n} n即可,因为,如果有大于 n \sqrt{n} n的合数 m m m,他必定有小于 m \sqrt{m} m的素因子,因为 m ≤ n \sqrt{m} \leq \sqrt{n} m≤n,因此合数 m m m一定能够被关掉。
int primes[100000];
void eratosthenes(int n)
{
int i = 0;
for (; i <= n; i++)
primes[i] = 1;
primes[0] = primes[1] = 0;
for (int i = 2; i * i <= n; i++)
{
if (primes[i])
{
int j = i;
for (; i * j <= n; j++)
primes[i * j] = 0;
}
}
}
奇数优化
除了偶素数2之外的偶数都是合数,因此跳过偶数只处理奇数。
因为素数一定是奇数,我们优化一下埃氏筛:
- 只在奇数范围内寻找素数
- 奇数乘以偶数还是偶数,因此没有必要标记,因为在第一步就已经排除了,所以我们只标记奇合数就好了。
var countPrimes = function(n) {
var isCom= new Int8Array(n), b = Math.sqrt(n), r = n > 2 ? 1 : 0, t
for(var i = 3; i < n; i += 2) {
if(!isCom[i]) {
r++
if (i <= b) for(var j = i; (t = i * j) < n; j += 2) isCom[t] = 1
}
}
return r
};
可以证明,埃氏筛法的时间复杂度为 O ( n log log n ) O(n \log \log n) O(nloglogn)。
线性筛
在上面的埃氏筛中,一个合数可能被关闭多次,我们就有了优化的方向,让一个合数只被关闭一次即可。
根据唯一素因子分解定理,一个合数等于一个素约数乘以另外一个约数,那么这么说,我们可以建立一个线性表,里面储存了我们到目前为止遇到的素数,然后迭代途中,遇到一个数,就依次关掉这个线性表中的元素与这个数的乘积即可。
可以证明,线性表primes中的元素一定是线性递增的。我们可以证明这个算法的正确性:
- 如果遇到的是素数,那么只能被他自己和1分解,而在这之前没有理由关掉这个素数,因此这个数一定是素数。
- 如果遇到一个合数,那么这个合数一定可以写作是,一个小于这个数的质数 p p p,和大于这个质数 p p p但小于这个数的一个数 k k k的乘积 p k pk pk。(有时候也可以是小于这个质数 p p p的 k k k,如果 k k k是质数,那么 p p p和 k k k互相交换一下就好,如果 k k k是合数,那么一定能在 k k k中找到一个素数,把它当做 p p p就好,总之,这两个数 p p p和 k k k一定能找到) p p p一定是在 k k k访问之前就已经被添加到素数列表中了,之后再访问 k k k的时候,就会关掉合数 p k pk pk了。
class Solution {
public:
int countPrimes(int n) {
vector<int> primes;
vector<int> isPrime(n, 1);
for (int i = 2; i < n; ++i) {
if (isPrime[i]) {
primes.push_back(i);
}
for (int j = 0; j < primes.size() && i * primes[j] < n; ++j) {
isPrime[i * primes[j]] = 0;
}
}
return primes.size();
}
};
但是反过来想一想,这么做一个合数还是可能被访问多次啊,比如 40 = 20 × 2 = 5 × 8 40 = 20 \times 2 = 5 \times 8 40=20×2=5×8,当访问到20的时候,40会被关掉一次,当访问到8的时候,40又被关掉了一次,根本没有做到优化,能不能让一个合数被唯一一个素数确定呢?
答案是肯定的,就是最小素因子,即一个合数,一定存在一个最小素因子 p p p,和一个数 k k k的乘积 p k pk pk,并且在 k k k的分解中,没有比 p p p还小的素因子了。此时,一个合数就能被唯一一个素数所确定。
那我们怎么改进代码呢?添加一句话:
class Solution {
public:
int countPrimes(int n) {
vector<int> primes;
vector<int> isPrime(n, 1);
for (int i = 2; i < n; ++i) {
if (isPrime[i]) {
primes.push_back(i);
}
for (int j = 0; j < primes.size() && i * primes[j] < n; ++j) {
isPrime[i * primes[j]] = 0;
if (i % primes[j] == 0) {
break;
}
}
}
return primes.size();
}
};
其中
if (i % primes[j] == 0) {
break;
}
也就是说,如果遇到一个素数可以整除i了,那么就停止标记。这是为什么呢?假如继续标记 i ∗ p r i m e s [ j + 1 ] i * primes[j+1] i∗primes[j+1]这个数,又因为 i i i可以被 p r i m e s [ j ] primes[j] primes[j]整除,也就是说 i ∗ p r i m e s [ j + 1 ] = p r i m e s [ j ] ∗ n ∗ p r i m e s [ j + 1 ] i * primes[j+1] = primes[j] * n * primes[j+1] i∗primes[j+1]=primes[j]∗n∗primes[j+1],其中 n = i / p r i m e s [ j ] n = i / primes[j] n=i/primes[j],并且 n n n不能再被分解为比 p r i m e s [ j ] primes[j] primes[j]更小的素数了(如果能分解,那么 i i i还能被比 p r i m e s [ j ] primes[j] primes[j]更小的质数整除了)因此,数 p r i m e s [ j ] ∗ n ∗ p r i m e s [ j + 1 ] primes[j] * n * primes[j+1] primes[j]∗n∗primes[j+1]分解中 p r i m e s [ j ] primes[j] primes[j]就是最小的素因子了,而 p r i m e s [ j + 1 ] primes[j+1] primes[j+1]不是,所以不能再这里标记,应该在数 n ∗ p r i m e s [ j + 1 ] n * primes[j + 1] n∗primes[j+1]的时候标记(这个数大于 i i i,肯定能在之后访问到)。其他数也同理。
总结一句话就是,针对被关闭的数 k = i ∗ p r i m e s [ j ] k = i * primes[j] k=i∗primes[j],必须令 p r i m e s [ j ] primes[j] primes[j]是其最小的素因子。
时间复杂度 O ( n ) O(n) O(n)。
筛法求数论函数
利用筛法,我们可以做很多事情,其中最有用的就是求欧拉函数和莫比乌斯函数。
求欧拉函数
首先我们要明白一点,欧拉函数是积性函数。
先看代码:
void init() {
phi[1] = 1;
for (int i = 2; i < MAXN; ++i) {
if (!vis[i]) {
phi[i] = i - 1;
pri[cnt++] = i;
}
for (int j = 0; j < cnt; ++j) {
if (1ll * i * pri[j] >= MAXN) break;
vis[i * pri[j]] = 1;
if (i % pri[j]) {
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
} else {
// i % pri[j] == 0
// 换言之,i 之前被 pri[j] 筛过了
// 由于 pri 里面质数是从小到大的,所以 i 乘上其他的质数的结果一定也是
// pri[j] 的倍数 它们都被筛过了,就不需要再筛了,所以这里直接 break
// 掉就好了
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
}
}
}
对于: i i i是素数的情况,很显然 φ ( i ) = i − 1 \varphi(i)=i-1 φ(i)=i−1。
对于:phi[i * pri[j]] = phi[i] * (pri[j] - 1);
来说,
i
i
i与
p
r
i
[
j
]
pri[j]
pri[j]互质,因此利用积性函数的性质求解即可。
对于:phi[i * pri[j]] = phi[i] * pri[j];
来说,
p
r
i
[
j
]
pri[j]
pri[j]是
i
i
i的一个因子,我们把数进行质因数分解
i
∗
p
r
i
[
j
]
=
p
r
i
[
j
]
m
k
i * pri[j] = pri[j]^{m}k
i∗pri[j]=pri[j]mk,因为两者互质,因此
φ ( p r i [ j ] m k ) = φ ( p r i [ j ] m ) φ ( k ) = ( p r i [ j ] m − p r i [ j ] m − 1 ) φ ( k ) = p r i [ j ] ( p r i [ j ] m − 1 − p r i [ j ] m − 2 ) φ ( k ) = p r i [ j ] φ ( i ) \varphi(pri[j]^{m}k) = \varphi(pri[j]^{m}) \varphi(k) = (pri[j]^{m} - pri[j]^{m-1}) \varphi(k) = pri[j] (pri[j]^{m-1} - pri[j]^{m-2}) \varphi(k) = pri[j] \varphi(i) φ(pri[j]mk)=φ(pri[j]m)φ(k)=(pri[j]m−pri[j]m−1)φ(k)=pri[j](pri[j]m−1−pri[j]m−2)φ(k)=pri[j]φ(i)
因此该算法是正确的。
求莫比乌斯函数
根据莫比乌斯函数的定义:
当 n = 1 n=1 n=1的时候 μ ( n ) = 1 \mu(n)=1 μ(n)=1;当 n = p 1 p 2 … p m n=p_{1}p_{2} \ldots p_{m} n=p1p2…pm的时候,即任意多个质数相乘并且没有幂次质数, μ ( n ) = ( − 1 ) m \mu(n)=(-1)^{m} μ(n)=(−1)m;其他情况 μ ( n ) = 0 \mu(n)=0 μ(n)=0
#define MAXN 1000
int phi[MAXN];
bool vis[MAXN];
int pri[MAXN];
int mu[MAXN];
int cnt = 0;
void euler()
{
phi[1] = 1;
mu[1] = 1;
for (int i = 2; i < MAXN; ++i)
{
if (!vis[i])
{
phi[i] = i - 1;
// 唯一一个质数为-1
mu[i] = -1;
pri[cnt++] = i;
}
for (int j = 0; j < cnt; ++j)
{
if (1ll * i * pri[j] >= MAXN)
break;
vis[i * pri[j]] = 1;
if (i % pri[j])
{
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
// 利用莫比乌斯函数的积性
mu[i * pri[j]] = -mu[i];
}
else
{
// i % pri[j] == 0
// 换言之,i 之前被 pri[j] 筛过了
// 由于 pri 里面质数是从小到大的,所以 i 乘上其他的质数的结果一定也是
// pri[j] 的倍数 它们都被筛过了,就不需要再筛了,所以这里直接 break
// 掉就好了
phi[i * pri[j]] = phi[i] * pri[j];
// pri[j] 是 i的倍数,存在幂次质数,函数值为0
mu[i * pri[j]] = 0;
break;
}
}
}
}
求约数个数
除法函数:
σ k ( n ) = ∑ d ∣ n d k \sigma_{k}(n)=\sum_{d \mid n}d^{k} σk(n)=d∣n∑dk
为积性函数。
则当 k = 0 k=0 k=0的时候,除法函数退化为因数个数函数。仍为积性函数。
#define MAXN 1000
bool vis[MAXN];
int pri[MAXN];
int pnxt[MAXN]; // 表示去掉所有质因子之后,剩下的质因子的乘积
int sigma[MAXN];
int cnt = 0;
void euler()
{
pnxt[1] = 1;
sigma[1] = 1;
for (int i = 2; i < MAXN; ++i)
{
if (!vis[i])
{
// 质数的约数只有自己和1
sigma[i] = 2;
// 质数的最小质因子的数量为1,去掉后乘积为1
pnxt[i] = 1;
pri[cnt++] = i;
}
for (int j = 0; j < cnt; ++j)
{
if (1ll * i * pri[j] >= MAXN)
break;
vis[i * pri[j]] = 1;
if (i % pri[j])
{
// 积性
sigma[i * pri[j]] = sigma[i] * 2;
// i * pri[j] 只有一个最小质数因子为 pri[j]
pnxt[i * pri[j]] = i;
}
else
{
// 添加了一个最小质因子
pnxt[i * pri[j]] = pnxt[i];
// pri[j] 是 i 的一个素约数
sigma[i * pri[j]] = sigma[i] + sigma[pnxt[i]];
break;
}
}
}
}
对于sigma[i * pri[j]] = sigma[i] + sigma[pnxt[i]];
来说:
σ ( i ∗ p ) = σ ( p m k ) = σ ( p m ) σ ( k ) = ( m + 1 ) σ ( k ) = ( σ ( p m − 1 ) + 1 ) ) σ ( k ) = σ ( i ) + σ ( k ) \sigma(i * p) = \sigma(p^{m}k) = \sigma(p^{m}) \sigma(k) = (m + 1) \sigma(k) = (\sigma(p^{m-1}) + 1)) \sigma(k) = \sigma(i) + \sigma(k) σ(i∗p)=σ(pmk)=σ(pm)σ(k)=(m+1)σ(k)=(σ(pm−1)+1))σ(k)=σ(i)+σ(k)
其中 k k k 为去掉所有最小质因子后剩下的乘积。
求约数和
#define MAXN 1000
bool vis[MAXN];
int pri[MAXN];
int pnxt[MAXN]; // 表示去掉所有质因子之后,剩下的质因子的乘积
int psum[MAXN]; // 表示最下质因子的幂次和
int sigma1[MAXN];
int cnt = 0;
void euler()
{
pnxt[1] = 1;
sigma1[1] = 1;
for (int i = 2; i < MAXN; ++i)
{
if (!vis[i])
{
// 质数的约数只有自己和1
sigma1[i] = 1 + i;
// 质数的最小质因子的数量为1,去掉后乘积为1
pnxt[i] = 1;
// p0 + p1
psum[i] = 1 + i;
pri[cnt++] = i;
}
for (int j = 0; j < cnt; ++j)
{
if (1ll * i * pri[j] >= MAXN)
break;
vis[i * pri[j]] = 1;
if (i % pri[j])
{
// 积性
sigma1[i * pri[j]] = sigma1[i] * (1 + pri[j]);
// i * pri[j] 只有一个最小质数因子为 pri[j]
pnxt[i * pri[j]] = i;
psum[i * pri[j]] = pri[j] + 1;
}
else
{
// 添加了一个最小质因子
pnxt[i * pri[j]] = pnxt[i];
// 将所有项都递推一步 并 加上 p0
psum[i * pri[j]] = psum[i] * pri[j] + 1;
// pri[j] 是 i 的一个素约数
sigma1[i * pri[j]] = psum[i * pri[j]] * sigma1[pnxt[i]];
break;
}
}
}
}
对于sigma1[i * pri[j]] = psum[i * pri[j]] * sigma1[pnxt[i]];
来说:
σ 1 ( i ∗ p ) = σ 1 ( p m ) σ 1 ( k ) = σ 1 ( k ) ∑ i = 0 m p i \sigma_{1}(i * p) = \sigma_{1}(p^{m}) \sigma_{1}(k) = \sigma_{1}(k) \sum_{i = 0}^{m}p^{i} σ1(i∗p)=σ1(pm)σ1(k)=σ1(k)i=0∑mpi
字母的含义同上。
求最大质因子
void euler()
{
for (int i = 2; i < 5000005; i++)
{
if (!isNprime[i])
{
smf[i] = i;
// 素数的最大质因子是其本身
primes[tot++] = i;
}
for (int j = 0; j < tot && i * primes[j] < 5000005; j++)
{
int val = i * primes[j];
isNprime[val] = true;
// 我们知道,线性筛的质数是从小到大添加的,因此最大质因子一定是原来数i的最大质因子
smf[val] = smf[i];
if (i % primes[j] == 0)
{
break;
}
}
}
}
求最小质因子
void euler()
{
for (int i = 2; i < 5000005; i++)
{
if (!isNprime[i])
{
smf[i] = i;
// 素数的最小质因子是其本身
primes[tot++] = i;
}
for (int j = 0; j < tot && i * primes[j] < 5000005; j++)
{
int val = i * primes[j];
isNprime[val] = true;
// 我们知道,线性筛的质数是从小到大添加的,因此最小质因子一定是新添加的素数primes[j]
smf[val] = primes[j];
if (i % primes[j] == 0)
{
break;
}
}
}
}