扩展欧几里得算法

其实我最近在想要不要把写堆排序的计划取消掉……毕竟最近在忙GDKOI，GDKOI出堆排的可能性又实在不大…………
好吧，扯远了。
N周后DWJED大神要将扩展欧几里得算法，说实话，其实我之前也学过这个了，不过由于DWJED大神要讲我还是决定复习一下这个，毕竟NOIP之前复习（详见我NOIP2015的总结）的
4000
时候发现快把这事忘光了……
其实扩展欧几里得算法（下简称扩欧）是在辗转相除法求最大公约数的基础上求形如ax+by=c（a,b,c∈Z）的二元一次方程（注：不是方程组）中变量x,y（x∈Z，y∈Z）的值。当然求出方程的整数解的关键是方程本身要有整数解（众人：废话……）。问题是怎么确定有没有整数解呢？
引入一个叫贝祖等式的东西：
ax+by=c有整数解，当且仅当gcd(a,b)∣c（a,b,c∈Z）
证明的话我就不给了，毕竟我也没怎么看懂证明……
然后的话我们来推导一下怎么求解。
设axA+ byA = gcd(a,b),bxB + (a mod b)yB = gcd(b, a mod b)
由gcd(a,b)= gcd(b,a mod b)得:
a*xA + b*yA = b*xB + (a mod b)*yB
= b*xB + (a - a div b*b)*yB
= b*xB + ayB - a div b*b*yB
= a*yB + b(xB-a div b*yB)
于是我们得到了扩展欧几里得算法的核心所在：
设ax+ by=gcd(a,b),bx0 + (a mod b)y0=gcd(b, a mod b)
则x=y0，y=x0-adiv b*y0
所以我们选择递归求解。当然，需要借用一下gcd求最大公约数的方法的架子。
具体实现代码如下：
 

procedure exgcd(a,b:longint;var x,y:longint);
var t:longint;
begin
if b=0 then
begin
x:=1;
y:=0;
exit;
end
else
begin
exgcd(b,a mod b,x,y);
t:=x;
x:=y;
y:=t-(a div b)*y;
end;
end;

几个小小的推广：
1.     同余方程
这里所讲的同余方程是指形如ax≡b(mod n)的方程
其实根据同余方程的定义我们可以得到这一方程等价于ax+ny=b
用上面所讲的方法求出x即可，y没什么（luan）用。（虽说你在求解的同时还是要把y存起来……）
2.     求方程通解
这个算法有一点很明显，它只能求一组解，不过在求出这组解以后我们就可以求出别的解。
设某一组解为x0,y0，则
x:=x0+t*b/gcd(a,b)
y:=y0-t*a/gcd(a,b)
其中t∈Z
3.     最小正整数解
x0=x mod (b/gcd(a,b))
y0=y mod (a/gcd(a,b))
好的，然后我们来看一道NOI2002的题目，用上exgcd就简单很多。
NOI2002荒岛野人
岛上有排列成环行的M个山洞。这些山洞顺时针编号为1,2,…,M。岛上住着N个野人，一开始依次住在山洞C1,C2,…,CN中，以后每年，第i个野人会沿顺时针向前走Pi个洞住下来。每个野人i有一个寿命值Li，即生存的年数。奇怪的是，虽然野人有很多，但没有任何两个野人在有生之年处在同一个山洞中，使得小岛一直保持和平与宁静，这让科学家们很是惊奇。他们想知道，至少有多少个山洞，才能维持岛上的和平。
原题是有图的，但是我还是不放出来了。
其实我们只要从野人开始时居住的山洞的最大编号max(c[i])开始枚举m到10^6（m的最大值）然后我们要使对于任意的i,j∈[1,n]，都有同余方程ci+x*pi≡cj+x*pj(modm)无解或最小正整数解大于两个野人的寿命。
以上方程可以化作(pi-pj)x≡cj-ci(modm)然后按上面推广1的通式进行求解。
好吧，上代码：

var n,i,j,x,y,maxn,m:longint;
c,p,l:array[1..15]of longint;
function exgcd(a,b:longint;var x,y:longint):longint;
var t:longint;
begin
if b=0 then
begin
exgcd:=a;
x:=1;
y:=0;
end
else
begin
exgcd:=exgcd(b,a mod b,x,y);
t:=x;
x:=y;
y:=t-(a div b)*y;
end;
end;
function max(a,b:longint):longint;
begin
if a>b then exit(a) else exit(b);
end;
function min(a,b:longint):longint;
begin
if a>b then exit(b) else exit(a);
end;
function solution(a,b,n:longint):boolean;
var d,x,y,e,minx:longint;
begin
d:=exgcd(a,n,x,y);
if b mod d<>0 then exit(false);
minx:=x*b div d;
minx:=minx mod(n div d);
if minx<0 then inc(minx,abs(n div d));
if (minx<=l[i])and(minx<=l[j]) then exit(true);
exit(false);
end;
function judge(m:longint):boolean;
begin
for i:=1 to n-1 do
for j:=i+1 to n do
if solution(p[i]-p[j],c[j]-c[i],m) then exit(false);
exit(true);
end;
begin
readln(n);
maxn:=0;
for i:=1 to n do
begin
readln(c[i],p[i],l[i]);
maxn:=max(c[i],maxn);
end;
for m:=maxn to 1000000 do
if judge(m) then
begin
writeln(m);
halt;
end;
end.