poj 2429 GCD & LCM Inverse

缘起

日常浪费生命 poj 2429 GCD & LCM Inverse

分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
两个正整数a,b, 已知gcd(a,b)和lcm(a,b) 求a和b

【输入】
多样例. 每个样例两个正整数，分别是gcd和lcm,并且gcd和lcm都<=2^63

【输出】
按照a<=b的顺序输出a和b,如果有多组解答, 输出a+b最小的那一组

【样例输入】
3 60

【样例输出】
12 15

【限制】
2s

令(a,b) = gcd, 则 a = gcd*a1, b = gcd*b1, lcm = gcd*a1*b1

所以a和b，知道了任何一个，另一个就确定了. 而lcm/gcd 得到的是a1*b1

所以我们需要将lcm/gcd的结果分解为2个互质的数的乘积. 一旦知道了a1,则b1也就知道了. 则乘以gcd,两个正整数也就知道了. 所以此题的解决之道就是将lcm/gcd分解为两个互质正整数的乘积.

而要求a+b=gcd*(a1+b1) 最小的那一组. 所以根据数学理论

1
两数乘积固定, 两数越接近则和越小.

所以要将 lcm/gcd 进行分解素因数. 然后因为要求出互质的两个数相乘, 所以其实是每种质因子只能二择一. 什么意思? 例如 108=2^2*3^3, 则2为4, 3为27. 则只能一个为4，一个为27. 不能两个都既有2因子又有3因子, 因为那样的话, 就不互质了。

所以可以dfs. 怎么讲，例如 lcm/gcd得到的分解素因数为

4*27*5

则考虑构成a1, 这三个数只能考虑取与不取. 而dfs的时候要考虑剪枝. 因为题目的范围到达了long long, 素因子很可能到达20位，如果不剪枝的话，一定T. 所以dfs的顺序也很重要——要尽快的提高标准，让不符合的搜索顺序被尽快pass掉. 下面的代码不懂的，可以参考【1】和【2】

然后辛辛苦苦写了200+行的代码，结果被T掉了

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
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
//#include "stdafx.h"

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <algorithm>
#include <map>
#include <math.h>
using namespace std;
//#define LOCAL
typedef long long LL;

int num;  // num是素因子的个数
map<LL, int> primefactors; // 分解素因数的结果
LL primes[100],lim,ans_a1,times[100]; // 用于存储素因子,times[i]表示 primes[i,...,num-1]的乘积
const int test = 12, A = 240;
typedef map<LL,int>::iterator mit;

LL modu_mul(LL a, LL b, LL mod) // a*b 模掉mod的余数
{
	a %= mod, b%=mod;
	int ans = 0;
	while(b)
	{
		if (b&1)
		{
			ans = (ans+a)%mod;
		}
		if (b>1)
		{
			a = (a+a)%mod;
		}
		b>>=1;
	}
	return ans;
}

LL modu_pow(LL a, LL b, LL mod) // a^b模掉mod的余数
{
	LL ans = 1;
	a%=mod;
	while(b)
	{
		if (b&1)
		{
			ans = modu_mul(ans, a, mod);
		}
		if (b>1)
		{
			a = modu_mul(a,a, mod);
		}
		b>>=1;
	}
	return ans;
}

LL gcd(LL a, LL b) // 求a和b的最大公约数
{
	return b?gcd(b, a%b):a;
}

bool mr(LL p) // Miller-Rabin 素性测试
{
	if (p<2)
	{
		return false;
	}
	if (!(p&1))
	{
		return p==2;
	}
	LL s = 0, t = p-1;
	while(!(t&1))
	{
		s++;
		t>>=1;
	}
	for (int i = 0; i<test; i++)
	{
		LL a = rand()%(p-1)+1;
		LL x = modu_pow(a,t,p);
		for (int j = 0; j<s; j++)
		{
			LL y = modu_mul(x,x,p);
			if (y==1 && x!=1 && x!=p-1)
			{
				return false;
			}
			x = y;
		}
		if (x!=1)
		{
			return false;
		}
	}
	return true;
}

LL pollard_rho(LL n, LL a) // 返回n的一个非平凡因子
{
	int i = 1, k=2,x = rand()%n, y = x;
	while(1)
	{
		i++;
		x = (modu_mul(x,x,n)+a)%n;
		LL d = gcd(y-x+n, n); // 这是一种不用取abs的高效方法
		if (1<d&& d<n)
		{
			return d;
		}
		if (x==y)
		{
			return n;
		}
		if (i==k)
		{
			y = x;
			k<<=1;
		}
	}
}

