挑战程序设计解题报告 2.6.1辗转相除法

挑战程序设计解题报告 2.6.1辗转相除法

1.POJ 2429

        这道理的意思就是告诉你GCD和LCM让你求原数，但是有个问题，很明显有很多情况，所以题目要求输出满足a+b最小的情况，怎样才能最小呢，已知GCD(a, b) == n, lcm(a, b) == m,那么a * b == n / m（由GCD和LCM的定义可知），问题变成了找到两个数a,b使得a * b == n / m、且a+b最小，即f(x)
= a * b = a  +  (n / m) / a，这个函数是对勾函数，最小值是a等于分子的平方根，如图：

 


        但是盲目去找是不行的，所以我们需要将n/m分解成素因子，由于同一个素因子可能有很多，所以需要先合并在一起，然后从这些数字及其乘积的组合中去寻找最接近上图最低值的点，首先是大整数的分解需要用到Pollard_rho，然后是大素数判断Miller_Rabin，这些都是模板，假设我们现在得到数字以后，只需要使用挑战程序设计这本书的31页的穷遏搜索即可，也就是从一堆数字中选几个数字乘在一起使得它的值接近给定值。

        例如2 120这组数据，分解形成 2 2 3 5，合并形成3 4 5，搜索的时候检查3 4 5 12 15 20 60最终发现5更接近sqrt（120 / 2），所以结果就是 5 和 12 即为 10 24

代码如下：

const int S = 10;    //米勒罗宾素数测试

map<LL, LL>v;    //记录并合并多余的素因子
vector<LL>ans;    //记录最终需要进行搜索的数字
LL mult_mod(LL a, LL b, LL c)    //计算a * b % c
{
a %= c;
b %= c;
LL ret = 0;
while(b)
{
if(b & 1)
{
ret += a;
ret %= c;
}
a <<= 1;
if(a >= c)a %= c;
b >>= 1;
}
return ret;
}

LL pow_mod(LL x, LL n, LL mod)  //计算a^b%c
{
if(n == 1)
return x % mod;
x %= mod;
LL tmp = x;
LL ret = 1;
while(n)
{
if(n & 1) ret = mult_mod(ret, tmp, mod);
tmp = mult_mod(tmp, tmp, mod);
n >>= 1;
}
return ret;
}

bool Miller_Rabin(LL n)     //模板 米勒罗宾大素数测试
{
if(n == 2)
return true;
if(n < 2 || !(n & 1))
return false;
int t = 0;
const int S = 10;

LL a, x, y, u = n - 1;
while((u & 1) == 0)
{
t++;
u >>= 1;
}
for(int i = 0; i < S; i++)
{
a = rand() % (n - 1) + 1;
x = pow_mod(a, u, n);
for(int j = 0; j < t; j++)
{
y = mult_mod(x, x, n);
if(y == 1 && x != 1 && x != n - 1)
return false;
x = y;
}
if(x != 1)
return false;
}
return true;
}

LL gcd(LL a, LL b)    //gcd
{
if(a == 0)return 1;
if(a < 0) return gcd(-a, b);
while(b)
{
LL t = a % b;
a = b;
b = t;
}
return a;
}

LL Pollard_rho(LL x, LL c)    //得到其中一个因子，配合findfac一起找到所有的因子
{
LL i = 1, k = 2;
LL x0 = rand() % x;
LL y = x0;
while(1)
{
i++;
x0 = (mult_mod(x0, x0, x) + c) % x;
LL d = gcd(y - x0, x);
if(d != 1 && d != x) return d;
if(y == x0) return x;
if(i == k)
{
y = x0;
k += k;
}
}
}

void findfac(LL n)
{
if(n == 1)
return;
if(Miller_Rabin(n))
{
if(v
== 0)    //记录并合并素因子
v
= n;
else
v
*= n;
return;
}
LL p = n;
while(p >= n)
p = Pollard_rho(p, rand() % (n - 1) + 1);
findfac(p);    //利用二分的思想找其余的因子
findfac(n / p);    //
}

void bfs(LL i, LL sum, LL& key, LL p)
{
if(abs(key - p) > abs(sum - p))    //用key保存最接近p的值
key = sum;
if(i == ans.size())
return;     //递归出口
bfs(i + 1, sum, key, p);    //不选i+1这个数字
bfs(i + 1, sum * ans[i], key, p);    //选上i+1这个数字
}

int main()
{
LL N, M;
while(cin >> N >> M)
{
v.clear();
ans.clear();
LL m = M / N;
findfac(m);   //找到因子
for(map<LL, LL>::iterator it = v.begin(); it != v.end(); ++it)
{
ans.push_back(it->second);    //提取到vector中
}
LL a = 1;
bfs(0, 1, a, sqrt(m));    //找到最接近M / N的值
LL b = m / a;
if(a < b)
cout << a * N << " " << b * N << endl;
else
cout << b * N << " " << a * N << endl;
}
return 0;
}

        这道题为了解决素数检测，因式分解。花了很多时间，虽然是模板但是这两个算法的说明，我将会在另一篇文章中指出。