POJ 2429 GCD & LCM Inverse 大素数分解

题目链接:http://poj.org/problem?id=2429

题意

给出两个数的最大公因数和最小公倍数(数据范围2^64),求这两个数(存在多组数时输出和最小的一组数)

题解

我们很容易得到以下方程: (a/gcd)*(b/gcd)=(lcm/gcd)。因为(a/gcd)和(b/gcd)一定是互质的(如果不互质,gcd就要改变),这样我们就可以看成是将(lcm/gcd)分解成互质的两个数。
使用Pollard-Rho算法算出大整数的素因子表,为了保证分解成的两个数是互质的,将素质因子表中相同的数进行相乘,可以证明这样得到的数组内的元素之间还是互质的。
最后只需要对表中的元素分成两组就行了,用一个简单的DFS就可以搞定。

代码(处于Runtime error )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
#include <iostream>
#include<time.h>
#include<stdlib.h>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <string>
#include <vector>
#include <cmath>
#include <set>
#include <queue>
#include <utility>

using namespace std;
#define INF 0x3f3f3f3f
typedef long long ll;

//Miller_Rabin 算法进行素数测试
//速度快,而且可以判断 <2^63的数

const int S=20; //随机算法判定次数

//计算 (a*b)%c 加法快速幂
ll mul_mod(ll a,ll b,ll 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;
}

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

//以a为基,n-1=x*2^t a^(n-1)=1(mod n) 验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(ll a,ll n,ll x,ll t)
{
ll ret=pow_mod(a,x,n);
ll last=ret;
for(int i=1;i<=t;i++)
{
ret=mul_mod(ret,ret,n);
if(ret==1 && last!=1 && last!=n-1) return true;
last=ret;
}
if(ret!=1) return true;
return false;
}

// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数,但概率极小)
//合数返回false;
bool Miller_Rabin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if((n&1)==0) return false; //偶数
ll x=n-1,t=0;
while(!(x&1))
{
x>>=1;
t++;
}
for(int i=0;i<S;i++)
{
ll a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;
}
return true;
}


// pollard_rho 算法进行质因数分解
ll factor[10000]; //分解结果
int tol; //分解个数

ll gcd(ll a,ll b)
{
if(a==0) return 1;
if(a<0) return gcd(-a,b);
return b==0?a:gcd(b,a%b);
}

ll Pollard_rho(ll x,ll c)
{
ll i=1,k=2;
ll x0=rand()%x;
ll y=x0;
while(1)
{
i++;
x0=(mul_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;}
}
}

//对n进行素因子分解
void findfac(ll n)
{
if(Miller_Rabin(n))
{
factor[tol++]=n;
return;
}
ll p=n;
while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
findfac(p);
findfac(n/p);
}

//DFS 求数
ll min_sum; //和的最小值设置为全局变量
ll ansa,ansb;
ll factor_new[10000],num1; //将重复素因子相乘以后的表

void dfs(ll tmpa,ll tmpb,ll pos,const ll len)
{
if(pos==len)
{
if(tmpa+tmpb<=min_sum)
ansa=tmpa,ansb=tmpb;
return;
}
else
{
dfs(tmpa*factor_new[pos],tmpb,pos+1,len);
dfs(tmpa,tmpb*factor_new[pos],pos+1,len);
}
}

void solve()
{
ll a,b;
while(scanf("%lld %lld",&a,&b)!=EOF)
{
num1=0;
b/=a;
findfac(b);
sort(factor,factor+tol);
factor_new[0]=factor[0];
for(int i=1;i<tol;i++) //除去重复素因子
{
if(factor[i]==factor[i-1])
factor_new[num1]*=factor[i];
else
factor_new[++num1]=factor[i];
}
min_sum=factor_new[0]+b/factor_new[0];
dfs(1,1,0,num1+1); //dfs 枚举结果
if(ansa>ansb) swap(ansa,ansb);
printf("%lld %lld\n",ansa*a,ansb*a);
}
}

int main()
{
freopen("input.txt","r",stdin);

//test();
solve();
return 0;
}