Продолжение про Евклида

Собственно, не рекурсивный алгоритм:

function ExtEuclid(X, Y : Int64; var A, B : Int64) : Int64;
// возвращает наибольший общий делитель X и Y, X >= Y
var
  B1, B1p, XS, YS : Int64;
  Q, T : Int64;
begin
  XS := X;
  YS := Y;
  B1p := 0;
  B1 := 1;
  while Y <> 0 do
  begin
    Q := X div Y;
    T := Y;
    Y := X - Q * Y;
    X := T;
    if Y <> 0 then
    begin
      T := B1;
      B1 := B1p - Q * B1;
      B1p := T;
    end;
  end;
  Result := X;
  B := B1;
  A := (X - B1*YS) div XS;
end;

Ну и его ассемблерный вариант:

function ExtEuclidASM(X, Y : UInt64; var A, B : Int64) : UInt64;
// возвращает наибольший общий делитель X и Y, X >= Y
asm
  push RDI; // B1Lo
  push RSI; // B1Hi
  push RBX; // Y;
  push R12; //XS;
  push R13; //YS
  push R14; // Временный
  push R15; // Временный
  mov R12, RCX;
  mov R13, RDX;
  mov RBX, RDX;
  xor R10, R10;
  xor R11, R11;
  xor RSI, RSI;
  xor RDI, RDI;
  inc DI;
@Loop:
  xor RDX, RDX;
  mov RAX, RCX;
  div RBX;
  mov RCX, RBX;
  mov RBX, RDX;
  test RDX, RDX;
  jz @FinLoop;
  mov R14, RAX
  mul RDI;
  xchg R14, RAX;
  mov R15, RDX;
  imul RSI;
  add R15, RAX;
  sub R10, R14;
  sbb R11, R15;
  xchg RDI, R10;
  xchg RSI, R11;
  jmp @Loop;
@FinLoop:
  mov R14, B;
  mov [R14], RDI;
  mov RAX, RDI;
  mul R13;
  mov R14, RAX;
  mov R15, RDX;
  mov RAX, RSI;
  imul R13;
  add R15, RAX;
  mov RAX, RCX;
  xor RDX, RDX;
  sub RAX, R14;
  sbb RDX, R15;
  idiv R12;
  mov R14, A;
  mov [R14], RAX;
  mov RAX, RCX;
  pop R15;
  pop R14;
  pop R13;
  pop R12;
  pop RBX;
  pop RSI;
  pop RDI;
end;

Бинарный вообще не стал переделывать по следующим соображениям: для моих задач нужна версия, которая работает на полном беззнаковым 64-битном целом, что требует 128-битных вычислений. Соответственно, отказ от второго if во внешнем цикле с обменом уравнений (с целью упростить условие в цикле repeat) потребует пять регистровых обменов, а в результате условие упроститься всего на 4 команды, так что особого выигрыша не будет.
Ресурсы по использованию conditional move исчерпаны: я использую все 16 регистров в теле функции, использовать же память для хранения временных значений, которые, возможно, понадобятся,  очень плохая идея с точки зрения производительности. А если вспомнить, что я использую 128-битную арифметику ‒ то просто катастрофичная.

Теперь результаты ассемблерных вариантов расширенного алгоритма Евклида (паскалевские чуть медленнее, да и сравнивать их не корректно: они используют более короткую 64-битную арифметику):

• рекурсивный ‒ 7.05 млн./сек;
• обычный ‒ 9.7 млн./сек;
• бинарный ‒ 1.47 млн./сек.
  

Видимо, что бы бинарный работал быстрее обычного, нужно что бы одна итерация бинарного выполнялась примерно в 4 раза быстрее.
Так как полный вариант деления и умножения и так уже почти полностью 128-битный у процессоров x64, то тут дополнительных издержек у обычного алгоритма нет, чего не скажешь про бинарный.
Но если потребуется выполнять те же действия с более длинными данными, например, 256-битными, то, видимо, бинарный алгоритм снова будет впереди.

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

Ситуация сложилась так, что временно три компьютера мне не доступны для теста производительности на решении СЛАУ методом Гаусса-Жордана. Поэтому пока я жду их доступности, что бы получить результаты, и тогда проведу детальный анализ.

Во время ожидания у меня возникла необходимость посчитать обратный элемент по модулю. Это можно делать либо с помощью расширенного алгоритма Евклида, либо по теореме Эйлера с помощью быстрого возведения в степень.

Сам-то я не очень силен в модулярной арифметике, поэтому принялся  искать готовое описание, годное для реализации, расширенного алгоритма Евклида. Не скажу, что бы прям сильно глубоко искал, но потратил в районе часа.
И выяснил в результате, что в русскоязычном сегменте Инета нет годного внятного описания алгоритма. В англоязычном искал существенно меньше, но и там, похоже, ситуация не лучше. То, что нашел, работает лишь частично верно.

