16.11.2020

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

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

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-битными, то, видимо, бинарный алгоритм снова будет впереди.

4 комментария:

  1. 9.7 млн/сек - это круто, молодцы! Я почему-то считал что длинная арифметика займет куда как больше инструкций, а оказалось все можно написать очень компактно.

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

    Я сначала чуть-чуть потерялся сравнивая pascal с asm, пока не понял что asm он еще и с длинной арифметикой :)

    > паскалевские чуть медленнее

    Я думаю там замедление из-за не умения pascal (да и java) выполнить деление и вычисление остатка одной операцией. В итоге плюс одно умножение и вычитание блокирующие начало следующего деления, что выливается в где-то +32% на piledriver. Длинная же арифметика она выполняется параллельно делению и не блокирует его.

    > Бинарный

    Про бинарный согласен, при каких-то числах он точно будет быстрее, но не при 64-битах.

    ОтветитьУдалить
    Ответы
    1. > Я почему-то считал что длинная арифметика займет куда как больше инструкций
      Если не брать сложение и вычитание, в полноценном виде длинная арифметика побольше команд потребует. У меня ситуация проще - необходимо выполнять деление и умножение 128-битных на 64-битные.
      Полноценное 128-битное умножение будет чуть сложнее, а вот с делением получится реальная засада. Поэтому есть предположение, что на аргументах, выходящих за рамки 64 бит, бинарный снова станет быстрее обычного.

      Удалить
    2. Ага, согласен, на сколько я помню, время деления пропорционально квадрату от числа бит, получается полное время: для классического пропорционально кубу, а для бинарного только квадрату, но константа выше. И где-то при 64 или 128 бит он должен обогнать :)))

      Удалить
    3. Ох, опечатался, при 64 он конечно не обгоняет, хотел написать при 128 или 256 битах.

      Удалить