【蒟蒻数学】素数相关算法

今天的网络赛再一次告诉我——数学真的很重要——队友推了两小时公式O(1)一发AC。所以今天就让我们来学习一下一些经典的素数相关算法吧!

素数

在程序设计竞赛中,我们很少碰到单纯关于素数的题目,但经常会在题目的某一部分用到关于素数的知识。而所谓素数,就是值自然数当中,除了1之外,只能被1和该数本身整除的数,下面我们就来一同学习求解素数的相关算法吧。

筛法求素数

朴素筛法求素数

我们先来看一下代码:

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#define SIZE 10000000
using namespace std;
int check[SIZE] = {0};//元素值为0代表是素数
int prime[SIZE] = {0};

int main() {
    int pos = 0;
    int flag;
    for (int i = 2 ; i < SIZE ; i++) {
        if (!check[i]) {
            //如果是素数
            prime[++pos] = i;
        }
        for (int j = 2 * i ; j <= SIZE ; j += i) {
            check[j] = 1;
        }
    }
    printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
    //for(int i = 1; i <= pos; i++)printf("%d ", prime[i]);
    return 0;
}

//Output:1.17

这种方法比较好理解,初始时,我们假设全部数都是素数,当找到一个素数时,即check[i]==0时,我们就把它加入素数表内,并且筛去掉它的倍数。显然这个素数乘上另外一个大于等于二的正整数之后得到的数——即它的倍数都是合数,把这些合数都筛掉,就首先了素数的筛选,这也是算法名字——“筛选法”的由来。而这种简单粗暴的方法在1e7的数据规模下平均用时约为1.17s。

其实就算是朴素的筛法也有优化:

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#define SIZE 10000000
using namespace std;
int check[SIZE] = {0};//元素值为0代表是素数
int prime[SIZE] = {0};

int main() {
    int pos = 0, num = int(sqrt(SIZE));
    int flag;
    for (int i = 2 ; i < SIZE ; i++) {
        if (!check[i]) {
            //如果是素数
            prime[++pos] = i;
        }
        if(i <= num)
        for (int j = i * i ; j <= SIZE ; j += i) {
            //i*i的优化快了不少 
            check[j] = 1;
        }
    }
    printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
    //for(int i = 1; i <= pos; i++)printf("%d ", prime[i]);
    return 0;
}

//Output:0.66

我们把删去素数倍数的操作从i2变成了ii,就快了近一倍,为什么可以这样优化呢?你想,我们假定i为大于2的一个素数,那么按照我们之前的筛法i2就在处理到素数2的时候筛了一遍,处理到素数i的时候又处理了一遍,造成了不必要的时间浪费,而改成ii,则可以从大于等于当前素数的倍数开始,从而避免了这个问题。

线性筛法求素数

我们先来看一下代码:

