质数的判定

最新推荐文章于 2022-05-05 16:24:13 发布

逃夭丶

最新推荐文章于 2022-05-05 16:24:13 发布

阅读量298

点赞数

分类专栏： # 数论训练

本文链接：https://blog.csdn.net/qq_43724031/article/details/100145417

版权

数论训练专栏收录该内容

45 篇文章 0 订阅

订阅专栏

试除法

bool is_prime(int n)	//O(√n)
{
    if(n<2) return false;
    for(int i=2;i*i<=n;i++)
        if(n%i==0)  return false;
    return true;
}

下文参考于：https://www.cnblogs.com/p-z-y/p/10603061.html
下文参考于：https://www.cnblogs.com/forward777/p/10561003.html

Miller-Robbin算法（蒙特卡罗算法）

有极小的概率把合数错误判定为质数，但是多次判断合起来的错误概率趋近于零。

前置定理

费马小定理
· 内容：若 φ( p ) = p−1，p>1，则a^p ≡ a (mod p) 或 a^p−1 ≡ 1 (mod p)，(a<p)
· 证明：戳这里
二次探测定理：
· 内容：如果 φ( p )=p−1,p>1,p>X ，且X²≡1(modp)，那么X=1 or p−1
· 证明：
∵ X2≡1(modp)
∴ p|X2−1
∴ p|(X+1)(X−1)
∵ p是大于X的质数
∴ p=X+1 or p≡X−1(modp)，即X=1 or p−1。

算法原理

由费马小定理，我们可以有一个大胆的想法：满足 a^p−1 ≡ 1(mod p) 的数字 p 是一个质数。
可惜，这样的猜想是错误的，可以举出大量反例，如：2³⁴⁰ ≡ 1(mod 341)，然而 341 = 31∗11 。
所以，我们可以取不同的 a 多验证几次，不过，∀ a < 561,a⁵⁶⁰ ≡ 1 (mod 561)，然而 561 = 51∗11 。
这时，二次探测就有很大的用途了。结合费马小定理，正确率就相当高了。
这里推荐几个 ai 的值： 2，3，5，7，11，61，13，17。用了这几个 ai，就连那个被称为强伪素数的 46856248255981 都能被除去。
当 N<4,759,123,141选取 a=2,7,61即可确保算法得出正确结果。
当 N<3,825,123,056,546,413,051≈3∗10^18选取 a=2,3,5,7,11,13,17,19,23即可确保算法得出正确结果。
当 N<18,446,744,073,709,551,616=2^64选取 a=2,3,5,7,11,13,17,19,23,29,31,37即可确保算法得出正确结果

主要步骤

§ rabin函数（ bool rabin(ll x,ll y) ）：

首先利用费马小定理，排除 x^y-1 % y != 1的情况，因为该情况下，y一定不是质数，返回false；
当 x^y-1 % y == 1时，由上述反例，我们不能确定 y 一定是是质数，将 y-1 提取出所有 2 的因数，进行二次探测。

§miller函数（ bool miller(ll x) ）：
枚举所有a[i]的rabin结果（rabin(a[i],x)），得出 x 是否为质数

#define ll long long
#define ld long double
#define ull unsigned long long
ll a[] = {2,3,5,7,11,61,13,17};
ll mul(ll x,ll y,ll z)	//快速乘优化
{
    ll sm = (ld)x/z*y;
    return ((ull)x*y-(ull)sm*z+z)%z;
}
ll ksm(ll x,ll y,ll z)	//快速幂
{
  ll ans = 1;
  for(ll i=y;i;i>>=1){
        if(i&1) ans=mul(ans,x,z);
        x=mul(x,x,z);
    }
    return ans%z;
}
bool rabin(ll x,ll y)
{
    if(ksm(x,y-1,y)!=1)	//费马小定理
        return 0;
    ll z = y-1,sm;
    while(!(z&1)){
        z >>= 1;
        sm = ksm(x,z,y);
        if(sm!=1&&sm!=y-1)  return 0;	//二次探测
        if(sm==y-1) return 1;
    }
    return 1;
}
bool miller(ll x)
{
    bool ans = true;
    for(int i=0;i<8;i++){
        if(a[i]==x) return 1;
        ans &= rabin(a[i],x);
    }
    return ans;
}

Pollard Rho （大数质因数分解算法）

大致流程

