浅谈 Miller-Robbin 与 Pollard Rho

时间  2019-11-12
标签 浅谈 miller robbin pollard rho 繁體版
原文   原文链接

前言

$Miller-Robbin$ 与 $Pollard Rho$ 虽然都是随机算法，不过用起来是真的爽。html

$Miller Rabin$ 算法是一种高效的质数判断方法。虽然是一种不肯定的质数判断法，可是在选择多种底数的状况下，正确率是能够接受的。算法

$Pollard Rho$ 是一个很是玄学的方式，用于在 $O(n^{1/4})$ 的指望时间复杂度内计算合数$n$的某个非平凡因子。优化

事实上算法导论给出的是 $O(\sqrt p)$ ， $p$ 是 $n$ 的某个最小因子，知足 $p$ 与 $\frac{n}{p}$ 互质。ui

可是这些都是指望，未必符合实际。但事实上 $Pollard Rho$ 算法在实际环境中运行的至关不错。spa

注：以上摘自洛谷。code

Miller-Robbin

前置芝士

1. 费马小定理htm

  • 内容：若 \(\varphi(p)=p-1,\,p>1\)，则\(a^{p}\equiv a\pmod{p}\) 或 \(a^{p-1}\equiv 1\pmod{p},\,(a<p)\)blog

  • 证实：戳这里递归

2. 二次探测定理get

  • 内容：若是 \(\varphi(p)=p-1,\,p>1,\,p>X\) ，且\(X^2\equiv 1\pmod{p}\)，那么\(X=1\ or\ p-1\)

  • 证实：

    \(\because\) \(X^2\equiv 1\pmod{p}\)

    \(\therefore\) \(p|X^2-1\)

    \(\therefore\) \(p|(X+1)(X-1)\)

    \(\because\) \(p\)是大于\(X\)的质数

    \(\therefore\) \(p=X+1\ or\ p\equiv X-1\pmod{p}\)，即\(X=1\ or\ p-1\)。

算法原理

由费马小定理，咱们能够有一个大胆的想法：知足 \(a^{p-1}\equiv 1\pmod{p}\) 的数字 \(p\) 是一个质数。

惋惜，这样的猜测是错误的，能够举出大量反例，如：\(2^{340}\equiv 1\pmod{341}\)，然鹅 \(341=31*11\) 。

因此，咱们能够取不一样的 \(a\) 多验证几回，不过，\(\forall a<561,\,a^{560}\equiv 1\pmod{561}\)，然鹅 \(561=51*11\) 。

这时，二次探测就有很大的用途了。结合费马小定理，正确率就至关高了。

这里推荐几个 \(a_i\) 的值： \(2，3，5，7，11，61，13，17\)。用了这几个 \(a_i\)，就连那个被称为强伪素数的 \(46856248255981\) 都能被除去。

主要步骤

  • .将 \(p-1\) 提取出全部 \(2\) 的因数，假设有\(s\) 个。设剩下的部分为 \(d\)（这里提取全部\(2\)的因数，是为了下面应用二次探测定理） 。

  • 枚举一个底数 \(a_i\) 并计算 \(x=a_i^{d}\pmod p\)。

  • 令 \(y=x^{2}\pmod p\)，若是没有出现过 \(p-1\)，那么 \(p\) 未经过二次探测，不是质数。

  • 不然，若底数已经足够，则跳出；不然回到第二步。

简易代码

#define ll long long
ll p,a[]={2,3,5,7,11,61,13,17};
inline ll mul(ll a,ll b,ll mod)
{
    ll ans=0;
    for(ll i=b;i;i>>=1)
    {
        if(i&1) ans=(ans+a)%p;
        a=(a<<1)%p;
    }
    return ans%p;
}
inline ll Pow(ll a,ll b,ll p)
{
    ll ans=1;
    for(ll i=b;i;i>>=1)
    {
        if(i&1) ans=mul(ans,a,p);
        a=mul(a,a,p);
    }
    return ans%p;
}
bool Miller_Robbin(ll p)
{
    if(p==2) return true;
    if(p==1 || !(p%2)) return false;
    ll d=p-1;int s=0;
    while(!(d&1)) d>>=1,++s;
    for(int i=0;i<=8 && a[i]<p;++i)
    {
        ll x=Pow(a[i],d,p),y=0;
        for(int j=0;j<s;++j)
        {
            y=mul(x,x,p);
            if(y==1 && x!=1 && x!=(p-1)) return false;
            x=y;
        }
        if(y!=1) return false;
    }
    return true;
}

Pollard Rho

