【数学】素数筛法、判断法与数的分解

更新于 分类于 ， 阅读次数： 本文字数： 4.8k 阅读时长 ≈ 4 分钟

素数筛法（埃式筛、欧拉筛、区间筛）、素数判断法（朴素、模6、Miller-Rabin及改进）、数的分解（朴素、Pollard-rho）

素数筛法

埃式筛法

O(nloglogn)O(nloglogn)

1
2
3
4
5
6
7
8
9
10
11
12
bool b[N];
bool prime[N];
void solve(int x)
{
    for (int i = 0; i <= x; i++)
        prime[i] = true;
    b[0] = b[1] = 0;
    for (int i = 2; i <= x; i++)
        if (b[i])
            for (int j = i << 1; j <= x; j += i)
                prime[j] = false;
}

欧拉筛法

O(n)O(n)

过程：以check[i]=0记i为素数。筛ii 时，已知质数prime[0...tot]，则i×prime[j]i\times prime[j] 一定不是质数，如果超出了要求的范围，或prime[j]∣iprime[j]|i 时jj 不再增加。

停止条件的原因：对于prime[j]∣i,i=prime[j]×kprime[j]|i,i=prime[j]\times k ，未来一定会筛到x=prime[j+1]×i=prime[j]×k×prime[j+1]x=prime[j+1]\times i=prime[j] \times k\times prime[j+1] 。

未来ii 一定会增到i′=k×prime[j+1]>k×prime[j]=ii'=k\times prime[j+1]>k\times prime[j]=i ，这一轮一定会筛到i×prime[j+1]i\times prime[j+1] 筛去。

跳出循环时prime[j]prime[j] 正是prime[j]×iprime[j]\times i 的最小质因子。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
const int MAX_L = (1e7) + 2;
int a, b;
bool check[MAX_L]; //检查某个数是否是素数，初始化置false
int prime[MAX_L];
int tot = 0;
void sieve(int x)
{
    for (int i = 2; i <= x; i++)
    {
        if (!check[i])
            prime[tot++] = i; //i是素数，添加到prime中
        for (int j = 0; j < tot; j++)
        {
            if (prime[j] * i > x)
                break;
            check[prime[j] * i] = true;
            if (i % prime[j] == 0)
                break; //已被筛过，可以不往下筛了
        }
    }
}

区间筛素数

[a,b)[a,b) ，小数区间推大数区间

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
typedef long long ll;
const int MAX_L = (1e8) + 2;
const int MAX_SQRTB = (1e4) + 2;
ll a, b;
bool isprime[MAX_L];           //[a,b)区间，用isprime[i-a]=true表示i为素数
bool isprime_small[MAX_SQRTB]; //[2,sqrt(b))区间
void segment_sieve(ll x, ll y)
{
    for (ll i = 0; (ll)i * i < y; i++)
        isprime_small[i] = true;
    for (ll i = 0; i < y - x; i++)
        isprime[i] = true;
    for (int i = 2; (ll)i * i < y; i++)
        if (isprime_small[i])
        {
            for (ll j = 2 * i; (ll)j * j < y; j += i)
                isprime_small[j] = false; //从2开始筛[2,sqrt(b))区间
            for (ll j = max((ll)2, (x + i - 1) / i) * i; j < y; j += i)
                //max找离x最近的i的倍数，如果没有就是2倍i
                isprime[j - a] = false;
        }
}

素数判断法

朴素判断法

  • 复杂度O(n)O(\sqrt n)。

  • 优化：22 之后，只有奇数有可能是素数，所以先判≤2\leq 2 的情况，然后从i=3i=3 每次循环i+=2i+=2

1
2
3
for (int i = 2; i <= sqrt(n); i++)
    if (x % i == 0)
        return false;

模6判断法

1
2
3
4
5
6
7
8
9
10
11
12
bool num[5] = {0, 0, 1, 1, 0};
bool prime(int x)
{
    if (x <= 4)
        return num[x]; //4以下直接判断
    if (x % 6 != 1 && x % 6 != 5)
        return false; //6k+2, 6k+3, 6k+4分别有因数2，3，2
    for (int i = 5; i * i <= x; i += 6)
        if (x % i == 0 || x % (i + 2) == 0)
            return false; //i从5开始，以6递增，判断是不是6k+5或6k+1的因数
    return true;          //都不是，就是素数
} //以及对于大数小数k*k<x和k<sqrt(x)差不多快

Miller-Rabin判断法