先判断当前数是否是素数（这里就可应用 Miller−Robbin ），如果是则直接返回
如果不是素数的话，试图找到当前数的一个因子（可以不是质因子）
递归对该因子和约去这个因子的另一个因子进行分解

实现

这里所写的Pollard-Rho算法不是最朴素的版本，是一个更快更好的版本

不是每次算出下一个随机数之后都算gcd，而是把算的这些数都乘起来(当然要%一下x)
累计了一定量的数之后再求一次gcd，这样就大大减少了求gcd的次数从而提高速度
这里选定的是127个数累计起来求一次gcd，为什么是127呢，因为它是个好数字（我也不知道）
上面的优化有局限性，就是很有可能环比较小，没到127个数就出现环，这样即使已经出现过含x因数的数也会跳出循环。
遇到这样的情况就会拖慢速度，甚至永远都算不出来。
这里可以用一个倍增的方法解决这个问题，分别在生成(1,2,4,8,16,32,64…)个数的时候算一次gcd。

ll Pollard(ll x)
{
    ll a,b,d,g,y,i,j;
    while(1){
        a = b = rand()%x;
        g = rand()%x;
        y=1;i=0;j=1;
        while(++i){
            a = (mul(a,a,x)+g) % x;
            y = mul(y,abs(a-b),x);
            if(a==b||!y)    break;
            if(i<127||i==j){
                d=gcd(x,y);
                if(d>1&&d!=x)   return d;
                if(i==j){
                    b=a;
                    j<<=1;
                }
            }
        }
    }
}
void Find(ll x)	//找最大质因子
{
    if(x<=ans)  return;
    if(miller(x)){
        if(x>ans)   ans = x;	
        return;
    }
    ll y = Pollard(x);
    while(x%y==0)   x /= y;
    Find(y);Find(x);
}

例题：https://www.luogu.org/problem/P4718
AC代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<vector>
#include<stack>
#include<queue>
#include<ctime>
#define ll long long
#define ld long double
#define ull unsigned long long
#pragma GCC optimize("Ofast","inline","-ffast-math")
#pragma GCC target("avx,sse2,sse3,sse4,mmx")
using namespace std;
int t;
ll a[] = {2,3,5,7,11,61,13,17};
ll n,ans;
ll mul(ll x,ll y,ll z)
{
    ll sm=(ld)x/z*y;
    return ((ull)x*y-(ull)sm*z+z)%z;
}
ll ksm(ll x,ll y,ll z)
{
  ll ans=1;
  for(ll i=y;i;i>>=1){
        if(i&1) ans=mul(ans,x,z);
        x=mul(x,x,z);
    }
    return ans%z;
}
ll gcd(ll x,ll y)
{
    return y==0 ? x : gcd(y,x%y);
}
bool rabin(ll x,ll y)
{
    if(ksm(x,y-1,y)!=1)
        return 0;
    ll z=y-1,sm;
    while(!(z&1)){
        z>>=1;
        sm=ksm(x,z,y);
        if(sm!=1&&sm!=y-1)  return 0;
        if(sm==y-1) return 1;
    }
    return 1;
}
bool miller(ll x)
{
    bool ans = true;
    for(int i=0;i<8;i++){
        if(a[i]==x) return 1;
        ans &= rabin(a[i],x);
    }
    return ans;
}
ll Pollard(ll x)
{
    ll a,b,d,g,y,i,j;
    while(1){
        a = b = rand()%x;
        g = rand()%x;
        y=1;i=0;j=1;
        while(++i){
            a = (mul(a,a,x)+g) % x;
            y = mul(y,abs(a-b),x);
            if(a==b||!y)    break;
            if(i<127||i==j){
                d=gcd(x,y);
                if(d>1&&d!=x)   return d;
                if(i==j){
                    b=a;
                    j<<=1;
                }
            }
        }
    }
}
void Find(ll x)
{
    if(x<=ans)  return;
    if(miller(x)){
        if(x>ans)   ans = x;
        return;
    }
    ll y=Pollard(x);
    while(x%y==0)   x/=y;
    Find(y);Find(x);
}
int main()
{
    scanf("%d",&t);
    while(t--){
        scanf("%lld",&n);
        ans=1;
        Find(n);
        if(ans==n){
            printf("Prime\n");
            continue;
        }
        printf("%lld\n",ans);
    }
    return 0;
}