Miller_Rabin素数测试与Pollard_Rho分解质因数

原创

BAJim·H 2017-11-08 21:48:20 ©著作权

文章标签 数论质数测试分解质因数费马小定理快速幂 文章分类 运维

©著作权归作者所有：来自51CTO博客作者BAJim·H的原创作品，请联系作者获取转载授权，否则将追究法律责任

Miller_Rabin测试

如果需要快速测试一个数是否是素数，有筛法与试除法

此处介绍的是一种基于费马小定理的不确定性算法，当然，这种算法的出错率是极其微小的，尤其当选择的测试数较多时，因此在OI比赛中成为一种实用的算法

好的进入正题。

我们知道费马小定理，若p为素数，那么对于a为任意正整数且a不为p的倍数时，有ap−1≡1(modp)

然而它的逆定理并不成立，即若ap−1≡1(modp)，p不一定是一个素数

我们称满足费马小定理的p为伪素数
它很大可能是一个素数

那么就有费马测试，多弄几个a，如果都符合，它就很有可能是素数

但是存在这样一些数，无论怎么取与这个数互质的a，它都满足费马小定理
定义这样的数叫做卡米切尔数
第一个卡米切尔数是561，=3*11*17

这样的数在1~108中有255个

这样的错误率是不能接受的

考虑强化这个算法。

我们知道若a2≡1(modp),(a,p)=1

那么(a+1)(a−1)≡0(modp)

由于(a,p)=1
那么要么a+1≡0(modp)，要么a−1≡0(modp)
那么同余方程的解为a1=1,a2=p−1

既然
ap−1≡1(modp)

考虑将p−1写成d∗2q的形式

那么就有了新的一种判断方法，将p-1不断除以2，判断它的结果是否仍是1或p-1
若是p-1则退出，若是1则继续除下去，直到除不了了或者出现了不是这两种结果的

这样的每除一次都需要跑一次快速幂，所以我们考虑优化。

可以一次直接除到底，然后再乘上来，看是否是出现了p-1或全部都是1

a可以取多几个，这样在1018内的错误率都是极小的，并且如果a随机，几乎是不可能卡掉的

附上代码

LL ti(LL x,LL y,LL m)
{
    if(y==0||x==0) return 0;
    LL s=ti(x,y/2,m);
    return (y%2)?((s+s)%m+x)%m:(s+s)%m;
}
bool miller_rabin(LL k)
{
    int d[5]={2,3,5,7,10007};
    fo(i,0,4)
    {
        if(k==d[i]) return 1;
        if(k%d[i]==0) return 0;
        LL t,n=k-1;
        while((n&1)==0) n>>=1;
        t=ksm(d[i],n,k);
        while(n<k-1&&t!=1&&t!=k-1) t=ti(t,t,k),n<<=1;//这里有可能两个小于10^18的数相乘会爆掉long long，所以要把乘改成类似快速幂的加法
        if(!((n&1)||t==k-1)) return 0;
    }   
    return 1;
}

Pollard_Rho算法

pollard_rho是一种基于miller_rabin测试的质因数分解算法

算法流程是若分解一个数k，先判断k是否是质数（Miller_rabin）
然后尝试找到k的一个非1非k的因子p，递归分解p和k/p

关键在于找这个因子p

首先介绍一个玄学的东西

生日“悖论”

一个班有k个同学，至少有两个同学的生日相同的概率P是多少？

k=1,P=0
k=2,P=1−1∗364365
k=3,P=1−1∗364365∗363365
。
。
。
令人难以置信的是
k=23,P=0.507
当k>55时，概率就很接近1了
当k=100,P=0.999999692751072

这与人们的感觉严重不符，因此称之为“悖论”。

事实上这个悖论可以扩展到一般情况
给出N个数，至少有一对数之差等于a的概率

这个增长同上面一样，也是非常快的。

回到正题
那么我们此时并不需要暴力随机因子，而是可以通过用两个数相减的方式计算因子。

定义一个随机函数f(x)=(x∗x+c)modk，c是这一次找因子随机出来的数

每次随机出新的x，与之前的x做差，与k求gcd，看这个数是否是合法的因子，不是则继续随机

这个算法有两点需要注意的地方。

之前有很多个x，取哪一个？
事实上，对于同一个c，x是会循环的

形如希腊字母ρ(rho)，又是pollard最先发明的，因此命名为pollard_rho算法

那么有两种做法。
一种是取一开始x=y，每次随机x=f(x),y=f(f(y))，y以两倍速度随机，如果最后x=y了，那么说明出现了环。

还有一种做法，我并不太知道这个做法的本质，但是据说这个做法非常快。

每次随机x=f(x)，在每第2的幂次随机时，更新y，即y=x，每次x与y做差，当某个时刻x=y了说明出现了环
这种做法具体可以看代码

LL ti(LL x,LL y,LL m)//乘法用类似快速幂的做法累加，防止爆long long
{
    if(y==0||x==0) return 0;
    LL s=ti(x,y/2,m);
    return (y%2)?((s+s)%m+x)%m:(s+s)%m;
}

LL get(LL k)
{
    LL c=rd(1e9)+1,r,x,y,i=1,n=2;
    x=y=rd(k-1)+1;
    while(1)
    {
        i++;
        x=(ti(x,x,k)+c)%k;//随机下一个x
        if(x==y) return k;//出现环
        r=gcd(abs(x-y),k);
        if(1<r&&r<k) 
        {
            return r;//找到因子
        }
        if(i==n) y=x,n<<=1;//n是2的幂，如果随机次数到了2的某次方则更新y
    }
}
void pollard_rho(LL k)
{
    if(k==1) return;
    if(miller_rabin(k))//这是判断是否是素数
    {
        pr[++pr[0]]=k;
        return;
    }
    LL p=k;
    while(p==k&&p!=1) p=get(k);
    pollard_rho(k/p),pollard_rho(p);
}