#include <ctime>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#define SIZE 10000000
using namespace std;
int prime[SIZE], isNotPrime[SIZE] = {1, 1};
int n, m, primeNum = 0, num;
int main() {
    for(int i = 2; i <= SIZE; i++) {
        if(!isNotPrime[i]) prime[++primeNum] = i;
        //key1 
        for(int j = 1; j <= primeNum && i * prime[j] <= SIZE; j++) {
            isNotPrime[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;//key2 
        }
    }
    printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
    //for(int i = 1; i <= primeNum; i++)printf("%d ", prime[i]); 
    return 0;
}
//Output:0.18

线性筛的最大优点就是每个数最多只被筛一次,也就是说时间复杂度是O(n)的。但其实线性筛素数的思想很简单,首先,我们要先明确一个条件,任何一个合数都能表示成一系列素数的积。

接着我们来看程序不管 i 是否是素数,都会执行到“key1”,这时我们分情况讨论:

①如果此时 i 是素数的话,我们进入j循环。那很简单,被筛选出来的合数就是一个大的素数 i 乘以不大于 i 的素数 prime[j],这样筛选出来的数跟之前的是不会重复的。此时筛出的数都是 N=p1p2的形式, p1,p2之间不相等且均为素数;例如35 = 5 7 只会在 i = 7 的时候被筛选出来;

②如果此时 i 是合数的话,我们进入j循环。此时 i 可以表示成递增素数相乘 i=p1p2...pn, pi都是素数(2<=i<=n),  pi<=pj  ( i<=j ),其中p1是最小的系数。根据“key2”的定义,当p1==prime[j] 的时候,筛除就终止了,也就是说,只能筛出不大于p1的质数i,这样筛出来的数也是之前没筛出来过的。我们可以直观地举个例子。i=235,此时能筛除 2i ,不能筛除 3i;如果能筛除3i 的话,当 i' 等于 i'=335 时,筛除2i' 就和前面重复了。

所以这样我们就在线性的时间里把素数筛选出来啦。

优化线性筛法求素数

我们先来看代码:

#include <ctime>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#define SIZE 10000000
using namespace std;
const int HALF = SIZE / 2; 
int prime[HALF] = {0, 2}, isNotPrime[HALF] = {1};
//isNotPrime[i] = 0 表示 i * 2 + 1 是素数 
int n, m, primeNum = 1, num;
int main() {
    for(int i = 1; i <= HALF; i++) {
        if(!isNotPrime[i]) prime[++primeNum] = 2 * i + 1;
        //key1 
        for(int j = 2; j <= primeNum && (2 * i + 1) * prime[j] <= SIZE; j++) {
            isNotPrime[((2 * i + 1) * prime[j] - 1) / 2] = 1;
            if((2 * i + 1) % prime[j] == 0) break;//key2 
        }
    }
    printf("%.2f", (double)clock()/CLOCKS_PER_SEC);
    //for(int i = 1; i <= primeNum; i++)printf("%d ", prime[i]); 
    return 0;
}
//Output:0.10

什么你觉得刚刚的线性筛还不够快?没事这里我想出了一种可以在约O(n/2)筛选出素数的线性筛。

其大致思路和普通的线性筛一样,只不过我们知道,除了2以外的偶数都是合数,所以这次我们只处理奇数。在key1后筛选合数的时候,偶数*质数 = 偶数 = 合数,所以我们直接去掉所有的偶数,只处理奇数也能筛选出全部的素数,这样我们就少处理了一半的数,实测运行时间也比之前的线性筛快了快一半。

素数的判定

朴素算法

素数的判定最原始的方法就是按照素数的定义进行判定,看看要判定的数存不存在除了1和它本身之外的其他因子:

bool isPrime(int x) {
    for(int i = 2; i < x; i++) {
        if(x % i == 0)return(false);
    }
    return(true);
}

改良算法

刚刚的那种算法时间复杂度是不太可观的。我们知道对于任意一个合数 x ,它一定可以分解成两个大于1的整数 a 和 b 的积,即 x = a * b;我们很容易发现,a 和 b 当中至少有一个数一定小于等于 x 的平方根。所以我们可以对算法进行如下改进:

bool isPrime(int x) {
    for(int i = 2; i * i <= x; i++) {
        if(x % i == 0)return(false);
    }
    return(true);
}

更优算法

在上面的算法中,其实我们也只需要找到 x 是否存在质因子即可,也就是说,从 2 到 x 的平方根的那些合数其实也是没有必要进行判断的。所以在要多次进行判断素数时,我们可以提前进行线性筛生成一个素数表,再进行判断,这样循环的次数还可以降低:

//int list[1100000] = {0, 2, 3, 5, ………………};
//如果需要判断的质数最大为 10^12, 只需要筛选√10^12 = 10^6 内的素数即可 
bool isPrime(int x) {
    for(int i = 1; list[i] * list[i] <= x; i++) {
        if(x % list[i] == 0)return(false);
    }
    return(true);
}

米勒-勒宾测试

上面的几种方法都是针对需要判定的数较小的情况,如果需要判定的数大到一定规模,上面的方法就不太适用了,这时可以考虑米勒-勒宾测试法,其描述如下:

令大奇数 n = 2^p * m + 1,其中 p 是非负整数,m 是正奇数。若 b^m ≡ 1(mod n)b^(2^j * m) ≡ 1(mod n),0 ≤ j ≤ p - 1,则称 n 通过以b为基的米勒-勒宾测试。若n为素数,b为整数,且 n 不能 整除 b,则 n 必然通过以b为基的米勒-勒宾测试。故代码如下:

#include<stdio.h>
#include<iostream>
#include<algorithm>
using namespace std;
long long mul(long long a, long long b, long long mod){//伪高精度乘法 
    a %= mod;
    b %= mod;
    long long c=(long double)a * b / mod;
    long long ans = a * b - c * mod;
    return (ans % mod + mod) % mod;
}
long long pow_mod(long long x, long long n, long long mod){//快速幂
    long long res=1;
    while(n){
        if(n & 1)
        res = mul(res, x, mod);
        x = mul(x, x, mod);
        n >>= 1;
    }
    return (res + mod) % mod;
}
bool Miller_Rabbin(long long a,long long n){

    //把n-1  转化成 (2^r)*s 
    long long s = n - 1, r = 0;
    while((s & 1) == 0){
        s >>= 1;
        r++;
    }

    //算出 2^d  存在 k 里
    long long k = pow_mod(a, s, n);

    //二次探测  看变化过程中是不是等于1 或 n-1
    if(k == 1)return true;
    for(int i = 0; i < r; i++,k = k * k % n){
        if(k == n - 1)return true;
    }
    return false;
}

而且据说按照下表取b的值  在一定范围内 可以使正确率达到百分之百,变成一个确定的算法:

if n < 1 373 653   b = 2 and 3.
if n < 9 080 191   b = 31 and 73.
if n < 4 759 123 141   b = 2, 7, and 61.
if n < 2 152 302 898 747      b = 2, 3, 5, 7, and 11.

结语

通过这篇BLOG,相信你对素数也有了一些初步的了解。最后希望你喜欢这篇BLOG!

Last modification:November 25th, 2019 at 05:24 pm
If you think my article is useful to you, please feel free to appreciate

Leave a Comment