1. 程式人生 > >大數因數分解pollard rho

大數因數分解pollard rho

  1. 快速乘&快速冪
  2. 費馬小定理
  3. 二次探測定理
  4. pollor rho
  5. miller rabin

在介紹pollard rho演算法之前,先普及一下快速乘法及快速冪,因為大數的乘法之類的可能會爆long long;

#include<cstdio>
#include<iostream>

using namespace std;

long long p,n;

int mul(long b,long a)//快速乘
{
    long long ret=0;
    while(b)
    {
        if(b&1)
        {
            ret=(ret+a)%p;
        }
        b>>=1
; a=(a+a)%p; } return ret; } int pow(long long a,long long n)//快速冪 { long long ret=1; while(n) { if(n&1) { ret=mul(ret,a); } n>>=1; a=mul(a,a); } return ret; } int main() { while(scanf("%I64d%I64d",&n,&p)) { printf
("%I64d",pow(2,n)%p); } }

費馬小定理

假如p是質數,且gcd(a,p)=1,那麼 a^(p-1)≡1(mod p)。
即:假如a是整數,p是質數,且a,p互質(即兩者只有一個公約數1),那麼a的 (p-1)次方除以p的餘數恆等於1。

證明

引理1.
  若a,b,c為任意3個整數,m為正整數,且(m,c)=1,則當ac≡bc(modm)時,有a≡b(modm)

證明:ac≡bc(mod m)可得ac–bc≡0(mod m)可得(a-b)c≡0(mod m)因為(m,c)=1即m,c互質,c可以約去,a– b≡0(mod m)可得a≡b(mod m)

引理2.
  設m是一個整數,且m>1,b是一個整數且(m,b)=1.如果a1,a2,a3,a4,…am是模m的一個完全剩餘系,則ba[1],ba[2],ba[3],ba[4],…ba[m]也構成模m的一個完全剩餘系.

證明:若存在2個整數ba和ba[j]同餘即ba≡ba[j](mod m),根據引理1則有a≡a[j](mod m).根據完全剩餘系的定義可知這是不可能的,因此不存在2個整數ba和ba[j]同餘.所以ba[1],ba[2],ba[3],ba[4],…ba[m]構成模m的一個完全剩餘系.

構造素數 p 的完全剩餘系
{1,2,3,4……,p-1}
因為 gcd(a,p)==1,由引理2可得
{a,2a,3a,4a,5a……(p-1)a}也為完全剩餘系;

則有完全剩餘系性質可得1*2*3*4*5…..*(p-1)==a*2a*3a*4a*5a……(p-1)a%p;
所以(p-1)!≡(p-1)!*a^(p-1) (mod p);
a^(p-1)≡1(mod p)

二次探測定理:

如果n是一個素數,且n>x>0,則方程x²%n=1的解為:x=1或 x=n-1.
證明:
除x=1這個解外,剩下的x滿足 n>x>√n , n²>x²>n;
那麼則有x²=n+1 2n+1 3n+1……(n-2)n+1 …(n-1)n+1
由於n是素數,所以n不可分解,所以只有(n-2)n+1可化為平方形式,即x=n-1;

pollard rho演算法帶有一定的概率性,但也算是比較靠譜的一種板子了,錯判率大概只有4^(-s),s值是自定的;
該演算法主要是將大數不斷分解為小數;
費馬小定理偽質數推斷可知,如果n是一個正整數,如果存在和n互素的正整數a滿足 a^(n-1)≡ 1(mod n),我們說n是基於a的偽素數。如果一個數是偽素數,那麼它幾乎肯定是素數,多次查詢,則誤判機率能降到極低;
所以不斷構造隨機數x1,找到這樣一個因子p,使p=gcd(x1-x2,n);若p==1則構造失敗,基於x1不斷調整x2,否則p就是a的一個因數,遞迴繼續找p,及n/p;

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<ctime>
using namespace std;
long long factor[1000005];
int tot;
const int S=50;

long long mult(long long a,long long b,long long c)
{
    a%=c;
    b%=c;
    long long ret=0;
    while(b)
    {
        if(b&1)
        {
            ret=(ret+a)%c;
        }
        b>>=1;
        a<<=1;
        if(a>=c)a%=c;
    }
    return ret;
}

long long pow(long long x,long long n,long long mod)
{
    if(n==1)return x%mod;
    x%=mod;
    long long tmp=x;
    long long ret=1;
    while(n)
    {
        if(n&1)
        {
            ret=mult(ret,tmp,mod);
        }
        tmp=mult(tmp,tmp,mod);
        n>>=1;
    }
    return ret;
}

bool check(long long a,long long n,long long x,long long t)
{
    long long ret=pow(a,x,n);
    long long last=ret;
    for(int i=1;i<=t;i++)
    {
        ret=mult(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1)return true;//是合數
        last=ret;
    }
    if(ret!=1)return true;
    return false;
}

bool miller_rabin(long long n)//判素數
{
    if(n<2)return false;
    if(n==2) return true;
    if((n&1)==0)return false;
    long long x=n-1;
    long long t=0;
    while((x&1)==0)
    {
        x>>=1;
        t++;
    }
    for(int i=0;i<S;i++)
    {
        long long a=rand()%(n-1)+1;
        if(check(a,n,x,t))//如果檢查出來是合數
        return false;
    }
    return true;
}

long long gcd(long long a,long long b)
{
    if(a==0)return 1;
    if(a<0)return gcd(-a,b);
    while(b)
    {
        long long t=a%b;
        a=b;
        b=t;
    }
    return a;
}

long long pollard_rho(long x,long long c)
{
    long long i=1,k=2;
    long long x0=rand()%x;
    long long y=x0;
    while(1)
    {
        i++;
        x0=(mult(x0,x0,x)+c)%x;
        long long d=gcd(y-x0,x);
        if(d!=1&&d!=x)return d;
        if(y==x0)return x;
        if(i==k)
        {
            y=x0;
            k+=k;
        }
    }
}
    void findphi(long long n)
    {
        if(miller_rabin(n))
        {
            factor[tot++]=n;
            return;
        }
        long long p=n;
        while(p>=n)
        {
            p=pollard_rho(p,rand()%(n-1)+1);

        }
        findphi(p);
        findphi(n/p);
    }

    int main()
    {
        long long n;
        while(scanf("%I64d",&n)!=EOF)
        {
            tot=0;
            findphi(n);
            for(int i=0;i<tot;i++)
            printf("%I64d",factor[i]),printf("\n");
            if(miller_rabin(n))printf("yes\n");
            else printf("no\n");
        }
        return 0;
    }