Pallord rho大整数分解

Pallord rho算法

算法原理

这里简要阐述一下，我们这里要分解的整数是n，我们将存下来的质因数存在fact数组中。这里我们在1~(n - 1)当中选取一个随机数x，然后定义了一个函数g(x) = x * x + c(c一般是1，为了得到不同结果可以不断改变c的值)，如果gcd(x - g(x), n) = d 的结果不是1也不是n的话那么d就是n的一个约数。那么这个数列不断计算出来就会出现一个长度为l的循环节，当出现循环节的时候我们就我们就在当前c的参数下没有找到这个约数。所以我们计算的时候可以按照这个思路，当我们计算gcd((x - y + n) % n, n)的时候可以按照这个顺序

x g(2) g(3) g(4) g(5) g(6) g(7) g(8) g(9) g(10) g(11)
y g(1) g(2) g(2) g(4) g(4) g(4) g(4) g(8) g(8)  g(8)

这样我们看到x和y序号的差值不断在增长，所以我们总会找到这个循环节，那么我们就要在找到循环节之前找到一对x， y使得计算出我们的约数。以上的方法一般就是我们常用的。

最后fact里面就是质因数，但是并不是有序的，需要再排序。

里面涉及到Miller_Rabin素数测试，详见Miller_Rabin素数测试

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<cstdlib>
#include<vector>
#define ll __int64
#define MAX 5500

const int S = 20;//素数测试误判率=1 / 2^S   S越大误差率越小

using namespace std;


ll gcd(ll a, ll b){
    return b == 0 ? a : gcd(b, a % b);
}


ll multi_mod(ll a, ll b, ll p){
    ll ret = 0, q = a % p;
    while(b){
        if(b & 1){
            ret = (ret + q) % p;
        }
        q = (q + q) % p;
        b >>= 1;
    }
    return ret % p;
}

ll pow_mod(ll a, ll b, ll p){
    ll ret = 1, q = a % p;
    while(b){
        if(b & 1){
            ret = multi_mod(ret, q, p);
        }
        q = multi_mod(q, q, p);
        b >>= 1;
    }
    return ret % p;
}

bool Miller_Rabin(ll p){
    int i, j, k;
    ll u, t, x, y;
    if(p == 2)  return true;
    if(p % 2 == 0 || p == 1) return false;
    u = p - 1, t = 0;
    while(u % 2 == 0){
        t++;
        u >>= 1;
    }
    for(int i = 0; i < S; i++){
        x = rand() % (p - 1) + 1;
        x = pow_mod(x, u, p);
        ll y = 0;
        for(j = 0; j < t; j++){
            y = multi_mod(x, x, p);
            if(y == 1 && x != 1 && x != p - 1){
                return false;
            }
            x = y;
        }
        if(y != 1)  return false;
    }
    return true;
}



ll fact[MAX];//质因数数组
ll tot = 0;

ll pollard_rho(ll n, ll c){
    ll i = 1, k = 2;
    ll x = rand() % (n - 1) + 1;
    ll y = x;
    while(true)
    {
        i++;
        x = (multi_mod(x, x, n) + c) % n;
        ll d = gcd((y - x + n) % n, n);
        if(1 < d && d < n) return d;
        if(y == x) return n;
        if(i == k)
        {
            y = x;
            k <<= 1;
        }
    }
}

void find(ll n, ll c){
    if(n == 1){
        return;
    }
    if(Miller_Rabin(n)){
        fact[tot++] = n;
        return;
    }
    ll d = n;
    ll k = c;
    while(d >= n){
        d = pollard_rho(n, c--);
    }
    find(d, k);
    find(n / d, k);
}