【poj 2429】GCD & LCM Inverse (Miller-Rabin素数测试和Pol

#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,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,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,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,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 %lldn",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 %lldn",rp * a,rq * a);
    }
}