poj2429 GCD & LCM Inverse

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

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

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <algorithm>
  4 #include <cmath>
  5 using namespace std;
  6 typedef __int64 LL;
  7 const int maxn = 100;
  8 const int times = 10;
  9 LL prime[maxn], k;
 10 int cnt[maxn];
 11 LL c, d, pm, mid;
 12
 13 void dfs(int pos, LL t){
 14     if(pos == k){
 15         if(pm < t) pm = t;
 16         return;
 17     }
 18     LL tem = 1;
 19     int c = cnt[pos];
 20     while(c--) tem *= prime[pos];
 21     dfs(pos + 1, t);
 22     if(t * tem <= mid) dfs(pos + 1, t * tem);
 23 }
 24
 25 LL random(LL n){
 26     return (double)rand() / RAND_MAX * n + 0.5;
 27 }
 28
 29 LL multi(LL a, LL b, LL mod){
 30     a %= mod, b %= mod;
 31     LL ans = 0;
 32     while(b){
 33         if(b & 1) ans += a, ans %= mod;
 34         b >>= 1;
 35         a <<= 1;
 36         a %= mod;
 37     }
 38     return ans;
 39 }
 40
 41 LL power(LL a, LL p, LL mod){
 42     a %= mod;
 43     LL ans = 1;
 44     while(p){
 45         if(p & 1) ans = multi(ans, a, mod), ans %= mod;
 46         p >>= 1;
 47         a = multi(a, a, mod);
 48         a %= mod;
 49     }
 50     return ans;
 51 }
 52
 53 LL gcd(LL a, LL b){
 54     if(!b) return a;
 55     return gcd(b, a % b);
 56 }
 57
 58 bool witness(LL a, LL n){
 59     LL u = n - 1;
 60     while(!(u & 1)) u >>= 1;
 61     LL t = power(a, u, n);
 62     while(u != n - 1 && t != 1 && t != n - 1){
 63         t = multi(t, t, n);
 64         u <<= 1;
 65     }
 66     return t == n - 1 || u & 1;
 67 }
 68
 69 bool miller_rabin(LL n){
 70     if(n == 2) return 1;
 71     if(n < 2 || !(n & 1)) return 0;
 72     for(int i = 0; i < times; i++){
 73         LL p = random(n - 2) + 1;
 74         if(!witness(p, n)) return 0;
 75     }
 76     return 1;
 77 }
 78
 79 LL pollard_rho(LL n, LL t){
 80     LL x = random(n - 2) + 1;
 81     LL y = x, i = 1, k = 2, d;
 82     while(1){
 83         ++i;
 84         x = (multi(x, x, n) + t) % n;
 85         d = gcd(y - x, n);
 86         if(1 < d && d < n) return d;
 87         if(x == y) return n;
 88         if(i == k){
 89             y = x;
 90             k <<= 1;
 91         }
 92     }
 93 }
 94
 95 void solve(){
 96     LL m = d / c;
 97     LL m1 = m;
 98     mid = (LL)sqrt(m);
 99     k = 0;
100     memset(cnt, 0, sizeof cnt);
101     if(m % 2 == 0){
102         prime[k++] = 2;
103         while(m % 2 == 0) m >>= 1, ++cnt[k - 1];
104     }
105     while(!miller_rabin(m) && m > 1){
106         int seed = 12312;
107         LL p1 = m;
108         while(p1 >= m || !miller_rabin(p1)) p1 = pollard_rho(m, seed--);
109         prime[k++] = p1;
110         while(m % p1 == 0) m /= p1, ++cnt[k - 1];
111     }
112     if(m != 1) prime[k++] = m, ++cnt[k - 1];
113     pm = 1;
114     dfs(0, 1);
115     LL qm = m1 / pm;
116     printf("%I64d %I64d\n", pm * c, qm * c);
117 }
118
119 int main(){
120     //freopen("in.txt", "r", stdin);
121     //freopen("out.txt", "w", stdout);
122     while(~scanf("%I64d%I64d", &c, &d)) solve();
123     return 0;
124 }