LG的数学计划----分解质因数(Pollard-Rho算法)

1.对于我们朴素的求解质因数, 暴枚真是个好算法

好吧一样的就不给出代码了,

2.对于另一种神奇的算法Pollard-Rho算法

随机化算法, 与Miller robin有着密切联系, 可以先看一看两种算法都不难, 只是很神奇。。
这个算法也是一步步分解, 再递归求解, 但是它实现分解的过程非常奇妙, 我们如果在p中rand()分解, 这无疑是很愚蠢的。。。但是p又一定可以分解, 所以我们要while(1),
接下来是一个神奇的结论,
上面说过一个个试还不如暴枚, 那么我们就要优化,
取两个数, x1, x2, 看gcd(x2-x1, p)是否大于1, 是的话递归求解,
那么现在随中的几率更高了, 但是对于大数还是很慢, 所以我们要再优化,
那么对于x1, x2的求法就要讲究一下, 我们可以发现
采用这个求法

x2 = (x1*x1+c) % p;
这样是不是快了很多, 因为x2-x1的差重复的几率减小了
但是还是可能出现很背的情况, 你不换个c这破电脑不给你出来, 那我们就要想办法解决这组种子找不出解跳出来的情况,
我们发现 ,这样算会求出一个%p的循环节, 正是我们所需要的,它会走出一个神奇的圈!
就像这样

我们的目的是找到这个重复显然一个个存下来很麻烦, 我们用一种倍增的方法来解决这个问题(正确性我不会证感觉这样就好了), 每次*2当一次倍增完成后, 覆盖一次x1, x2
代码

#include <ctime>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

#define rep(i, s, t) for(int i = s; i <= t; ++i)
typedef long long ll;

ll H;
ll pls(ll a, ll b, ll p) {
    ll res = 0;
    for(; b; b >>= 1, a = (a+a)%p)
        if(b & 1) res = (res+a)%p;
    return res%p;
}

ll pow(ll a, ll b, ll p) {
    ll res = 1;
    for(; b; b >>= 1, a = pls(a, a, p))
        if(b & 1) res = pls(res, a, p)%p;
    return res % p;
}

bool M(ll p) {
    ll x[60] = {0}, s = 20;
    ll rest = p-1, t = 0;
    while(!(rest%2)) {
        t ++;
        rest >>= 1;
    }

    while(s--) {
        ll a = rand()%(p-1)+1;
        x[0] = pow(a, rest, p);
        rep(i, 1, t) {
            x[i] = pow(x[i-1], 2, p)%p;
            if(x[i] == 1) 
                if((x[i-1] != 1) && (x[i-1] != p-1)) return false;
        }
        if(x[t] ^ 1) return false;
    }
    return true;
}

ll gcd(ll a, ll b) {
    while(b) {
        ll t = a%b;
        a = b;
        b = t;
    }
    return a;
}

ll P(ll p) {
    ll c = rand()%(p-1)+1, x1 = rand()%(p-1)+1, x2 = x1;
    for(ll i = 2, k = 2; true; ++i) {
        x1 = (pls(x1, x1, p) + c)%p;
        ll G = gcd(p, (x2-x1+p)%p);
        if(G > 1 && G < p) return G;
        if(x2 == x1) return p;
        if(i == k) x2 = x1, k <<= 1;
    }
}

void solve(long long n) {
    if(n == 1) return;
    if(M(n)) {
        H = min(H , n);
        return;
    }
    long long p = n;
    while(p == n) p = P(p);
    solve(p);
    solve(n/p);
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("data.in", "r", stdin);
    freopen("res.out", "w", stdout);
#endif
    int _;
    scanf("%d", &_);
    while(_--) {
        ll p;
        H = 1LL << 54;
        scanf("%lld", &p);
        solve(p);
        if(H ^ p)
            printf("%lld
", H);
        else puts("Prime");
    }
    return 0;
}

其实这是一个裸题代码poj1811, 会写一篇code解释

原文地址:https://www.cnblogs.com/pbvrvnq/p/8530160.html