扩展欧几里得算法

其实我最近在想要不要把写堆排序的计划取消掉……毕竟最近在忙GDKOI，GDKOI出堆排的可能性又实在不大…………

好吧，扯远了。

N周后DWJED大神要将扩展欧几里得算法，说实话，其实我之前也学过这个了，不过由于DWJED大神要讲我还是决定复习一下这个，毕竟NOIP之前复习（详见我NOIP2015的总结）的时候发现快把这事忘光了……

其实扩展欧几里得算法（下简称扩欧）是在辗转相除法求最大公约数的基础上求形如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.