算法心传·素数判定及素数筛

一、素数(质数)

1.引入

素数又称质数,是大于1的自然数中只能被1和它本身整除的数。相对的,如果一个大于1的自然数能被除了1和它本身以外的整数整除的数叫做合数。特别的,1既不是质数也不是合数。

素数有着包括以下几条有趣的事实:

  • 素数的数量有无限多。
  • 素数的分布随着整数的增大分布越来越稀疏。
  • 对于任意正整数n,至少存在n个连续的正合数。

同时,有大量关于素数的猜想,以下是几个著名的猜想:

  • 波特兰猜想:对于任意一个给定的正整数n>1,存在一个素数p,使得n<p<2n。已经证明。
  • 孪生素数猜想:存在无穷多个形如p和p+2的素数对。尚未证明。
  • 素数等差数列猜想:对于任意正整数n>2,有一个由素数组成的长度为n的等差数列。已证明存在任意长的素数等差数列。
  • 哥德巴赫猜想:每个大于2的正偶数可以写成两个素数的和。最有名的素数猜想,尚未证明。

二、单素数判定

1.素数的判定

素数的判定方法主要有三种,第一种是试除法,我们可以根据素数的定义从2逐个罗列,试着去整除n,如果n可以被其中任何一个数整除就说明其不是素数。很容易发现,我们可以把范围从[2,n-1]缩小到[2,sqrt(n)]。下面给出解决问题的函数:

// ll 代替 long long 详见竞赛模板
bool prime(ll a){
    if(a < 2)
        return false;
    for (int i = 2; i <= sqrt(a);++i)
        if(a % i == 0)
            return false;
    return true;
}

代码很好理解,这里不过多解释。不难发现,当a小于1e12时,开方的复杂度是可以运行不超时的,但是随着a的进一步变大,素数的判定用试除法就有点不够用了,这时候就要引出下面的两种方法了。

2.Fermat 素性测试

Fermat 素性测试中文名叫费马素性测试,他基于费马小定理。

费马小定理:设n是素数,a是正整数且与n互素,那么有an-1≡1(mod n)。费马小定理的逆命题也几乎成立,Fermat 素性测试就是基于费马小定理的逆命题,下面详细介绍方法。

为了测试n是否为素数,在1~n之间任选一个随机的基值a,a并不一定要和n互素。

  1. 如果an-1≡1(mod n)不成立,那么n肯定不是素数。
  2. 如果an-1≡1(mod n)成立,那么n大概率是素数,尝试的a越多,n是素数的概率越大,我们称n是一个基于a的伪素数。

但是,可惜的是从第二种情况可以看出Fermat 素性测试并不是完全正确的。对于某个a的值,总有一些合数被误判通过测试。不同的a值,被误判的合数都不太一样,但是特别的,总有一部分合数无论a怎么选,都能通过测试,我们称其为Carminchsel数,中文名叫做卡迈克尔数。前三个Carminchsel数是561,1105和1729。Carminchsel数很少,前1亿的数里面也就只有255个,并且当n趋向于无穷时,Carminchsel数的分布极其稀疏,Fermat 素性测试几乎不会出错,所以它一种相当不错的方法,下面给出它的实现函数:

ll mul(ll a, ll b, ll mod)
{
    // return (ll)((__int128)a * b % mod);  // 环境允许选择这一条语句
    // 先将 a 和 b 对 mod 取模,缩小范围(但若 mod 很大,仍可能溢出)
    a %= mod;
    b %= mod;

    ll res = 0; // 累加结果
    while (b > 0)
    {
        if (b & 1)
        {                          // 如果当前二进制位为1,则累加当前的 a
            res = (res + a) % mod; // 注意:res + a 可能溢出(当 mod 接近 LLONG_MAX)
        }
        a = (a + a) % mod; // a 翻倍,对应二进制位的权重
        b >>= 1;           // 处理下一位
    }
    return res;
}

ll qpow(ll a, ll b, ll mod)
{
    ll res = 1;
    a %= mod;
    while (b)
    {
        if (b & 1)
            res = mul(res, a, mod); // 当前位为1时累乘
        a = mul(a, a, mod);         // 底数平方
        b >>= 1;                    // 指数右移
    }
    return res;
}

bool fermat(ll n, int k = 50)
{
    if (n < 2)
        return false;
    for (int i = 0; i < k; ++i)
    {
        ll a = rand() % (n - 1) + 1; // 随机选取随机数,范围为1到n-1
        if (qpow(a, n - 1, n) != 1)
            return false;
    }
    return true;
}