由费马小定理，n是素数,p∤a⇒an−1≡1(mod  n)n是素数,p\nmid a \Rightarrow a^{n-1}\equiv 1(\mod n) ，反过来不一定。但是概率小，重复操作约30次差不多能保证判断正确。复杂度O(log⁡n)O(\log n) 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// qpow(x,n,m): 计算x^n mod m
// 记得srand
bool Miller_Rabin(int x)
{
    if (x < 2)
        return false;
    if (x == 2)
        return true;
    if ((x & 1) == 0)
        return false;
    for (int i = 1; i <= 30; i++)
    {
        int base = rand() % (x - 1) + 1; // 随机[1,x)的数
        if (qpow(base, x - 1, x) != 1)
            return false;
    }
    return true;
}

改进Miller-Rabin判断法

二次探测定理：如果 pp 是一个素数, 且 0<x<p0<x<p, 则方程 x2%p=1x^{2} \% p=1 的解为 x=1x=1 或 x=p−1x=p-1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
// 通过 a^(n−1)=1(mod n)来判断 n 是不是素数，并进行二次判断
// 合数返回true
bool check(ll a, ll n, ll x, ll t)
{
    ll ret = qpow(a, x, n); // 快速幂
    ll last = ret;
    for (int i = 1; i <= t; i++)
    {
        ret = multi_add(ret, ret, n); // 乘法改加法防止溢出
        if (ret == 1 && last != 1 && last != n - 1)
            return true; //合数
        last = ret;
    }
    if (ret != 1)
        return true;
    return false;
}
bool Miller_Rabin(ll n)
{
    if (n < 2)
        return false;
    if (n == 2)
        return true;
    if ((n & 1) == 0)
        return false;
    ll x = n - 1, t = 0;
    while ((x & 1) == 0)
        x >>= 1, t++;
    for (int i = 0; i < 8; i++) // 8次差不多了
        if (check(rand() % (n - 1) + 1, n, x, t))
            return false;
    return true;
}

数的分解

普通整数分解

  • 枚举因子试除
  • 素数打表试除（kuangbin的模板）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// prime[i]中存的是第i个素数
ll factor[100][2]; // factor[i][0]表示给定的数第i个素因子，[1]表示系数
int fatCnt;
int getFactors(ll x) // 返回素因子个数
{
    fatCnt = 0;
    ll tmp = x;
    for (int i = 1; prime[i] <= tmp / prime[i]; i++)
    {
        factor[fatCnt][1] = 0;
        if (tmp % prime[i] == 0)
        {
            factor[fatCnt][0] = prime[i];
            while (tmp % prime[i] == 0)
            {
                factor[fatCnt][1]++;
                tmp /= prime[i];
            }
            fatCnt++;
        }
    }
    if (tmp != 1)
    {
        factor[fatCnt][0] = tmp;
        factor[fatCnt++][1] = 1;
    }
    return fatCnt;
}

大整数分解：Pollard-rho

本质：随机因子检查

Pollard的伪随机数生成器：f(x)=(x2+c)mod  Nf(x)=(x^2+c)\mod N ，随机数序列：{x,f(x),f(f(x))⋯ }\{x,f(x),f(f(x))\cdots\} 。可以修改参数打表画图得到，这个序列会进入循环，常类似ρ\rho 型，即经过一些数才进入循环

（这是kuangbin的模板）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
ll factor[100]; //质因素分解结果（刚返回时时无序的）
int tol;        //质因素的个数，编号 0∼tol-1
//找出一个因子
ll pollard_rho(ll x, ll c)
{
    ll i = 1, k = 2;
    srand(time(NULL));
    ll x0 = rand() % (x - 1) + 1;
    ll y = x0;
    while (1)
    {
        i++;
        x0 = (multi_add(x0, x0, x) + c) % x;
        ll d = __gcd(y - x0, x);
        if (d != 1 && d != x)
            return d;
        if (y == x0)
            return x;
        if (i == k)
        {
            y = x0;
            k += k;
        }
    }
}
//对 n 进行素因子分解，存入 factor; k 设置为 107 左右即可
void findfac(ll n, int k)
{
    if (n == 1)
        return;
    if (Miller_Rabin(n))
    {
        factor[tol++] = n;
        return;
    }
    ll p = n;
    int c = k;
    while (p >= n)
        p = pollard_rho(p, c--); //值变化，防止死循环 k findfac(p,k);
    findfac(p, k);
    findfac(n / p, k);
}