大体流程

  • 先判断当前数是不是素数（这里就可应用 \(Miller-Robbin\) ），若是是则直接返回

  • 若是不是素数的话，试图找到当前数的一个因子（能够不是质因子）

  • 递归对该因子和约去这个因子的另外一个因子进行分解

如何找因子

一个一个试确定是不行的。而这个算法的发明者采起了一种清奇的思路。（即采起随机化算法）

  • 咱们假设要找的因子为p

  • 随机取一个 \(x、y\)，不断调整 \(x\) ，具体的办法一般是 \(x=x*x+c\)（c是随机的，也能够本身定）。

  • 取 \(p=gcd(y-x,n)\) ，若\(p \in \left(1,n\right)\) ，则找到了一个因子，就返回。

  • 若是出现 \(x=y\) 的循环，就说明出现了循环，并不断在这个环上生成之前生成过一次的数，因此咱们必须写点东西来判环：咱们能够用倍增的思想，让\(y\)记住\(x\)的位置，而后\(x\)再跑当前跑过次数的一倍的次数。这样不断让\(y\)记住\(x\)的位置，x再往下跑，由于倍增因此当\(x\)跑到\(y\)时，已经跑完一个圈。

  • 同时最开始设定两个执行次数\(i=一、k=2\)（即倍增的时候用） ，每次取 \(gcd\) 时 \(++i\) ；若是 \(i==k\) ，则令 \(y=x\) ，并将 \(k\) 翻倍。

完整代码

#include<cstdio>
#include<ctime>
#include<map>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
    TP f=1;x=0;register char c=getchar();
    for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
    for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
    x*=f;
}
template<typename TP>V print(TP x)
{
    if(x<0) x=-x,putchar('-');
    if(x>9) print(x/10);
    putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17,24251};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
template<typename TP>inline ll mul(TP a,TP b,TP p)
{
    ll ans=0;
    for(TP i=b;i;i>>=1)
    {
        if(i&1) ans=(ans+a)%p;
        a=(a<<1)%p;
    }
    return ans%p;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
    ll ans=1;
    for(TP i=b;i;i>>=1)
    {
        if(i&1) ans=mul(ans,a,p);
        a=mul(a,a,p);
    }
    return ans%p;
}
B Miller_Robbin(ll n)
{
    if(n<2) return false;
    if(n==2) return true;
    if(n%2==0) return false;
    ll d=n-1;int s=0;
    while(!(d&1)) d>>=1,++s;
    for(rg i=0;i<=8 && a[i]<n;++i)
    {
        ll x=Pow(a[i],d,n),y=0;
        F1(j,0,s-1)
        {
            y=mul(x,x,n);
            if(y==1 && x!=1 && x!=(n-1)) return false;
            x=y;
        }
        if(y!=1) return false;
    }
    return true;
}
inline ll Pollard_Rho(ll n)
{
    ll x,y,c,i,k;
    while(true)
    {
        ll x=rand()%(n-2)+1;
        ll y=rand()%(n-2)+1;
        ll c=rand()%(n-2)+1;
        i=0,k=1;
        while(++i)
        {
            x=(mul(x,x,n)+c)%n;
            if(x==y) break;
            ll d=Gcd(abs(y-x),n);
            if(d>1) return d;
            if(i==k) {y=x;k<<=1;}
        }
    }
}
V Find(ll n)
{
    if(n==1) return;
    if(Miller_Robbin(n)) {ans=max(ans,n);return;}
    ll p=n;
    while(n<=p) p=Pollard_Rho(p);
    Find(p);Find(n/p);
}
int main()
{
    read(T);srand(time(0));
    while(T--)
    {
        ans=0;
        read(n);Find(n);
        if(ans==n) puts("Prime");
        else print(ans),ed;
    }
    return 0;
}

emmmm \(\cdots\)


这数据也太毒瘤了吧！！


看来要疯狂卡常了

优化1（不如叫作卡常？）

蛋定的分析一波，咱们发现除了 $Pollard-Rho$ 是 $O(n^{1/4})$ 的指望时间复杂度外， $gcd$ 和龟速乘都是 $O(\log N)$ 的。

虽然这种复杂度已经很优秀了，可对于本题的数据（\(T≤350\) 、 \(1≤n≤10^{18}\)），仍是太 \(\cdots\)

因此咱们要果断摒弃这种很 $low$ 的龟速乘，改用一种暴力溢出的快速乘：

  • 简单原理： \(a \times b \ \ mod \ \ p=a \times b−\left \lfloor \frac{a \times b}{p} \right \rfloor\)

  • 用 long double 来处理这个 \(\left \lfloor \frac{a \times b}{p} \right \rfloor\)

  • 而后处理一下浮点偏差就能够了。

  • 模数较大时可能会出锅。

  • 不过出锅几率很小 \(\cdots\)

