zcmimi's blog
avatar
zcmimi
2020-08-05 10:38:00
  • 本文总阅读量

Miller Rabin

判断一个数是否为素数O(\log n)

前置知识:

  1. 费马小定理

    p为素数,a^{p-1}\equiv 1 \pmod p

  2. 二次探测定理

    p为质数x^2\equiv 1 \pmod p,那么x\equiv 1\pmod px\equiv -1\equiv p-1\pmod p

p为素数,那么p-1可以分解为2^k\cdot t的形式

随机取一个数a,先算出a^t,然后再不断自乘,如果当前数\not \equiv \pm 1 \pmod p而自乘后的数\equiv 1 \pmod p,那么p不是素数

看上去有些玄学

p通过一次测试,则p不是素数的概率为25\%

那么测试x次,则p不是素数的概率为\displaystyle \frac 1{4^x}

我们可以保证非常高的正确率

HDU 2138 How many prime numbers

#include<cstdio>
int pw(int x,int b,int p){
    int res=1;
    while(b){
        if(b&1)res=1ll*res*x%p;
        b>>=1;x=1ll*x*x%p;
    }
    return res;
}
bool chk(int p){
    if(p==2)return 1;
    if(p<2||!(p&1))return 0;
    static int test[]={2,3,5,7,11,13,17,19,23,29,31,37,41};
    int t=p-1,k=0;
    while(!(t&1))++k,t>>=1;
    for(int i=0;i<4&&test[i]<p;++i){
        int a=pw(test[i],t,p),nxt;
        for(int j=1;j<=k;++j){
            nxt=1ll*a*a%p;
            if(nxt==1&&a!=1&&a!=p-1)return 0;
            a=nxt;
        }
        if(a!=1)return 0;
    }
    return 1;
}
int main(){
    int T;
    while(~scanf("%d",&T)){
        int x,ans=0;
        while(T--)scanf("%d",&x),ans+=chk(x);
        printf("%d\n",ans);
    }
}

Pollard-Rho

Pollard Rho是一个著名的大数质因数分解算法

Pollard Rho可以在O(n^{\frac 14})内找出n的一个因子

我们需要在n个数中猜出n的一个因数,成功率\dfrac 1n

我们知道因数是成对的,我们可以在[1,\sqrt n]内猜,成功率\dfrac 1{\sqrt n}

但对于n\ge10^{18},这个成功率还是不够

生日攻击

考虑生日攻击_维基,生日攻击_百度

若老师选择特定日期(例如9月16日),则至少有一名学生在那天出生的概率是\displaystyle 1-(364/365)^{30},约为7.9\%。但是,与我们的直觉相反的是,至少一名学生和另外任意一名学生有着相同生日的概率大约为70.63\%n = 30时)

生日悖论的实质是:由于采用"组合随机采样"的方法,满足答案的组合比单个个体要多一些.这样可以提高正确率.

Pollard Rho的随机函数

我们不放选取一组数x_1,x_2,x_3,\dots,x_nx,若有\gcd(|x_i-x_j|,N)>1,则称\gcd(|x_i-x_j|,N)N的一个因子

考虑这么一个序列x_i={x_{i-1}}^2+c \pmod N,其中c是一个随机的常数.

这样生成的数列为伪随机数

可以证明期望n^{\frac 14}个数后会开始循环

判环

极少数情况下,该算法甚至在遍历整个数列后仍然无法找到解

我们需要判环并更换c后重新尝试

4718 【模板】Pollard-Rho算法

n是质数,则直接返回

Pollard Rho找出n的一个因子p,除去n中所有p得到n',继续分析pn,可进行适当剪枝

考虑路径倍增优化

由于一直取\gcd会让算法变得很慢,我们可以通过乘法累积来减少求gcd的次数,若\gcd(a,b)>1,那么\gcd(ac,b)>1,我们可以考虑将一段路径乘起来\pmod n,再做一次\gcd

我们不妨考虑倍增的思想:每次在路径\rho上倍增地取一段[2^{k-1},2^k]的区间.将2^{k-1}记为l,2^k记为r.取而代之的,我们每次取的\gcd测试样本为|x_i-x_l|,其中l \le i \le r.我们每次积累的样本个数就是2^{k-1}个,是倍增的了.这样可以在一定程度上避免在环上停留过久,或者取gcd的次数过繁的弊端.

我们选择2^7=127将"迭代7次"作为上界是实验得出的较优方案

#include<bits/stdc++.h>
typedef long long ll;
typedef __int128 i128;
ll pw(ll x,ll b,const ll&p){
    ll res=1;
    while(b){
        if(b&1)res=(i128)res*x%p;
        b>>=1;x=(i128)x*x%p;
    }
    return res;
}
bool chk(const ll&p){
    if(p==2)return 1;
    if(!(p&1)||p<2)return 0;
    ll t=p-1;int k=0;
    while(!(t&1))t>>=1,++k;
    static int test[]={2,3,5,7,11,13,17,19,23,29};
    for(int i=0;i<10&&test[i]<p;++i){
        ll a=pw(test[i],t,p),nxt;
        for(int j=1;j<=k;++j){
            nxt=(i128)a*a%p;
            if(nxt==1&&a!=1&&a!=p-1)return 0;
            a=nxt;
        }
        if(a!=1)return 0;
    }
    return 1;
}
ll mxf;
ll f(const ll&x,const ll&c,const ll&n){
    return ((i128)x*x+c)%n;
}
ll abs(const ll&x){return x<0?-x:x;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll pr(const ll&x){
    ll s=0,t=0,c=1ll*rand()%(x-1)+1;
    int stp=0,goal=1;
    ll v=1;
    for(goal=1;;goal<<=1,s=t,v=1){
        for(stp=1;stp<=goal;++stp){
            t=f(t,c,x);
            v=(i128)v*abs(t-s)%x;
            if(stp%127==0){
                ll d=gcd(v,x);
                if(d>1)return d;
            }
        }
        ll d=gcd(v,x);
        if(d>1)return d;
    }
}
void fact(ll x){
    if(x<=mxf||x<2)return;
    if(chk(x)){
        mxf=mxf>x?mxf:x;
        return;
    }
    ll p=x;
    while(p>=x)p=pr(x);
    while(x%p==0)x/=p;
    fact(x),fact(p);
}
int main(){
    int T;scanf("%d",&T);
    while(T--){
        srand(time(0));
        ll n;scanf("%lld",&n);
        mxf=0;
        fact(n);
        if(mxf==n)printf("Prime\n");
        else printf("%lld\n",mxf);
    }
}

参考:

Miller Rabin算法与Pollard Rho算法
comment评论
Search
search