【poj 2429】GCD & LCM Inverse (Miller-Rabin素数测试和Pollard_Rho_因数分解)

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
using namespace std;
#define ll long long

int top;
const int S = 10;
ll sta[101], Mul[101];

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

ll mul(ll a, ll b, ll mod) {
    ll r = 0;
    while (b) {
        if (b & 1) {
            r = r - mod + a;
            if (r < 0) r += mod;
            //r = (r + a) % mod;
        }
        a = a - mod + a;
        if (a < 0) a += mod;
        //a = (a + a) % mod;
        b >>= 1;
    }
    return r;
}  // 64位数相乘

ll ksm(ll a, ll n, ll mod) {
    ll r = 1;
    while (n) {
        if (n & 1) r = mul(r, a, mod);
        a = mul(a, a, mod);
        n >>= 1;
    }
    return r;
}  // 64位数乘方

bool isprime(ll n) {
    if (n == 2) return true;
    if (n < 2 || (n & 1) == 0) return false;
    ll a, x, y, u = n - 1;
    int t = 0;
    while ((u & 1) == 0) {
        t++;
        u >>= 1;
    }
    for (int i = 0; i < S; i++) {
        a = rand() % (n - 1) + 1;  //[1,n-1]
        x = ksm(a, u, n);
        for (int j = 0; j < t; j++) {
            y = mul(x, x, n);
            if (y == 1 && x != 1 && x != n - 1) return false;
            x = y;
        }
        if (x != 1) return false;
    }
    return true;
}

ll Abs(ll x) {
    if (x >= 0) return x;
    return -x;
}

void rho(ll n) {  //注意:不能处理1，否则会运行到对n-1=0取模
    if (isprime(n)) {
        for (int i = 1; i <= top; i++) {
            if (n == sta[i]) {
                Mul[i] *= n;
                return;
            }
        }
        top++;
        sta[top] = Mul[top] = n;
        return;
    }
    ll x, y, z, c, d;
    while (true) {
        x = y = rand() * rand() % (n - 1) + 1;
        c = rand() * rand() % (n - 1) + 1;  // c!=0&&c!=-2
        for (int i = 2, j = 2;; i++) {
            x = mul(x, x, n) - n + c;
            if (x < 0) x += n;  // 64位数相加
            d = gcd(Abs(x - y), n);
            if (d > 1 && d < n) {
                rho(d);
                rho(n / d);
                return;
            }
            if (x == y) break;
            if (i == j) y = x, j <<= 1;
        }
    }
}

int main() {
    ll a, b;
    while (scanf("%lld%lld", &a, &b) != EOF) {
        top = 0;
        memset(sta, 0, sizeof(sta));
        if (b / a == 1) {
            printf("%lld %lld\n", a, a);
            continue;
        }
        rho(b / a);
        ll p, q, res = 0, rp, rq;
        for (int i = 0; i < 1 << top; i++) {
            p = q = 1;
            for (int j = 0; j < top; j++) {
                if (i & (1 << j))
                    p *= Mul[j + 1];
                else
                    q *= Mul[j + 1];
            }
            if (p + q <= res || res == 0) {
                res = p + q;
                rp = p;
                rq = q;
            }
        }
        if (rp > rq) swap(rp, rq);
        printf("%lld %lld\n", rp * a, rq * a);
    }
}