超浅谈素数

bbbzzx · 2021-11-10 20:38:59 · 算法·理论

前言

一些关于素数的东西，有误请指正。

Update

2021-114-514：开坑。

2024-9-26：增加前言，改名关于素数，增加欧拉筛。

2025-2-6：乱改。

2025-2-7：又学了一坨。

2025-2-10：改名超浅谈素数。

定义

素数（prime number），又称质数，其定义为只能被 1 和自己整除的正整数，素数集用 P 表示。

性质

存在无穷多个素数。
素数分布随整数增大越来越稀疏，随机整数 x 是素数的概率为 \frac{1}{\log_2 x}（证明：波特兰定理）。
素数定理：定义 \pi(x) 小于等于 x 的素数个数，随 x 的增长，\pi(x) 与 \frac{x}{\ln x} 的比值趋于 1。
威尔逊定理：若 p 为素数，有 (p-1) ! \equiv -1\pmod p，非素数则不成立。
素性检验

试除法

朴素算法，有一些小技巧，只需判断到 \sqrt{n} 即可，时间复杂度为 O(\sqrt{n})。
```
bool isprime(int n)
{
if(n==1) return 0;
for(int i=2;i*i<=n;i++)
    if(n%i==0)
        return 0;
return 1;
}
```
费马素性检验
费马小定理：设 p \in P，对于所有正整数 a，且 \gcd(a,p)=1，都有 a^{p-1} \equiv 1\pmod p

其逆定理几乎成立，而费马素性检验就是基于费马小定理的逆定理的随机概率法素性测试。

过程：

对于正整数 n，随机选取正整数 a,1 < a \leq n 。

若不满足费马小定理，则一定不是素数。
若满足，则很有可能是素数，继续尝试其他 a 值，尝试次数越多则可能性越大，称 n 为基于 a 的可能素数。

缺陷：存在一类合数，对于所有满足条件的 a 都能满足条件，称这类数为 Carmichael 数，其分布密度很小，前三个为 561、1105、1729。

设尝试次数为 k，此算法的时间复杂度为 O(k \log n)，在 n \geq 10^{12} 时比试除法要快得多。

bool Fermat(int n)
{
    if(n<3) return n==2;
    for(int i=1;i<=8;i++){
        int a=rand()%(n-2)+2;
        if(qpow(a,n-1,n)!=1) return 0;
    }
    return 1;
}

Miller-Rabin 素性检验

显然 Carmichael 数的存在使得费马测试的正确性大打折扣，而 Miller-Rabin 素性检验则可以排除 Carmichael 数。

二次探测定理：对于奇素数 p ，方程 x^2 \equiv1 \pmod p(x \in [1,p]) 仅有两个解 x_1=1 和 x_2=p-1。

证明：

移项可得 x^2-1 \equiv 0 \pmod {p}，即 (x+1)(x-1) \equiv 0 \pmod p，于是有 x \equiv 1 \pmod p 或 x \equiv p-1 \pmod p。

Miller-Rabin 算法就应用了这个定理。

步骤：

对于正整数 n，随机选取 a,1 < a \leq n，若 a|n 则 n 不是素数，反之则用费马小定理与二次探测定理结合起来验证：

设 n-1=u \times 2^t，于是只需验证 (a^{u})^{2^t}，设 v=a^u \mod n，然后对 v 做 t 次二次探测，若 v^2 \mod n=1，而 v \neq n-1 足以证明 n 是合数；若不存在这种情况，认为 n 是一个可能素数。

此算法时间复杂度为 O(k \log^3 n)。

bool Miller_Rabin(int n)
{
    if(n==1)
        return 0;
    int k=0;
    int p=n-1;
    while(!(p&1))
    {
        p>>=1;
        k++;
    }
    int t=8;
    for(int i=0;i<t;i++)
    {
        int a=rand()%(n-1)+1;
        int b=qpow(a,p,n);
        bool flag=0;
        if(b==1)
            continue;
        for(int j=0;j<k;j++)
        {
            if((b+1)%n==0)
            {
                flag=1;
                break;
            }
            else
                b=(b*b)%n;
        }
        if(flag)
            continue;
        return 0;
    }
    return 1;
}

筛法

对于一段连续区间的素数查找，考虑使用筛法来求。

埃氏筛

Eratosthenes 筛，主要的思想就是 1 不是素数，然后再从 2 开始，先把第一个数标记为素数，这个数的倍数都标记为合数。

时间复杂度为 O(n \log\log n)。

证明：

观察过程，易得筛的次数为：

\sum_{i=1}^{\pi(n)} \lfloor \frac{n}{p_i} \rfloor

=O(n \sum_{i=2}^{\pi(n) } \frac{1}{p_i} )

=O(n \int_{2}^{\pi(n)} \frac{\mathrm dx}{p_i})

由素数定理，\pi(n)=\Theta(\frac{n}{\ln n})，有 \pi(p_i)=i=\Theta(\frac{p_i}{\ln p_i})，于是 p_n=O(n \ln n)。

所以，

