<模板>(Miller-Rabin和Pollard_rho算法)poj 2429 GCD & LCM Inverse （数论）

本文出自http://blog.sina.com.cn/s/blog_69c3f0410100uac0.html

ps：我就把这个题当作(Miller-Rabin和Pollard_rho）算法的模板来储存了，回头再细学这两个算法。

题意：给出两个整数 m、n ( m =< n < 2^63)，求出两个整数 a、b ，使他们的最大公约数为 m ,最小公倍数为 n ，如果有多组，取 a + b 值最小的那一组。
思路：由题意，则令 p = a / m , q = b / m , s = n / m,则有GCD(p , q) = 1 ,LCM(p , q) = p * q = s 。
由以上结论可知此题是让我们把 s 分解成两个互质数相乘的形式。我们可以把 pollard_rho 和 miller_rabin 算法相结合将 s 分解成素因子相乘的形式，由于GCD(p , q) = 1，所以多个相同的素数不能同时分到 p 和 q 里，因此，我们多个相同的素数当作一个整体处理。又由于 s 在 2^63 之内，不会超过14个不同的素数相乘，所以可以将这些数枚举相乘，找到最接近s的平方根且不大于 s 平方根的数即为 p ，q = s / p.
例如 m = 2, n = 120,则 s = 60 = 2 * 2 * 3 * 5 = 4 * 3 * 5 ,从 3 、4 、5中挑出任意个数相乘得到 3 , 4 , 5 , 12 , 15 , 20 , 60 ，其中最靠近根号s且不大于根号s的就是 5 ，即 p =5 , q = s / p = 12 ,所以最后结果为 a = p * m = 10 , b = q * m = 24 

Get新知识：（1）关于文中p,q的lcm的求解。

  （2）大数的分解和大素数判定。

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
__int64 pr[100],f[100],t,ans;
__int64 mult(__int64 a,__int64 b,__int64 n)// a*b % n
{
    __int64 s=0;
    while(b){
        if(b&1)s=(s+a)%n;
        a<<=1;b>>=1;
        if(a>n)a=a%n;
    }
    return s;
}
__int64 miller_rabin(__int64 n,__int64 i)//素数测定，错误机率为1/4
{
    __int64 j,k,b=0,a=1;
    j=(n-1)/(1<<i);
    while(b%n==0)b=abs(rand())%n;
    while(j){
        if(j&1)a=mult(a,b,n);
        b=mult(b,b,n);j>>=1;
    }
    for(k=1;k<=i;k++){
        if(k<i-1&&mult(a,a,n)==1&&a!=1&&a!=n-1)return 0;
        a=mult(a,a,n);
    }
    if(a==1)return 1;
    return 0;
}
__int64 prime(__int64 n)
{//调用miller_rabin素数测试10次,正确概率为1 - (1/4)^10,几乎为1
    __int64 m,k,i,s=1;
    if(n==1)return 0;
    for(k=0,m=n-1;m%2==0;m/=2,k++);
    for(i=0,s=1;i<10;i++)s*=miller_rabin(n,k);
    return s;
}
__int64 gcd(__int64 a,__int64 b)//求最大公约数
{
    return b?gcd(b,a%b):a;
}
__int64 pollard_rho(__int64 n)pollard_rho随机算法求n的一个约数
{
    __int64 i=1,a,b,p,t=2,c=0;
    a=b=abs(rand())%(n-1)+1;
    while(c==0||c==2)c=abs(rand())%(n-1)+1;
    do{
        i++;
        p=gcd(n+b-a,n);
        if(p>1&&p<n)return p;
        if(i==t){
            b=a;t<<=1;
        }
        a=(mult(a,a,n)+n-c)%n;
    }while(a!=b);
    return n;
}
void add(__int64 n)//将素因子添加到素因子表中
{
    __int64 i;
    for(i=0;i<t;i++)
        if(pr[i]==n)break;
    if(i<t)f[i]++;//第i个素因子个数加1
    else{//如果之前没有添加
        pr[t]=n;//添加到表尾
        f[t++]=1;//设素因子个数初始化为1
    }
}
void find(__int64 n)
{
    __int64 p;
    if(n<2)return;
    if(prime(n))add(n);//如果n为素数，添加到素因子表中
    else{//如果n不为素数，继续拆分n
        p=pollard_rho(n);
        find(p);//尝试拆分n的一个约数p
        find(n/p);//尝试拆分n的一个约数n/p
    }
}
void findx(__int64 i,__int64 x,__int64 q)
{
    if(i==t)return;
    if(x>ans&&x<=q)ans=x;//如果约数x比ans更接近q，则更新ans
    findx(i+1,x,q);
    x*=f[i];
    if(x>ans&&x<=q)ans=x;//如果约数x比ans更接近q，则更新ans
    findx(i+1,x,q);
}
int main()
{
    __int64 i,j,k,m,n,q;
    srand((unsigned)time(NULL));
    while(~scanf("%I64d%I64d",&m,&n)){
        t=0;n/=m;
        find(n);//分解质因数
        for(i=0;i<t;i++){//将相同素因子作为一个整体处理
            for(j=1,k=pr[i];j<f[i];j++)k*=pr[i];
            f[i]=k;
        }
        q=(__int64)sqrt(n*1.0);
        ans=1;
        findx(0,1,q);//枚举n的因子,找出小于q且最接近q的因子
        printf("%I64d %I64d\n",ans*m,n/ans*m);
    }
    return 0;
}