poj2429 GCD & LCM Inverse

用miller_rabin 和 pollard_rho对大数因式分解，再用dfs寻找答案即可。

http://poj.org/problem?id=2429

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef __int64 LL;
const int maxn = 100;
const int times = 10;
LL prime[maxn], k;
int cnt[maxn];
LL c, d, pm, mid;

void dfs(int pos, LL t){
    if(pos == k){
        if(pm < t) pm = t;
        return;
    }
    LL tem = 1;
    int c = cnt[pos];
    while(c--) tem *= prime[pos];
    dfs(pos + 1, t);
    if(t * tem <= mid) dfs(pos + 1, t * tem);
}

LL random(LL n){
    return (double)rand() / RAND_MAX * n + 0.5;
}

LL multi(LL a, LL b, LL mod){
    a %= mod, b %= mod;
    LL ans = 0;
    while(b){
        if(b & 1) ans += a, ans %= mod;
        b >>= 1;
        a <<= 1;
        a %= mod;
    }
    return ans;
}

LL power(LL a, LL p, LL mod){
    a %= mod;
    LL ans = 1;
    while(p){
        if(p & 1) ans = multi(ans, a, mod), ans %= mod;
        p >>= 1;
        a = multi(a, a, mod);
        a %= mod;
    }
    return ans;
}

LL gcd(LL a, LL b){
    if(!b) return a;
    return gcd(b, a % b);
}

bool witness(LL a, LL n){
    LL u = n - 1;
    while(!(u & 1)) u >>= 1;
    LL t = power(a, u, n);
    while(u != n - 1 && t != 1 && t != n - 1){
        t = multi(t, t, n);
        u <<= 1;
    }
    return t == n - 1 || u & 1;
}

bool miller_rabin(LL n){
    if(n == 2) return 1;
    if(n < 2 || !(n & 1)) return 0;
    for(int i = 0; i < times; i++){
        LL p = random(n - 2) + 1;
        if(!witness(p, n)) return 0;
    }
    return 1;
}

LL pollard_rho(LL n, LL t){
    LL x = random(n - 2) + 1;
    LL y = x, i = 1, k = 2, d;
    while(1){
        ++i;
        x = (multi(x, x, n) + t) % n;
        d = gcd(y - x, n);
        if(1 < d && d < n) return d;
        if(x == y) return n;
        if(i == k){
            y = x;
            k <<= 1;
        }
    }
}

void solve(){
    LL m = d / c;
    LL m1 = m;
    mid = (LL)sqrt(m);
    k = 0;
    memset(cnt, 0, sizeof cnt);
    if(m % 2 == 0){
        prime[k++] = 2;
        while(m % 2 == 0) m >>= 1, ++cnt[k - 1];
    }
    while(!miller_rabin(m) && m > 1){
        int seed = 12312;
        LL p1 = m;
        while(p1 >= m || !miller_rabin(p1)) p1 = pollard_rho(m, seed--);
        prime[k++] = p1;
        while(m % p1 == 0) m /= p1, ++cnt[k - 1];
    }
    if(m != 1) prime[k++] = m, ++cnt[k - 1];
    pm = 1;
    dfs(0, 1);
    LL qm = m1 / pm;
    printf("%I64d %I64d
", pm * c, qm * c);
}

int main(){
    //freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    while(~scanf("%I64d%I64d", &c, &d)) solve();
    return 0;
}