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

Posted 2022-03-07 hs-zlq

本题涉及的算法个人无法完全理解，在此提供两个比较好的参考。

 

 原理

 代码实现

 

个人改编的AC代码：

#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);