O(n \int_{2}^{\pi(n)} \frac{\mathrm dx}{p_i})

=O(n \int_{1}^{n}\frac{\mathrm dx}{x \ln x})

=O(n \int_{1}^{n} \frac{\mathrm d\ln x}{lnx})

=O(n \ln \ln n)

void getprime(int n)
{
    for(int i=2;i<=n;i++)//1不是素数，所以从2开始
    {
        if(isprime[i]==false)
        {
            prime[++cnt]=i;//先标记为素数
            for(int j=i;i*j<=n;j++)
                isprime[i*j]=true;//倍数标记为合数
        }
    }
}

欧拉筛

对埃氏筛进行改进，注意到在埃氏筛中，一个合数被筛了多次，欧拉筛的改进就是使每个合数只被筛一次，那么用哪个因子筛呢？显然是用最小素因子。

时间复杂度为线性的 O(n)，因此得名线性筛。

void Prime()
{
    for(int i=2;i<=n;i++)
    {
        if(!a[i])
            prime.push_back(i);
        for(int j=0;j<prime.size();j++)
        {
            if(i*prime[j]>=n)//判越界
                break;
            a[i*prime[j]]=1;//用最小素因子去筛
            if(i%prime[j]==0)//若不是最小素因子，说明已经被筛过了
                break;
        }
    }
}

线性筛求 \pi(x)

void getprime(int n)
{
    memset(isprime,1,sizeof isprime);
    isprime[1]=0;
    pi[1]=0;
    for(int i=2;i<=n;i++)
    {
        if(isprime[i])
        {
            prime[++cnt]=i;
            pi[i]=cnt;  
        }
        else
            pi[i]=pi[i-1];
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
        {
            isprime[i*prime[j]]=0;
            if(i%prime[j]==0)
                break;
        }
    }
}

质因数分解

整数唯一分解定理（算术基本定理）：\forall n \in N_+，都有 n=\prod_{p_i \in P,c_i \in N_+} p_i^{c_i}。

线性筛求最小素因子

只需将筛数组用于记录最小素因子即可。

void getprime(int n)
{
    memset(vis,0,sizeof vis);
    int cnt=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prime[++cnt]=i;
            vis[i]=i;   
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
        {
            vis[i*prime[j]]=prime[j];
            if(i%prime[j]==0)
                break;
        }
    }
}

试除法

同样，考虑暴力。

对于 n，求其质因数可以先找到最小素因子，于是考虑用试除法找最小素因子 p，然后令 n_1=\frac{n}{p}，继续找 n_1，根据算术基本定理，素因子单调不降，所以时间复杂度为 O(\sqrt n)，比较低效。

void factor(int n){
    int p[25],c[25];
    int cnt=0;
    for(int i=2;i<=sqrt(n);i++){
        if(n%i==0){
            p[++cnt]=i,c[cnt]=0;
            while(n%i==0) n/=i,c[cnt]++;
        }
    }
    if(n>1) p[++cnt]=n,c[cnt]=1;
}

Pollard-Rho

Pollard-Rho 是一种启发式随机化算法。

考虑找出 n 的一个因子 p，然后就可以递归找 p 和 \frac{n}{p}。

考虑如何找 p，试除法对 [2,\sqrt n] 内的数一个个试除导致效率低下，为 O( \sqrt n)，考虑用随机数寻找，说不定运气好能碰上一个因子。

对于 n，考虑构造一个数列 \{x_n\}，有以下递推关系 x_i=(x_{i-1}^2+c) \mod n，而其中的 x_1 和 c 则用随机数生成初始值。

按此规律迭代下去，在后面会形成一个循环节，很好证明，由于在模 n 意义下，产生的数列取值一定是在 [1,n-1]，这是个有限的范围，于是在往复迭代下总会产生重复的取值，于是就有了循环节。

生日悖论：

假设每年都是 365 天，问至少多少人中，存在两个人生日在同一天的概率超过 50 \% ？

设有 k 个人，根据概率公式，有以下不等式成立：

\prod_{i=0}^{k-1} \frac{n-i}{n} \leq \frac{1}{2}

因为 1+x \leq e^x ，

于是

\prod_{i=0}^{k-1} e^{-\frac{i}{n}} \leq \frac{1}{2}

\frac{k(k-1)}{2n} \geq \ln 2

解得 k=O(\sqrt n)。

生日悖论带给我们一个启发，随机选取一列数，出现重复数字的抽样规模期望是 O(\sqrt n)，于是 \{x_i\} 循环节的期望长度为 O(\sqrt n)。

