Расширенный алгоритм Евклида

Алгоритм Евклида можно расширить так, что он не только даст НОД(a,b)=d, но и найдет целые числа x и y, такие что ax + by = d.

Псевдокод.

НА ВХОДЕ: два неотрицательных числа a и b: a>=b
НА ВЫХОДЕ: d=НОД(a,b) и целые x,y: ax + by = d.

1. Если b=0 положить d:=a, x:=1, y:=0 и возвратить (d,x,y)
2. Положить x2:=1, x1:=0, y2:=0, y1:=1
3. Пока b>0
    3.1 q:=[a/b], r:=a-qb, x:=x2-qx1, y:=y2-qy1
    3.2 a:=b, b:=r, x2:=x1, x1:=x, y2:=y1, y1:=y
4. Положить d:=a, x:=x2, y:=y2 и возвратить (d,x,y)

Исходник на Си.

/* Author:  Pate Williams (c) 1997 */

#include <stdio.h>
#define DEBUG

// расширенный алгоритм Евклида
void extended_euclid(long a, long b, long *x, long *y, long *d)
/* calculates a * *x + b * *y = gcd(a, b) = *d */
{
  long q, r, x1, x2, y1, y2;
  if (b == 0) {
    *d = a, *x = 1, *y = 0;
    return;
  }

  x2 = 1, x1 = 0, y2 = 0, y1 = 1;

  #ifdef DEBUG
  printf("------------------------------");
  printf("-------------------\n");
  printf("q    r    x    y    a    b    ");
  printf("x2   x1   y2   y1\n");
  printf("------------------------------");
  printf("-------------------\n");
  #endif

  while (b > 0) {

    q = a / b, r = a - q * b;

    *x = x2 - q * x1, *y = y2 - q * y1;

    a = b, b = r;

    x2 = x1, x1 = *x, y2 = y1, y1 = *y;

    #ifdef DEBUG
    printf("%4ld %4ld %4ld %4ld ", q, r, *x, *y);
    printf("%4ld %4ld %4ld %4ld ", a, b, x2, x1);
    printf("%4ld %4ld\n", y2, y1);
    #endif

  }

  *d = a, *x = x2, *y = y2;

  #ifdef DEBUG
  printf("-------------------\n");
  #endif
}

// главная фукнция
int main(void)
{
  long a = 4864, b = 3458, d, x, y;
  extended_euclid(a, b, &x, &y, &d);
  printf("x = %ld y = %ld d = %ld\n", x, y, d);
  return 0;
}

Алгоритм работает за O(log²n) операций.