void factorize(LL n, LL a)
{
	if (mr(n))
	{
		if (primefactors[n])
		{
			primefactors[n]*=n;
		}
		else
		{
			primefactors[n] = n;
		}
		return;
	}
	LL d = n;
	while(d==n)
	{
		d = pollard_rho(n, a--);
	}
	factorize(d,a);
	factorize(n/d, a);
}

bool cmp(LL &a, LL &b)
{
	return a>b;
}

void dfs(int i, LL a1) // i表示当前决定的是primes[i]的取与不取, a1是当前已经达到的乘积
{
	if (i==num) // 说明已经达到了一组解
	{
		if (ans_a1<a1)
		{
			ans_a1 = a1;
		}
		return;
	}
	if (a1>lim) // 已经逾界了,剪枝!
	{
		return;
	}
	if (a1*times[i]<=ans_a1) // 如果就算后面所有的都选了也没办法达到ans_a1当前的高度的话,剪枝!
	{
		return;
	}
	dfs(i+1, a1*primes[i]); // 先考虑选primes[i], 再考虑不选primes[i],而不是反过来考虑,这也是为了剪枝的效率! 因为我们想先快速得到一个近似解用于高效剪枝
	dfs(i+1, a1);
}

int main()
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
	//freopen("d:\\my.out", "w", stdout);
#endif
	LL x,y;
	while(~scanf("%lld%lld", &x, &y))
	{
		if (x==y) // 如果x和y相等
		{
			printf("%lld %lld\n", x,y);
			continue;
		}
		ans_a1 = 1;
		primefactors.clear();
		num = 0;
		factorize(y/x, A); // 对y/x 进行分解素因数,结果保存在 primefactors 中.
		lim = (int)(sqrt(1.0*y/x)+0.5); // 防止精度丢失
		for (mit i = primefactors.begin(); i!=primefactors.end(); i++)
		{
			if (i->second<=lim) // i->second都是素因子的幂次, 只收集不超过平方根的这种幂次即可
			{
				primes[num++] = i->second;
			}
		} // 收集, 例如 20 = 2^2*5, 则num=2(表示2种素因子), primes[0] = 4, primes[1] = 5
		sort(primes, primes+num, cmp); // 从大到小排序——注意, 这是为了dfs的效率——为了能快速得到一个近似解从而实现大量剪枝
		primes[num++] = 1; //加个1
		times[num-1] = 1;
		for (int i = num-2; ~i; i--)
		{
			times[i] = times[i+1]*primes[i];
		} // 预处理得到times, 这个times数组是为了剪枝用
		dfs(0,1LL); // 求出a=x*a1,b=x*b1 的 a1来, 它是不超过 lim=sqrt(y/x)的y/x的最大因子
		printf("%lld %lld\n", x*ans_a1, y/ans_a1);
	}
	return 0;
}

结果是

Status Time Limit Exceeded
Length 3274
Lang C++
Submitted 2019-08-31 18:15:36
Shared
RemoteRunId 20821941

后来求得一组测试数据

1
2
3
3 67553994410557440
答案是
15 13510798882111488

上述测试数据我洗了个澡回来都没有跑出结果来. 原因很显然———分解超大整数的素因数算法超时. 那是pollard-rho 算法效率不够么? 不是（这里的pollard-rho的板子是经过【1】和【2】考验的, 效率是有保证的）. 而是 factorize 效率太低了. 为什么? 因为factorize方法是递归，而且找到一个素因子不把n的该素因子都除掉. 所以factorize的递归写法只是好看，简洁，易懂, 但是效率不够. 正确的做法显然是用pollard-rho算法找到一个素因子之后就立马将这个素因子除干净. 然后再找下一个. 所以给出如下ac代码

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
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
//#include "stdafx.h"

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <algorithm>
#include <math.h>
using namespace std;
//#define LOCAL
typedef long long LL;

int num;  // num是素因子的个数
LL primes[100],lim,ans_a1,times[100]; // 用于存储素因子,times[i]表示 primes[i,...,num-1]的乘积
const int test = 12, A = 240;

LL modu_mul(LL a, LL b, LL mod) // a*b 模掉mod的余数
{
	a %= mod, b%=mod;
	LL ans = 0;
	while(b)
	{
		if (b&1)
		{
			ans = (ans+a)%mod;
		}
		if (b>1)
		{
			a = (a<<1)%mod;
		}
		b>>=1;
	}
	return ans;
}

LL modu_pow(LL a, LL b, LL mod) // a^b模掉mod的余数
{
	LL ans = 1;
	a%=mod;
	while(b)
	{
		if (b&1)
		{
			ans = modu_mul(ans, a, mod);
		}
		if (b>1)
		{
			a = modu_mul(a,a, mod);
		}
		b>>=1;
	}
	return ans;
}