再考虑构造另一个数列 \{x_i \mod p\}，记为 \{x_i' \}，有，

x_i'=(x_{i-1}^2+c)\mod n\mod p

由于 p|n，原式又等于：

x_i'=(x_{i-1}^2+c)\mod p

所以它的循环节期望长度为 O(\sqrt p)。

考虑在数列 \{x_i'\} 中找到因子，如果存在 x_i' \equiv x_j' \pmod p，那么就存在因子 p=\gcd(|x_i-x_j|,n)，但是如果出现 p=n 的情况说明 x_i=x_j，此时应该重新随机初始值。

考虑判断循环节以确定寻找因子的次数上限，判环有两种方式，Floyd 和 Brent，Floyd 判环的思想是让 j=2i，判断 x_i 是否等于 x_j；Brent判环基于倍增思想，令 y=x_{2^k}，依次判断 x_{2^k+1} 到 x_{2^{k+1}} 是否与 y 相等。

不使用优化的期望时间复杂度为 O(\sqrt p \log n)=O(n^{\frac{1}{4}} \log n)。

基于 Brent 判环的 Pollard_Rho:

inline int pollard_rho(int n){
    int i=1,k=2;
    int c=rand()%(n-1)+1;
    int x=rand()%n;
    int y=x;
    while(1){
        i++;
        x=(qmul(x,x,n)+c)%n;
        int d=gcd(abs(y-x),n);
        if(d!=1&&d!=n) return d;
        if(y==x) return n;
        if(i==k) y=x,k<<=1;
    }
}
void fac(int n){
    if(n<=max_factor||n==1) return;
    if(miller_rabin(n)){
        max_factor=max(max_factor,n);
        return;
    }
    int p=n;
    while(p>=n) p=pollard_rho(p);
    while(n%p==0) n/=p;
    fac(p); fac(n);
}

这个复杂度只能算小优化，我们考虑把 log 消掉，虽然可以用 Stein 算法优化 gcd，但是收效甚微，考虑使用倍增优化 pollard_rho 函数，而 Brent 判环为我们提供了倍增优化的条件，我们考虑把每次 |y-x_i| 乘起来以减少求 gcd 的次数。

证明：

如果 \gcd(a,n)>1，那么 \gcd(ab \mod n,n)>1。

所以考虑只计算 \gcd(\prod |y-x_i| \mod n,n)，若其大于 1 说明成功找到因子，若出现 \prod |y-x_i| \mod n=0，则相当于判环了。

那我们考虑取几个 |y-x_i| 做一次 gcd，设取 k 个，期望时间复杂度为 O(\sqrt p+k^{-1}\sqrt p \log n)，要使这两项同阶，则要使 k 和 \log n 同阶，这就可以用 Brent 判环的倍增优化，于是我们成功把复杂度降到了 O(n^{\frac{1}{4}})。

使用这个优化可以过模板题，实测很快。

Code

#include<bits/stdc++.h>
#define int __int128
#define inf 1e18
#define db double
#define gc getchar
#define pc putchar
#define rg register
#define IOS ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
using namespace std;
inline int read(){
    int x=0,f=1;
    char ch=gc();
    while(!isdigit(ch)){
        if(ch=='-') f=-f;
        ch=gc();
    }
    while(isdigit(ch)) x=(x<<3)+(x<<1)+ch-'0',ch=gc();
    return x*f;
}
inline void write(int x){
    if(x<0) pc('-'),x=-x;
    if(x>9) write(x/10);
    pc(x%10+'0');
}
#define ctz __builtin_ctzll
int max_factor;
inline int gcd(int a, int b){
    if(!a||!b) return a|b;
    if(a==1||b==1) return 1;
    int az=ctz(a),bz=ctz(b),z=min(az,bz),t; b>>=bz;
    while(a) a>>=az,t=a-b,az=ctz(t),b=min(a,b),a=t<0?-t:t;
    return b<<z;
}
inline int qpow(int a, int b, int mod){
    int ans=1;
    while(b){
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
inline bool miller_rabin(int n)
{
    if(n==1) return 0;
    if(n==2||n==3) return 1;
    int k=0,p=n-1;
    while(!(p&1)) p>>=1,k++;
    for(int i=0;i<8;i++){
        int a=rand()%(n-1)+1;
        int b=qpow(a,p,n);
        if(b==1) continue;
        int j;
        for(j=0;j<k;j++){
            if(b==n-1) break;
            else b=b*b%n;
        }
        if(j==k) return 0;
    }
    return 1;
}
inline int pollard_rho(int n){
    int i=1,k=2;
    int c=rand()%(n-1)+1;
    int x=rand()%n;
    int y=x,s=1;
    while(1){
        i++;
        x=(x*x%n+c)%n,s=s*(y-x<0?x-y:y-x)%n;
        if(y==x||!s) return n;
        if(i==k){
            int d=gcd(s,n);
            if(d>1) return d;
            y=x,k<<=1;
        }
    }
}
void fac(int n){
    if(n<=max_factor||n==1) return;
    if(miller_rabin(n)){
        max_factor=max(max_factor,n);
        return;
    }
    int p=n;
    while(p>=n) p=pollard_rho(p);
    while(n%p==0) n/=p;
    fac(p); fac(n);
}
signed main(){
    int T=read();
    while(T--){
        srand(11451919);
        max_factor=0;
        int n=read(); 
        fac(n);
        if(max_factor==n) puts("Prime");
        else{
            write(max_factor);
            puts("");
        }
    }
    return 0;
}

End