Обратный элемент

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

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

Итак, наиболее часто встречаемый вариант на целых числах основан на решении итерационным методом диофантова уравнения первого порядка:

Ax + By = R

Здесь A и B ‒ числа, для которых мы ищем частное решение (если R=0, то мы также можем найти наибольший общий делитель, который будет равен значению R с предыдущей итерации метода Евклида, если R=1, то мы находим как раз обратный элемент в кольце вычетов). А x и y ‒ это неизвестные целые коэффициенты, которые нужно найти, что бы выполнялось равенство.
В прошлых постах я несколько нестандартно трактовал их в обратном порядке, когда x и y были числами, а A и B ‒ коэффициентами. Впрочем, так как они вполне симметричны, то суть дела это не меняет.

Для поиска обратного элемента удобнее в качестве A взять модуль, а в качестве B ‒ число, для которого мы ищем обратный элемент.
Тогда после решения диофантова уравнения Ax + By = 1 y будет обратным элементом к B.
Не к каждому B в кольце вычетом по модулю A может быть найден элемент. Это возможно только когда A и B является взаимно простыми числами. В противном случае в результате работы расширенного алгоритма Евклида мы скорее всего получим решение Ax + By = 0, т. е. R с предыдущей итерации будет НОД(A, B), а By mod A  = 0.

В алгоритме с целыми числами знак коэффициентов на каждой итерации меняется, поэтому если y получился меньше ноля, то в качестве обратного элемента к B мы должны взять A + y.

Вариант для натуральных чисел отличается от целых лишь уравнением:

Ax - By = R

Но если в классическом алгоритме на каждой итерации менялся знак элементов, то для натуральных чисел меняются сами элементы местами.
То есть если на первой итерации

Ax₁ - By₁ = R₁,
где  x₁  = 1; y₁ = Q₁ = A div B; R₁ = A mod B,

То на второй

By₂ - Ax₂ = R₂,

где y₂ и x₂ считаются аналогично первой итерации.
Здесь тоже есть тонкость: так как на каждой нечетной итерации коэффициент при B берется со знаком "минус", то обратный элемент должен считаться так: A - y, в то время как на четных итерациях он точно равен y.
Для расчета обратного элемента коэффициент x нам не нужен, поэтому его и вовсе можно не считать.
А для y можно использовать следующую рекуррентную формулу:

Кстати, x считается точно по такой же формуле, но начальные значения поменяны местами.

В результате на чистом паскале получается так:

function getInvElNat(Module, El : UInt64) : UInt64;
var
  b : byte;
  Y1, Y0 : UInt64;
  D, A, R : UInt64;
begin
  A := Module;
  b := 0;
  Y1 := 0;
  Y0 := 1;
  while El > 1 do
  begin
    D := A div El;
    R := A - El*D;
    A := El;
    El := R;
    Result := Y0*D + Y1;
    b := not(b);
    Y1 := Y0;
    Y0 := Result;
  end;
  if b <> 0 then
    Result := Module-Result;
end; 

По хорошему b должна быть типа boolean, но с точки зрения производительности разницы нет.

А вот так я переписал на ассемблере:

function getInvElNatASM(Module, El : UInt64) : UInt64;
// Module - RCX, El - RDX, Result - RAX
asm
  push RSI;
  mov RAX, Module;
  mov R8, El;
  mov R9, El;
  xor R10b, R10b; // b
  xor R11, R11; //Y0
  xor RSI, RSI; // Y1;
  inc R11;
  cmp R8, small 1;
  jbe @Finish;
@Loop:
  xor RDX, RDX;
  div R8;
  mov R9, R11;
  imul R9, RAX;
  add R9, RSI;
  not R10b;
  mov RAX, R8;
  mov R8, RDX;
  mov RSI, R11;
  mov R11, R9;
  cmp R8, small 1;
  ja @Loop;
@Finish:
  sub Module, R9;
  test R10b, R10b;
  cmovz RAX, R9;
  cmovnz RAX, Module;
  pop RSI;
end;

Получилось гораздо компактнее, чем целочисленный вариант с длинной арифметикой. Правда, производительность отличается не сильно:
‒ целочисленный вариант на asm ‒ 10.7 млн. оп./с;
‒ натуральный вариант на asm ‒ 11.4 млн. оп./с;
‒ натуральный вариант на Pascal ‒ 7.8 млн. оп./с.