以下

inline ll mul(ll a,ll b,ll mod)
{
    a%=mod,b%=mod;
    ll c=(long double)a*b/mod;
    ll ret=a*b-c*mod;
    if(ret<0) ret+=mod;
    else if(ret>=mod) ret-=mod;
    return ret;
}

实践证实，战果辉煌：$6pts -> 94pts$ ！！！

优化2（正解）

虽然关于龟速乘的 $O(\log n)$ 的恶劣影响获得了必定遏制，不过，我仍是好想 $AC$ 啊！

经过办法1 ：打表 \(\cdots\)

正确 $AC$ 姿式以下：

  • 咱们发如今 \(Pollard-Rho\) 中若是长时间随机化而得不到结果，\(gcd\)带来的 \(O(\log n)\) 仍是很伤肾的！！那有没有办法优化呢？答案是确定的。

  • 在生成\(x\)的操做中，龟速乘所模的数就是\(n\)，而要求的就是\(n\)的某一个约数，即如今的模数并非一个质数。

  • 根据取模的性质：若是模数和被模的数都含有一个公约数，那么此次模运算的结果必然也会是这个公约数的倍数。因此若是咱们将若干个\((y−x)\) 相乘，由于模数是 \(n\) ，因此若是若干个\((y−x)\)中有一个与\(n\)有公约数，最后的结果定然也会含有这个公约数。

  • 因此能够多算几回\((y−x)\)的乘积再来求\(gcd\) （通常连续算\(127\)次再求一次\(gcd\)）。

  • 不过\(\cdots\)，若是在不断尝试\(x\)的值时碰上一个环，就可能会还没算到\(127\)次就跳出这个环了，就没法得出答案；同时，可能\(127\)次计算以后，全部\((y−x)\)的乘积都变成了\(n\)的倍数（即\(\prod_{i=1}^{127} {(y-x)} \equiv 0 \pmod{n}\) ）

  • 因此咱们能够不只在每计算\(127\)次以后求\(gcd\)、还要在倍增时（即判环时）求\(gcd\)，这样既维护了其正确性，又判了环！！

  • 完整\(AC\)代码：

#include<cstdio>
#include<ctime>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
    TP f=1;x=0;register char c=getchar();
    for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
    for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
    x*=f;
}
template<typename TP>V print(TP x)
{
    if(x<0) x=-x,putchar('-');
    if(x>9) print(x/10);
    putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
inline ll mul(ll a,ll b,ll mod)
{
    a%=mod,b%=mod;
    ll c=(long double)a*b/mod;
    ll ret=a*b-c*mod;
    if(ret<0) ret+=mod;
    else if(ret>=mod) ret-=mod;
    return ret;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
    ll ans=1;
    for(TP i=b;i;i>>=1)
    {
        if(i&1) ans=mul(ans,a,p);
        a=mul(a,a,p);
    }
    return ans%p;
}
B Miller_Robbin(ll n)
{
    if(n<2) return false;
    if(n==2) return true;
    if(n%2==0) return false;
    ll d=n-1;int s=0;
    while(!(d&1)) d>>=1,++s;
    for(rg i=0;i<=8 && a[i]<n;++i)
    {
        ll x=Pow(a[i],d,n),y=0;
        F1(j,0,s-1)
        {
            y=mul(x,x,n);
            if(y==1 && x!=1 && x!=(n-1)) return false;
            x=y;
        }
        if(y!=1) return false;
    }
    return true;
}
inline ll Pollard_Rho(ll n)
{
    while(true)
    {
        ll x=rand()%(n-2)+1;
        ll y=rand()%(n-2)+1;
        ll c=rand()%(n-2)+1;
        ll i=0,k=1,b=1;
        while(++i)
        {
            x=(mul(x,x,n)+c)%n;
            b=mul(b,abs(y-x),n);
            if(x==y || !b) break;
            if(!(i%127) || i==k)
            {
                ll d=Gcd(b,n);
                if(d>1) return d;
                if(i==k) y=x,k<<=1;
            }
        }
    }
}
V Find(ll n)
{
    if(n<=ans) return;
    if(Miller_Robbin(n)) {ans=max(ans,n);return;}
    ll p=Pollard_Rho(n);
    while(n%p==0) n/=p;
    Find(p),Find(n);
}
int main()
{
    read(T);srand(time(0));
    while(T--)
    {
        ans=0;
        read(n);Find(n);
        if(ans==n) puts("Prime");
        else print(ans),ed;
    }
    return 0;
}
相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程