3.Miller-Rabin素性测试

Miller-Rabin素性测试中文名叫米勒-拉宾素性测试,可以看作是Fermat素性测试的改进版,可以用于判断排除Carminchsel数。Miller-Rabin素性测试用到了二次探测定理,如果p是一个奇素数,且e>=1,则方程x2≡1(mod pe)仅有两个解:x=1和x=-1。当e=1时,方程x2≡1(mod p)仅有两个解,x=1和x=p-1。我们把x=1和x=p-1称作“x对模p来说1的平凡平方根”。例如x=6,n=35,6是对模35来说的1的平凡平方根。

由二次探测定理得到一个推论:如果对模n存在1的非平凡平方根,那么n是合数。这个推论就是二次探测定理的逆否命题,即如果对n存在1的非平凡平方根,则n不可能是奇素数或者奇素数的幂。

这里对之前的Fermat素数测试的函数做出调整,如果an-1≡1(mod n)不成立,那么n肯定不是素数,由于n-1实在太大,显然直接计算an-1不可行,我们可以用到数学技巧:把n-1表示为幂的形式,然后借助快速幂来运算。

令n-1 = 2tu,其中u是奇数,t是正整数。如何求t和u呢?这里直接给出方案:n-1的二进制表示是奇数u的二进制表示后面加t个0。如255-1 = 25*7,n-1=255-1=244的二进制是11100000,它是u=7二进制111,后面加t个0。

故而公式转化为:

an-1≡(au)2t^{2^t}(mod n)

当Miller-Rabin素性测试次数达到50次时,出错概率已经小到可以忽视了,读者可以自行推导Carminchsel数561是否能被Miller-Rabin素性测试排除掉。

对于想要详细了解原理的读者可以自行去搜索学习证明过程,这里直接给出代码:

ll mul(ll a, ll b, ll mod)
{
    // return (ll)((__int128)a * b % mod);  // 环境允许选择这一条语句
    // 先将 a 和 b 对 mod 取模,缩小范围(但若 mod 很大,仍可能溢出)
    a %= mod;
    b %= mod;

    ll res = 0; // 累加结果
    while (b > 0)
    {
        if (b & 1)
        {                          // 如果当前二进制位为1,则累加当前的 a
            res = (res + a) % mod; // 注意:res + a 可能溢出(当 mod 接近 LLONG_MAX)
        }
        a = (a + a) % mod; // a 翻倍,对应二进制位的权重
        b >>= 1;           // 处理下一位
    }
    return res;
}

ll qpow(ll a, ll b, ll mod)
{
    ll res = 1;
    a %= mod;
    while (b)
    {
        if (b & 1)
            res = mul(res, a, mod); // 当前位为1时累乘
        a = mul(a, a, mod);         // 底数平方
        b >>= 1;                    // 指数右移
    }
    return res;
}

bool witness(ll a, ll n)
{
    // 将 n-1 分解为 u * 2^t,其中 u 为奇数
    ll u = n - 1;
    int t = 0;
    while ((u & 1) == 0) // 当 u 为偶数时
    {
        u = u >> 1; // 不断除以 2
        t++;        // 记录除以 2 的次数
    }

    // 计算 a^u mod n,作为初始值
    ll x1 = qpow(a, u, n); // x1 = a^u % n
    ll x2;                 // 用于存储平方后的值

    // 进行 t 次平方,依次检查是否出现非平凡平方根
    for (int i = 1; i <= t; ++i)
    {
        x2 = mul(x1, x1, n); // x2 = x1^2 % n
        if (x2 == 1 && x1 != 1 && x1 != n - 1)
            return true; // 发现非平凡平方根 → n 是合数
        x1 = x2;         // 更新 x1 继续下一轮平方
    }

    // 经过 t 次平方后,最终值应为 1(若 n 是素数且 a 与 n 互质)
    if (x1 != 1)
        return true; // 最终值不为 1 → n 是合数
    return false;    // 本次测试未发现合数证据,n 可能是素数
}

bool miller_rabin(ll n)
{
    int s = 50; // 测试次数,可根据需要调整
    if (n < 2)
        return false; // 小于 2 的数不是素数
    if (n == 2)
        return true; // 2 是素数
    if (n % 2 == 0)
        return false; // 偶数(除 2 外)不是素数

    // 进行 s 次独立测试
    for (int i = 0; i < s; ++i)
    {
        // 随机选取底数 a,范围 [1, n-1]
        ll a = rand() % (n - 1) + 1;
        if (witness(a, n)) // 若 witness 返回 true,说明 n 是合数
            return false;
    }
    return true; // 所有测试通过,n 很可能是素数
}