LL gcd(LL a, LL b) // 求a和b的最大公约数
{
	return b?gcd(b, a%b):a;
}

bool mr(LL p) // Miller-Rabin 素性测试
{
	if (p<2)
	{
		return false;
	}
	if (!(p&1))
	{
		return p==2;
	}
	LL s = 0, t = p-1;
	while(!(t&1))
	{
		s++;
		t>>=1;
	}
	for (int i = 0; i<test; i++)
	{
		LL a = rand()%(p-1)+1;
		LL x = modu_pow(a,t,p);
		for (int j = 0; j<s; j++)
		{
			LL y = modu_mul(x,x,p);
			if (y==1 && x!=1 && x!=p-1)
			{
				return false;
			}
			x = y;
		}
		if (x!=1)
		{
			return false;
		}
	}
	return true;
}

LL pollard_rho(LL n, LL a) // 返回n的一个非平凡因子
{
	LL i = 1, k=2,x = rand()%n, y = x;
	while(1)
	{
		i++;
		x = (modu_mul(x,x,n)+a)%n;
		LL d = gcd(y-x+n, n); // 这是一种不用取abs的高效方法
		if (1<d&& d<n)
		{
			return d;
		}
		if (x==y)
		{
			return n;
		}
		if (i==k)
		{
			y = x;
			k<<=1;
		}
	}
}

LL getprime(LL n, LL a) // 获取n的一个素因子或者1
{
	if (n==1)
	{
		return 1;
	}
	if (mr(n))
	{
		return n;
	} // 下面n既不是1也不是素数, 即是一个合数,所以得到的d一定是非平凡因子
	LL d = n;
	while(d==n)
	{
		d = pollard_rho(n, a--);
	}
	return getprime(d,a);
}

bool cmp(LL &a, LL &b)
{
	return a>b;
}

void dfs(int i, LL a1) // i表示当前决定的是primes[i]的取与不取, a1是当前已经达到的乘积
{
	if (i==num) // 说明已经达到了一组解
	{
		if (ans_a1<a1)
		{
			ans_a1 = a1;
		}
		return;
	}
	if (a1>lim) // 已经逾界了,剪枝!
	{
		return;
	}
	if (a1*times[i]<=ans_a1) // 如果就算后面所有的都选了也没办法达到ans_a1当前的高度的话,剪枝!
	{
		return;
	}
	dfs(i+1, a1*primes[i]); // 先考虑选primes[i], 再考虑不选primes[i],而不是反过来考虑,这也是为了剪枝的效率! 因为我们想先快速得到一个近似解用于高效剪枝
	dfs(i+1, a1);
}

int main()
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
	freopen("d:\\my.out", "w", stdout);
#endif
	LL x,y;
	while(~scanf("%lld%lld", &x, &y))
	{
		if (x==y) // 如果x和y相等
		{
			printf("%lld %lld\n", x,y);
			continue;
		}
		ans_a1 = 1;
		num = 0;
		LL t = y/x;
		lim = (LL)(sqrt(1.0*y/x)+0.5); // 防止精度丢失
		while(t>1)
		{
			LL p = getprime(t,A);
			bool flag = true;
			primes[num] = 1;
			while(!(t%p))
			{
				primes[num]*=p;
				if (primes[num]>lim)
				{
					flag = false;
				}
				t/=p;
			}
			if (flag)
			{
				num++;
			}
		}
		sort(primes, primes+num, cmp); // 从大到小排序——注意, 这是为了dfs的效率——为了能快速得到一个近似解从而实现大量剪枝
		primes[num++] = 1; //加个1
		times[num-1] = 1;
		for (int i = num-2; ~i; i--)
		{
			times[i] = times[i+1]*primes[i];
		} // 预处理得到times, 这个times数组是为了剪枝用
		dfs(0,1LL); // 求出a=x*a1,b=x*b1 的 a1来, 它是不超过 lim=sqrt(y/x)的y/x的最大因子
		printf("%lld %lld\n", x*ans_a1, y/ans_a1);
	}
	return 0;
}

ac情况 (水过，水过~ 竟然跑了719ms这么久!)

20822074 yfsyfs 2429 Accepted 140K 719MS C++ 3078B 2019-08-31 19:29:46

参考

【1】https://yfsyfs.github.io/2019/08/17/poj-1811-Prime-Test-Pollard-rho%E7%AE%97%E6%B3%95-%E5%88%86%E8%A7%A3%E7%B4%A0%E5%9B%A0%E6%95%B0-%E5%BF%AB%E9%80%9F%E5%B9%82/

【2】https://yfsyfs.github.io/2019/08/18/hdu-4344-Mark-the-Rope-pollard-rho/