Some algorithms for prime judgment (summary && contrast)

Source: Internet
Author: User
Tags fread log log cmath

Primality testing is one of the most common techniques in number theoretic problems. It can be very basic or very advanced ( philosophical ). This time, we mainly introduce the artifice of the prime judgment.

The judgment of prime number is divided into two main types: range Screening && single judging type

Let's start with the usual starting point of the scope filter type, which is tested using a template problem Luogu P3383 "template" linear sieve primes.

1. Ed Sieve

This is the most commonly used sieve, and the idea is simple: any multiples of a prime number are composite

Then we O (n) sweep again, while sifting the multiples of prime numbers

But some numbers, like 6, will be screened by 2 and 3, resulting in a waste of efficiency, so the complexity is proven to be **o (n log log n)

CODE

#include<cstdio>using namespace std;const int N=10000005;bool vis[N];int n,m,x;inline char tc(void){    static char fl[100000],*A=fl,*B=fl;    return A==B&&(B=(A=fl)+fread(fl,1,100000,stdin),A==B)?EOF:*A++;}inline void read(int &x){    x=0; char ch=tc();    while (ch<‘0‘||ch>‘9‘) ch=tc();    while (ch>=‘0‘&&ch<=‘9‘) x=x*10+ch-‘0‘,ch=tc();}inline void get_prime(int m){    register int i,j;    for (vis[1]=1,i=2;i<=m;++i)    if (!vis[i]) for (j=i<<1;j<=m;j+=i) vis[j]=1;}int main(){    //freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);    read(n); read(m); get_prime(n);    while (m--)    {        read(x);         puts(vis[x]?"No":"Yes");    }    return 0;}
2. Linear sieve (Euler sieve)

This is actually the optimization of the upper, we realize that a number should be only its minimum factorization deleted , so we can be the same as the number of screens simultaneously recorded prime, this is the real O (n) complexity

CODE

#include<cstdio>using namespace std;const int N=10000005;int prime[N],n,m,x,cnt;bool vis[N];inline char tc(void){    static char fl[100000],*A=fl,*B=fl;    return A==B&&(B=(A=fl)+fread(fl,1,100000,stdin),A==B)?EOF:*A++;}inline void read(int &x){    x=0; char ch=tc();    while (ch<‘0‘||ch>‘9‘) ch=tc();    while (ch>=‘0‘&&ch<=‘9‘) x=x*10+ch-‘0‘,ch=tc();}inline void Euler(int n){    register int i,j;    for (vis[1]=1,i=2;i<=n;++i)    {        if (!vis[i]) prime[++cnt]=i;        for (j=1;j<=cnt&&i*prime[j]<=n;++j)        {            vis[i*prime[j]]=1;            if (!(i%prime[j])) break;        }    }}int main(){    //freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);    read(n); read(m);    Euler(n);    while (m--)    {        read(x);        puts(vis[x]?"No":"Yes");    }    return 0;}

Note the above sentence:

if (!(i%prime[j])) break;

This ensures that the efficiency of the linear sieve does not result in duplication, because when i%prime[j]==0 the number is the number of subsequent deletions.

3. Basic primality test

This is the most basic method of prime number determination. From 2 to sqrt (X) enumeration whether there are numbers divisible by x

The proof is simple, because if the number is a prime, then its factor must be 1 and X, if its factor is greater than sqrt (x), then the square is greater than X, which is obviously not possible.

So we O (sqrt (x)) Judge once

CODE

#include<cstdio>#include<cmath>using namespace std;int n,m,x;inline char tc(void){    static char fl[100000],*A=fl,*B=fl;    return A==B&&(B=(A=fl)+fread(fl,1,100000,stdin),A==B)?EOF:*A++;}inline void read(int &x){    x=0; char ch=tc();    while (ch<‘0‘||ch>‘9‘) ch=tc();    while (ch>=‘0‘&&ch<=‘9‘) x=x*10+ch-‘0‘,ch=tc();}inline bool check(int x){    if (!(x^1)) return 0;    register int i; int bound=(int)sqrt(x);    for (i=2;i<=bound;++i)    if (!(x%i)) return 0;    return 1;}int main(){    //freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);    read(n); read(m);    while (m--)    {        read(x);        puts(check(x)?"Yes":"No");    }    return 0;}
4. Optimization of Algorithm 3

First we look at a conclusion:

Prime numbers greater than or equal to 5 must be adjacent to multiples of 6.

Proof and other references: Dalao ' s blog