素性检测能处理的最大变量long long约为1e19,再大就要自己处理高精度大数了,下面我们将深入探究如何大批量判断连续的一段数字中的素数的个数。

三、素数筛

1.引入

想想以下以下场景:给定你一个n,求2~n之间的所有素数。

不难发现,如果逐个判断,当n极大时,显然会超时,这时候就要用素数筛筛选所有整数,把非负数筛去。常见的素数筛有两种,埃氏筛和欧拉筛,欧拉筛的事件复杂度时O(n),达到最优,不可能更快了。

2.埃氏筛

埃氏筛(Eratosthenes Sieve)中文全称埃拉托斯特尼筛法,它用的是素数的定义,我们对初始2~n的队列(默认从小到大)进行筛选:

  1. 输出最小的素数2,然后筛掉所有2的倍数。
  2. 输出最小的素数3,然后筛掉所有3的倍数。
  3. 输出最小的素数5,然后筛掉所有5的倍数。
  4. ……

重复以上步骤,直到最后队列为空。

下面给出代码,方便读者理解:

int E_sieve(ll n)
{
    vector<bool> vis(n + 1, false);
    vector<int> prime;
    for (int i = 2; i <= sqrt(n); ++i) // 原理和试除法相似
    {
        if (!vis[i])
            for (int j = i * i; j <= n; j += i)
                vis[j] = true;
    }
    for (int i = 2; i <= n; ++i)
        if (!vis[i])
            prime[k++] = i;
    return prime.size();
}

要是读者想保存每一个素数可以将数组定义在全局里。相较于马上要介绍的欧拉筛,埃氏筛还是略有不足,比如12被2和3筛了两次,具体性能比较读者可自行查阅相关资料,在此不做详细分析。

3.欧拉筛

欧拉筛(Euler Sieve),是埃氏筛的改良版,达到了线性时间复杂度,因此又称线性筛。

欧拉筛的原理:一个合数肯定有一个最小质因数,让每个合数只被它的最小质因数筛选一次,以达到不重复筛选的目的。

具体原理如下:

  1. 逐一检查2~n的所有数,第一个检查的是2,它是第一个素数。
  2. 当检查到第i个数的时,利用已经求得的素数筛去对应的合数x,而且是用x的最小质因数去筛。

下面给出欧拉筛的代码:

vector<ll> euler_sieve(ll n)
{
    vector<bool> vis(n + 1, true); // 标记是否为素数,初始假设全是素数
    vector<ll> prime;              // 存储已发现的素数

    for (ll i = 2; i <= n; ++i)
    {
        if (vis[i])
        { // i 是素数,加入列表
            prime.push_back(i);
        }
        ll size = prime.size();
        // 用 i 乘以每个已发现的素数来筛除合数
        for (ll j = 0; j < size; ++j)
        {
            if (i > n / prime[j])
                break;
            vis[i * prime[j]] = false; // 标记为合数
            // 如果 primes[j] 是 i 的最小质因子,则停止循环
            if (i % prime[j] == 0)
                break;
        }
    }
    return prime;
}

虽然比赛可能用不到long long范围,但是于此知识为了敬告读者,警钟长鸣!!!

这里有一个小技巧,在一般算法竞赛环境中的空间限制下我们如果要使用静态数组,最大可以开到int prime[5800000]来储存N=1e8范围内的所有素数。

这里给到一个拓展,求1~n中每一个数的最小质因数,这里只需要简单的对标记数组修改变量类型为int然后每次储存的就是欧拉筛第一次筛到这个合数的素数值,对于素数,最小质因数就是它本身,这个问题就解决了,这里不给出代码,望读者能够自行思考,自己完成该代码然后去测试,如果还是有点不太理解我说的这一段话,那就代表还没明白欧拉筛的本质原理。

这一片是数论里素数的很重要的章节,望各位多动手勤加练习敲打,今天就到此结束了,good day!

文末附加内容
暂无评论

发送评论 编辑评论


				
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
(╯‵□′)╯︵┴─┴
 ̄﹃ ̄
(/ω\)
∠( ᐛ 」∠)_
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ`)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ( ̄∇ ̄o)
ヾ(´・ ・`。)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò。)
Σ(っ °Д °;)っ
( ,,´・ω・)ノ"(´っω・`。)
╮(╯▽╰)╭
o(*////▽////*)q
>﹏<
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
颜文字
Emoji
小恐龙
花!
上一篇
下一篇