扩展欧几里得算法

欧几里得算法中，计算 x, y 的最大公约数的方法是辗转相除，例如：
gcd (26, 15)
26 % 15 = 1 ... 11
15 % 11 = 1 ... 4
11 % 4 = 2 ... 3
4 % 3 = 1 ... 1
3 % 1 = 3 ... 0

可知，gcd (26, 15) = 1
如果 gcd(x, y) = r，那么有 ax + by = r，可以看出，上面的步骤实际上是可以直接得出 a, b 的：

26 % 15 = 1 ... 11 => 11 = 26 - 15 1 1 -1
15 % 11 = 1 ... 4 => 4 = 15 - 11 = 15 - (26 - 15) = -26 + 2*15 1 -1 2
11 % 4 = 2 ... 3 => 3 = 11 - 4*2 = (26 - 15) - (-26 + 15) * 2 = 3*26 - 5*15 2 3 -5
4 % 3 = 1 ... 1 => 1 = 4 - 3 = (-26 + 2*15) - (3*26 - 5*15) = -4*26 + 7*15 1 -4 7
3 % 1 = 3 ... 0

在每一轮，我们都可以得到一个模的表达式为：ri = aix + biy
如果不考虑第一轮和第二轮，那么ai 和 bi 可以表示为(qi 为每一轮得到的商)：
ai = ai-2 - qi * ai-1
bi = bi-2 - qi * bi-1

证明如下：
输入：x,y,则有如下：
x/y=q1…..r1 =>r1=x-q1*y
y/r1=q2…r2 =>r2=y-q2*r1=y-q2*(x-q1*y)=-q2*x+(1+q2*q1)*y
r1/r2=q3…r3 =>r3=r1-q3*r2=(x-q1*y)-q3*(-q2*x+(1+q2*q1)*y)
=(1+q2*q3)x+(-q1-q3*(1+q2*q3))*y
则可以看出有：a3=1+q3*q2=a1-q3*a2
B3=b1-q3*b2
由此可以推测出：
ai = ai-2 - qi * ai-1
bi = bi-2 - qi * bi-1
但是a1,b1,a2,b2比较特别：
a1=1=a-1 – q1*a0
b1=-q1=b-1 – q1*b0
由此我们可以知道a-1=1,a0=0,b-1=0,b0=1即可满足。所以这个是初试条件。
现在来考虑第一轮和第二轮，按照上面的公式，可以认为

a-1 = 1, b-1 = 0
a0 = 0, b0 = 1

有了这两对预设值，上面的两个公式就成立了

求逆元，由上面可以看出，gcd(x, y) = 1 的时候
如果 ax + by = 1，那么 ax = -by + 1 => ax = 1 (mod b)
这时候，a 即是 x 的逆元

#include <stdio.h>

  int gcd (int x, int y, int *a1, int *a2, int *b1, int *b2)
  {
    int q, r, a, b;

    q = x / y;
    r = x % y;
    a = *a2 - q*(*a1);
    b = *b2 - q*(*b1);

    if (0 == r)
      {
        return y;
      }

    *a2 = *a1;
    *b2 = *b1;

    *a1 = a;
    *b1 = b;

    return gcd (y, r, a1, a2, b1, b2);
  }

  int main ()
  {
   int x, y, a1, a2, b1, b2, r;

    scanf("%d %d", &x, &y);

    a1 = 0;
    a2 = 1;
    b1 = 1;
    b2 = 0;

    r = gcd (x, y, &a1, &a2, &b1, &b2);

    printf ("%d = (%d) * (%d) + (%d) * (%d)
",
            r, a1, x, b1, y);

    return 0;
  }


输入样例：
7 160
输出：
//7 的 mod 160 ==1 的逆元为 23