And then the same 3, we just fast into 6 units each time, and then the constant gets the indescribable optimization (the fastest algorithm to become the problem)

CODE

#include<cstdio>#include<cmath>using namespace std;int n,m,x;inline char tc(void){    static char fl[100000],*A=fl,*B=fl;    return A==B&&(B=(A=fl)+fread(fl,1,100000,stdin),A==B)?EOF:*A++;}inline void read(int &x){    x=0; char ch=tc();    while (ch<‘0‘||ch>‘9‘) ch=tc();    while (ch>=‘0‘&&ch<=‘9‘) x=x*10+ch-‘0‘,ch=tc();}inline bool check(int x){    if (!(x^1)) return 0;    if (!(x^2)||!(x^3)) return 1;    if ((x%6)^1&&(x%6)^5) return 0;    register int i; int bound=(int)sqrt(x);    for (i=5;i<=bound;i+=6)    if (!(x%i)||!(x%(i+2))) return 0;    return 1;}int main(){    //freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);    read(n); read(m);    while (m--)    {        read(x);        puts(check(x)?"Yes":"No");    }    return 0;}
5.miller-rabin algorithm

This is the fastest way to judge the number of primes in history (but the algorithm 4 amputated in this topic)

First, the algorithm is based on the Fermat theorem and the two-time detection theorem:

Two-time detection theorem: If P is an odd prime, then the solution for X2≡1 (MODP) is x = 1 or x = p-1 (mod p)

So we can turn X into the R*2^t form, where R is an odd

Then we combine two algorithms && fast power to stabilize O (log x) for a single judgment

But this algorithm is a non-perfect algorithm , every time it is 25% probability is wrong, so we can choose a few numbers to get a few times

But by chance I see a paragraph on the internet:

For the primality judgment of large numbers, the Miller-rabin algorithm is widely used at present. The general base is still randomly selected, but when the number to be measured is not too large, choosing a test base has some tricks. For example, if the number of tests is less than 4 759 123 141, then only the three base 2, 7 and 61 will be tested. Of course, the more you test, the greater the right range is sure to be. If you test each time with the first 7 primes (2, 3, 5, 7, 11, 13, and 17), all numbers not exceeding 341 550 071 728 320 are correct. If you use 2, 3, 7, 61, and 24251 as the base, then the only strong pseudo prime in 10^16 is 46 856 248 255 981. Some of these conclusions make the Miller-rabin algorithm very useful in Oi. It is generally believed that the correct rate of the Miller-rabin primality test is acceptable, and the error rate of a randomly selected K-base test algorithm is probably 4^ (-K).

So for this topic n=10000000 the scope of only need to choose 2,7,61 can

CODE

#include <cstdio>using namespace Std;typedef long long ll;const int prime[3]={2,7,61};int n,m,x;inline char tc (void    ) {static Char fl[100000],*a=fl,*b=fl; Return a==b&& (b= (A=FL) +fread (Fl,1,100000,stdin), a==b)? eof:*a++;}    inline void read (int &x) {x=0; char ch=tc (); while (ch< ' 0 ' | |    Ch> ' 9 ') CH=TC (); while (ch>= ' 0 ' &&ch<= ' 9 ') x=x*10+ch-' 0 ', CH=TC ();}    inline int quick_pow (int x,int p,int mod) {int tot=1;        while (p) {if (p&1) tot= ((LL) tot*x)%mod; x= (LL) x*x)%mod;    p>>=1; } return tot;} inline bool Miller_rabin (int x) {if (!) (    x^2)) return 1; if (x<2| |! (x&1))    return 0;    int t=0,u=x-1; while (!) (    u&1)) ++t,u>>=1; for (register int i=0;i<3;++i) {if (!) (        X^prime[i])) return 1; if (! (        X%prime[i])) return 0;        int Lst=quick_pow (PRIME[I],U,X);            for (register int j=1;j<=t;++j) {int now= (LL) lst*lst)%x; if (! ( NOW^1) &&lst^1&&lst^ (x-1)) return 0;        Lst=now;    } if (lst^1) return 0; } return 1;}    int main () {//freopen ("code.in", "R", stdin), Freopen ("Code.out", "w", stdout); Read (n);    Read (m);        while (m--) {read (x); Puts (Miller_rabin (x)? "    Yes ":" No "); } return 0;}

Finally, the results of 5 algorithms are given (no O2).

    1. Ed Sieve

    2. Linear sieve (Euler sieve)

    3. Basic primality testing

    4. Optimized algorithm 3

    5. Miller-rabin

Some algorithms for prime judgment (summary && contrast)

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.