В результате получил верно работающую лишь рекурсивную версию классического расширенного алгоритма Евклида. Догадываюсь, что можно перенести на итерации, но что-то туплю, и не могу понять, как это сделать. И самое главное: не ясно, будет ли с этого какой толк в плане производительности. Возможно, кому-то понадобится рабочий вариант, поэтому выкладываю его на Pascal:

function ExtEuclidRec(X, Y : Int64; var A, B : Int64) : Int64;
var
  OldA : Int64;
begin
  OldA := X mod Y;
  if OldA = 0 then
  begin
    Result := Y;
    A := 0;
    B := 1;
  end
  else
  begin
    Result := ExtEuclidRec(Y, OldA, A, B);
    OldA := A;
    A := B;
    B := OldA - B*(X div Y);
  end;
end;

Уточню обозначения. Расширенный алгоритм Евклида фактически решает линейное диофантово уравнение Ax + By = D.
Встречаются разные трактовки переменных в этом выражении. У меня x и у ‒ заданные значения, а A и B ‒ те константы, которые находит алгоритм Евклида, D же, соответственно ‒ наибольший общий делитель x и у.
Часто встречает обратный подход, когда A и B задано, а x и y необходимо найти.
Не тестил, но, возможно при значения x, y, близких к 2^63 будет работать неверно, из-за переполнения при операциях вычитания. При меньших значения работает всегда верно.

Есть еще бинарный вариант расширенного алгоритма Евклида. Вроде как в обычном варианте, по утверждению Кнута, бинарный вариант быстрее на 60% классического. Но с того момента много воды утекло, и операции деления, хоть и остаются медленными, но стали все же намного быстрее, чем это было 20 лет назад, поэтому это утверждение требует проверки на современных компьютерах. Тем не менее, вполне допускаю это для обычного алгоритма Евклида, так как он вполне себе компактен и красив.
Но вот та версия расширенного бинарного алгоритма, что мне удалось найти, не так как красива, как обычного. Приведу ее рабочий и протестированный вариант:

function ExtEuclidBin(X, Y : Int64; var A, B : Int64) : Int64;
var
  g : int64;
  u, v : Int64;
  A1, B1, C, D : Int64;
begin
  g := 1;
  while (X and 1 = 0) and (Y and 1 = 0) do
  begin
    X := X shr 1;
    Y := Y shr 1;
    g := g shl 1;
  end;
  u := X;
  v := Y;
  A1 := 1;
  B1 := 0;
  C := 0;
  D := 1;
  repeat
    while (u and 1 = 0) do
    begin
      u := u shr 1;
      if (A1 and 1 = 0) and (B1 and 1 = 0) then
      begin
        A1 := A1 div 2;
        B1 := B1 div 2;
      end
      else
      begin
        A1 := (A1+Y) div 2;
        B1 := (B1-X) div 2;
      end;
    end;
    while (v and 1 = 0) do
    begin
      v := v shr 1;
      if (C and 1 = 0) and (D and 1 = 0) then
      begin
        C := C div 2;
        D := D div 2;
      end
      else
      begin
        C := (C+Y) div 2;
        D := (D-X) div 2;
      end;
    end;
    if u > v then
    begin
      u := u-v;
      A1 := A1-C;
      B1 := B1-D;
    end
    else
    begin
      v := v-u;
      C := C-A1;
      D := D-B1;
    end;
  until (u = 0) or (v=0);
  if u = 0 then
  begin
    Result := v*g;
    A := C;
    B := D;
  end
  else
  begin
    Result := u*g;
    A := A1;
    B := B1;
  end;
end;

Здесь такое же обозначение переменных, как и в предыдущем примере.
Это просто жесть, конечно, по сравнению с компактной рекурсивной записью расширенного алгоритма Евклида в классическом варианте. Как думаете, будет ли такой алгоритм работать быстрее классического?
Кроме сложности алгоритма, у него есть и другие особенности. В частности, линейное диофантово уравнение допускает множество решений, и классический расширенный алгоритм Евклида дает минимальный вариант коэффициентов А и B, в то время как бинарный дает произвольное правильное решение, то есть найденные А и B могут быть не минимальными по абсолютному значению.
В связи с этим у приведенного бинарного алгоритма есть более серьезное ограничение, чем у классического варианта: он гарантированно корректно работает только при x и у < 2^32. Вызвано это возможным переполнением при вычислении A1, B1, C, D. Может быть, и есть вариант, у которого данное ограничение отсутствует, но мне про него не известно.

Кстати, любопытно, но последние версии Delphi не умеют оптимизировать деление знаковых целых на два, хотя, например, Delphi 7 делала это влет.