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

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

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

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

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