欧几里得算法和扩展欧几里得算法
欧几里得算法看上去高大上，其实质是辗转相除法，用于计算两个整数的a,b的最大公约数。
基本算法: 假设 $$ a=qb+r $$ 其中a,b,q,r都是整数, 则gcd(a,b) = gcd(b,r); 证明可以参考链接欧几里得与扩展欧几里得
扩展欧几里得描述的是: $$ap+b*q=c$$(a,b,c已知, p,q未知), 若存在解则c mod gcd(a,b) = 0, 并且p,q的其中一个解可以通过扩展的欧几里得算法求得。
首先求$$ap + bq = gcd(a,b)$$的解p,q
不妨假设 a > b, a, b不同时为0
1) 当b = 0, 则p=1, q = 0原等式化为 a = a
2) 当a,b都不为0的时候, $$ap1 + bq1 = gcd(a,b)$$
$$bp2 + (a mod b)q2 = gcd(b, a mod b)$$
推出:
$$ap1 + bq1 = bp2 + (a mod b)q2$$
=> $$p1 = q2$$
$$q1 = p2 - (a/b)*q2$$
即递归的边界在b=0时返回p=1, q = 0, 否则依次按上式递归
代码描述:1
2
3
4
5
6
7
8
9
10
11
12
13
14__int64 exgcd(__int64 a, __int64 b, __int64 &x, __int64 &y)
{
if (b==0)
{
x = 1;
y = 0;
return a;
}
__int64 x2,y2;
__int64 t = exgcd(b,a%b,x2,y2);
x = y2;
y = x2 - (a/b)*y2;
return t;
}
POJ1061
题目描述:
两只青蛙在网上相识了，它们聊得很开心，于是觉得很有必要见一面。它们很高兴地发现它们住在同一条纬度线上，于是它们约定各自朝西跳，直到碰面为止。可是它们出发之前忘记了一件很重要的事情，既没有问清楚对方的特征，也没有约定见面的具体位置。不过青蛙们都是很乐观的，它们觉得只要一直朝着某个方向跳下去，总能碰到对方的。但是除非这两只青蛙在同一时间跳到同一点上，不然是永远都不可能碰面的。为了帮助这两只乐观的青蛙，你被要求写一个程序来判断这两只青蛙是否能够碰面，会在什么时候碰面。
我们把这两只青蛙分别叫做青蛙A和青蛙B，并且规定纬度线上东经0度处为原点，由东往西为正方向，单位长度1米，这样我们就得到了一条首尾相接的数轴。设青蛙A的出发点坐标是x，青蛙B的出发点坐标是y。青蛙A一次能跳m米，青蛙B一次能跳n米，两只青蛙跳一次所花费的时间相同。纬度线总长L米。现在要你求出它们跳了几次以后才会碰面。
Input
输入只包括一行5个整数x，y，m，n，L，其中x≠y < 2000000000，0 < m、n < 2000000000，0 < L < 2100000000。
Output
输出碰面所需要的跳跃次数，如果永远不可能碰面则输出一行”Impossible”
假设进过次数t, k圈两只青蛙相遇:
$$x+mt - (y+nt) = kl$$
=> $$(n-m)t + lk = x - y$$
和扩展欧几里得公式对应起来
$$a = n-m, b = l, c = x - y,
p = t, q = k$$
那么调用扩展欧几里得公式后得到的解 $$(gcd(a,b)) = d = exgcd(a,b,p,q)$$
p,q是$$ap + bq = gcd(a,b)$$的解, 那么你需要实际的此时应该是$$xx = pc/gcd(a,b) = t(x-y)/d$$ 但是还有一个问题就是你的取值范围问题。
如下证明:
定理一：如果d = gcd(a, b)，则必能找到正的或负的整数k和l，使d = ax+ by。
定理二：若gcd(a, b) = 1，则方程ax ≡ c (mod b)在[0, b-1]上有唯一解。
定理三：若gcd(a, b) = d，则方程ax ≡ c (mod b)在[0, b/d - 1]上有唯一解。
证明：上述同余方程等价于ax + by = c，如果有解，两边同除以d，就有a/d x + b/d y = c/d，即a/d x ≡ c/d (mod b/d)，显然gcd(a/d, b/d) = 1，所以由定理二知道x在[0, b/d - 1]上有唯一解。所以ax + by = c的x在[0, b/d - 1]上有唯一解，即ax ≡ c (mod b)在[0, b/d - 1]上有唯一解。
如果得到ax ≡ c (mod b)的某一特解X，那么令r = b/gcd(a, b)，可知x在[0, r-1]上有唯一解，所以用x = (X % r + r) % r就可以求出最小非负整数解x了！（X % r可能是负值，此时保持在[-(r-1), 0]内，正值则保持在[0, r-1]内。加上r就保持在[1, 2r - 1]内，所以再模一下r就在[0, r-1]内了）。
r = b/gcd(a,b), 所以实际的次数是 (xx%r + r)%r
说的这么多该上代码了: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#include <iostream>
using namespace std;
__int64 exgcd(__int64 a, __int64 b, __int64 &x, __int64 &y)
{
if (b==0)
{
x = 1;
y = 0;
return a;
}
__int64 x2,y2;
__int64 t = exgcd(b,a%b,x2,y2);
x = y2;
y = x2 - (a/b)*y2;
return t;
}
int main()
{
__int64 x,y,m,n,l;
while (scanf("%I64d%I64d%I64d%I64d%I64d",&x,&y,&m,&n,&l)==5)
{
__int64 a,b,c,p,q,d,r,t;
a = n-m, b = l, c = x-y;
d = exgcd(a,b,p,q);
if (c%d!=0)
printf("Impossible\n");
else{
r = b/d;
t = p*(c/d);
printf("%I64d\n",(t%r+r)%r);
}
}
return 0;
}
POJ2262
问题描述: 每一个大于4的偶数都可以被写成两个奇素数相加
比如
8 = 3 + 5
20 = 3 + 17
42 = 5 + 37
输入一个n>=6, n < 100000寻找间隔最大的两个奇素数相加表示, 如果存在按照上面的格式输出, 否则输出Goldbach’s conjecture is wrong.
Sample Input
8
20
42
0
Sample Output
8 = 3 + 5
20 = 3 + 17
42 = 5 + 37
此题很简单, 只需要从最小的素数开始遍历就行, 如果差数也是素数就找到了。
代码如下: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#include <iostream>
using namespace std;
bool isprime(int x)
{
for (int i=2;i*i<=x;i++)
{
if (x%i==0)
return false;
}
return true;
}
int main()
{
int n;
while (scanf("%d",&n)==1,n)
{
bool flag = false;
for (int i=3;i<=n/2;i=i+2)
{
if (isprime(i) && isprime(n-i))
{
printf("%d = %d + %d\n",n,i,n-i);
flag = true;
break;
}
}
if (!flag)
printf("Goldbach's conjecture is wrong.\n");
}
return 0;
}
POJ1730
问题描述:
我们知道25 = 5^2, 1073741824 = 2^30, 等等，现在给你一个32位的int型数n, 你需要找到最大的指数i使得n = a^i
这个题有个很巧妙的解法:
首先我们来思考一下
1.32位的int型数是包括正数和负数的, 我们知道负数是不可能等于某个数的偶数次方的。
2.我们知道双精度的double表示是1, 11, 52, 那么我们将一个double x强制转化成__int64 a, 显然如果abs(a-x)< (1.0e-12)或者abs(a-x+1) < (1.0e-12)(这个地方你可以自己检验一下), 就可以保证x == a。
3.可以直接对n进行开32~1次方根
代码如下: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#include <iostream>
#include <math.h>
using namespace std;
#define esp 1.0e-12 // 1 11 52 double型需要12位
int main()
{
long long n;
while (scanf("%I64d",&n)==1 && n!=0)
{
bool flag = false;
if (n < 0)
{
flag = true;
n = -n;
}
for (int i=32;i >= 1;i--)
{
if (!flag || (flag && (i&1))) // 负数只能开奇数次方不能开偶数次方
{
double x = pow((double)n, 1.0/i); // 1 8 23 1 11 52
long long a = x;
if (abs(a-x) < esp || abs(a-x+1) < esp)
{
printf("%d\n",i);
break;
}
}
}
}
return 0;
}
POJ1183
问题描述:根据此题和反正切函数没有半毛钱关系
反正切等式: $$ arctan(p)+arctan(q)=arctan[(p+q)/(1-pq)] $$
我们将上述公式写成如下形式
$$ arctan(1/a)=arctan(1/b)+arctan(1/c) $$
其中a,b和c均为正整数。
我们的问题是：对于每一个给定的a(1<=a<=60000)，求b＋c的值。我们保证对于任意的a都存在整数解。如果有多个解，要求你给出b+c最小的解。
证明推导:
$$ 1/a = (1/b+1/c)/(1-1/b1/c) $$
=> $$a(b+c) = bc -1$$
利用均值不等式我们可以推出$$b+c >= 2(a+sqrt(aa+1))$$
不妨设$$b = a + m, c = a + n$$
$$a(b+c) = bc -1 => mn = aa + 1$$
因为需要求b+c的最小值也就是m,n的最小值因而很容易理解m从a开始往下搜, (aa+1)%m==0结束
最后$$b+c = 2*a + m + n$$
代码如下:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19#include <iostream>
using namespace std;
int main()
{
__int64 a,b,c;
while (scanf("%I64d",&a)==1)
{
__int64 m,n;
for (m=a;m>=1;m--)
{
if ((a*a+1) %m ==0)
break;
}
n = (a*a+1)/m;
printf("%I64d\n",2*a+m+n);
}
return 0;
}