src/core/arithmetic/gcd/_extended_euclidean_algorithm_loop.js
import assert from 'assert';
import increment from '../../../api/arithmetic/add/increment.js';
import _idivmod from '../../../api/arithmetic/div/_idivmod.js';
import mul from '../../../api/arithmetic/mul/mul.js';
import ge from '../../../api/compare/ge.js';
import _copy from '../../array/_copy.js';
import _reset from '../../array/_reset.js';
import _trim_positive from '../../convert/_trim_positive.js';
import _iadd from '../add/_iadd.js';
/**
* Extended Euclidean algorithm.
*
* @see https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
*
* Given two input integers a and b, a > b.
* Let r_0 = a, r_1 = b,
* s_0 = 1, s_1 = 0,
* t_0 = 0, t_1 = 1 (Fibonacci :winkemoji:).
*
* Let r_{i+1} = r_{i-1} % r_i (division remainder)
* Let q_i = (r_{i-1} - r_{i+1}) / r_i > 0 (division quotient)
*
* An alternative definition is
* r_{i+1} = r_{i-1} - q_i r_i (with 0 <= r_{i+1} < r_i)
*
* Then define
* s_{i+1} = s_{i-1} - q_i s_i
* t_{i+1} = t_{i-1} - q_i t_i
*
* Let k be such that r_i > 0 for all i <= k and r_{k+1} = 0.
*
* Since q_i > 0, if t_{i-1} > 0 and t_i < 0 then t_{i+1} > 0. On the other
* hand, if t_{i-1} < 0 and t_i > 0 then t_{i+1} < 0. Note that t_2 < 0, so the
* signs of the t_i alternate: t_0 = 0, t_1 = 1, t_2 < 0, t_3 > 0, t_4 < 0, t_5
* > 0, etc. The pattern for t is 0, 1, -, +, -, +, -, etc.
*
* The same holds for the s_i with flipped signs: note that s_2 = 1, then
* s_0 = 1, s_1 = 0, s_2 = 1, s_3 < 0, s_4 > 0, s_5 < 0, etc.
* The pattern for s is 1, 0, 1, -, +, -, +, etc.
*
* | 0 1 2 3 4 5 6
* ------------------------
* s | 1 0 1 - + - +
* t | 0 1 - + - + -
*
* With i >= 1, |t_{i+1}| >= |t_i| + |t_{i-1}| > |t_i| because q_i >= 1 and
* |t_i| > 0 and, since the signs are alternating, we have |t_{i-1} - q_i t_i|
* = |t_{i-1}| + |q_i t_i|. Same goes for s_{i+1} w.r.t s_i and s_{i-1}.
*
*
* Input
* -----
* - No leading zeroes
* - R0 >= R1
* - (two bullets above imply |R1| >= |R0|)
* - More at _extended_euclidean_algorithm_allocate.
*
*
* Implementation details
* ----------------------
*
* Here the implementation avoids computing with negative numbers and will only
* output the absolute value of s_{i} and t_{i}. The signs can be recovered
* from the number of steps of the algorithm.
*
* Note that s_i <= b and t_i <= a (proof?). Equality holds
* when i = k+1 and a and b are coprime. Thus s_i can sit in an array as large
* as b and same for t_i and a.
*
* Output
* ------
*
* Calling this function returns the number of steps that were executed.
*
*/
export default function _extended_euclidean_algorithm_loop(
r,
R0,
R1,
S0,
T0,
S1,
T1,
Q,
X,
) {
assert(r >= 2);
const m = R0.length;
const n = R1.length;
let R0i = 0;
const R0j = m;
let R1i = 0;
const R1j = n;
let S0i = S0.length - 1;
const S0j = S0.length;
let T0i = m;
const T0j = m;
let S1i = n;
const S1j = n;
let T1i = m - 1;
const T1j = m;
let Qi = 0;
const Qj = m;
let Xi = 0;
const Xj = 2 * m;
assert(R0i >= 0 && R0j <= R0.length);
assert(R0j - R0i <= 0 || R0[R0i] !== 0);
assert(R1i >= 0 && R1j <= R1.length);
assert(R1j - R1i <= 0 || R1[R1i] !== 0);
assert(ge(R0, R0i, R0j, R1, R1i, R1j));
// We handle the first two steps outside of loop because s_1 = t_0 = 0
// and s_1 = 0, s_2 = 1
// Invariants
// ----------
//
// 1. No leading zeros in R0
// 2. No leading zeros in R1
// 3. |Q| = |R0| (why ???)
assert(Qj - Qi === R0j - R0i);
// 4. s_0 = S0 > 0
// 5. s_1 = S1 < 0
// 6. t_0 = T0 < 0
// 7. t_1 = T1 > 0
if (R1i === R1j) return [R0i, S0i, T0i, S1i, T1i, 1];
// Q_1 = (r_0 - r_2) / r_1
// R0 is r_0 and becomes r_2
// R1 is r_1
// Q is q_1
_idivmod(r, R0, R0i, R0j, R1, R1i, R1j, Q, Qi, Qj);
// Remove leading zeros from Q
// since Q = R0 / R1 we have |R0| - |R1| <= |Q| <= |R0| - |R1| + 1
Qi = Qj - (R0j - R1j + 1); // R0i = R1i = 0
if (Q[Qi] === 0) ++Qi;
assert(Qi < Qj && Q[Qi] !== 0);
// Remove leading zeros from R0
// since R0 = R0 % R1 we have |R0| <= |R1|
R0i = _trim_positive(R0, R0j - (R1j - R1i), R0j);
// S_2 = s_0 - q_1 * s_1 = s_0
// S0 is s_0 and becomes s_2 i.e. NOTHING TO DO
// t_2 = t_0 - q_1 * t_1 = q_1
// T0 is t_0 and becomes t_2
T0i = T0j - (Qj - Qi);
_copy(Q, Qi, Qj, T0, T0i);
// Invariants
// ----------
//
// 1. No leading zeros in R0
assert(R0i >= 0 && R0j <= R0.length);
assert(R0j - R0i <= 0 || R0[R0i] !== 0);
// 2. No leading zeros in R1
assert(R1i >= 0 && R1j <= R1.length);
assert(R1j - R1i <= 0 || R1[R1i] !== 0);
// 3. |Q| = |R1| (why ???)
// assert(Qj - Qi === R1j - R1i); // NOT TRUE !
// 4. s_1 = S1 < 0
// 5. s_2 = S0 > 0
// 6. t_1 = T1 > 0
// 7. t_2 = T0 < 0
if (R0i === R0j) return [R1i, S0i, T0i, S1i, T1i, 2];
// Q_2 = (r_1 - r_{i+1}) / r_2
// R1 is r_1 and becomes r_3
// R0 is r_2
// Q is q_2
Qi = Qj - (R1j - R1i);
_reset(Q, Qi, Qj);
_idivmod(r, R1, R1i, R1j, R0, R0i, R0j, Q, Qi, Qj);
// Remove leading zeros from Q
// since Q = R1 / R0 we have |R1| - |R0| <= |Q| <= |R1| - |R0| + 1
Qi = Qj - (R1j - R0j + R0i + 1); // R1i = 0
if (Q[Qi] === 0) ++Qi;
assert(Qi < Qj && Q[Qi] !== 0);
// Remove leading zeros from R1
// since R1 = R1 % R0 we have |R1| <= |R0|
R1i = _trim_positive(R1, R1j - (R0j - R0i), R1j);
// S_3 = s_1 - q_2 * s_2 = -q_2
S1i = S1j - (Qj - Qi);
_copy(Q, Qi, Qj, S1, S1i);
// Q_2 * t_2
// since Q and T0 have no leading zeros then
// Q * T0 has |Q| + |T0| - 1 <= |Q*T0| <= |Q| + |T0| limbs with no leading zeros.
Xi = Xj - (Qj - Qi) - (T0j - T0i);
mul(r, T0, T0i, T0j, Q, Qi, Qj, X, Xi, Xj);
// T_3 = t_1 - q_2 * t_2 = 1 - q_2 * t_2
// T1 is t_1 and becomes t_3
increment(r, X, Xi, Xj);
Xi = _trim_positive(X, Xi, Xj);
T1i = T1j - (Xj - Xi);
_copy(X, Xi, Xj, T1, T1i);
let steps = 3;
while (true) {
// Invariants
// ----------
//
// 1. No leading zeros in R0
assert(R0i >= 0 && R0j <= R0.length);
assert(R0j - R0i <= 0 || R0[R0i] !== 0);
// 2. No leading zeros in R1
assert(R1i >= 0 && R1j <= R1.length);
assert(R1j - R1i <= 0 || R1[R1i] !== 0);
// 3. |Q| = |R0| (why ???)
// 4. s_{i-1} = S0 > 0
// 5. s_i = S1 < 0
// 6. t_{i-1} = T0 < 0
// 7. t_i = T1 > 0
if (R1i === R1j) return [R0i, S0i, T0i, S1i, T1i, steps];
++steps;
// Q_i = (r_{i-1} - r_{i+1}) / r_i
// R0 is r_{i-1} and becomes r_{i+1}
// R1 is r_i
// Q is q_i
Qi = Qj - (R0j - R0i);
_reset(Q, Qi, Qj);
_idivmod(r, R0, R0i, R0j, R1, R1i, R1j, Q, Qi, Qj);
// Remove leading zeros from Q
// since Q = R0 / R1 we have |R0| - |R1| <= |Q| <= |R0| - |R1| + 1
Qi = Qj - (R0j - R0i - R1j + R1i + 1);
if (Q[Qi] === 0) ++Qi;
assert(Qi < Qj && Q[Qi] !== 0);
// Remove leading zeros from R0
// since R0 = R0 % R1 we have |R0| <= |R1|
R0i = _trim_positive(R0, R0j - (R1j - R1i), R0j);
// Q_i * s_i
// since Q and S1 have no leading zeros then
// Q * S1 has |Q| + |S1| - 1 <= |Q*S1| <= |Q| + |S1| limbs with no leading zeros.
Xi = Xj - (Qj - Qi) - (S1j - S1i);
_reset(X, Xi, Xj);
mul(r, Q, Qi, Qj, S1, S1i, S1j, X, Xi, Xj);
if (X[Xi] === 0) ++Xi; // Remove leading zero if no carry
// s_{i+1} = s_{i-1} - q_i * s_i
// S0 is s_{i-1} and becomes s_{i+1}
S0i = S0j - (Xj - Xi + 1);
S0i = Math.max(0, S0i); // Next addition never overflows beyond bounds
_iadd(r, S0, S0i, S0j, X, Xi, Xj);
if (S0[S0i] === 0) ++S0i;
// Q_i * t_i
// since Q and T1 have no leading zeros then
// Q * T1 has |Q| + |T1| - 1 <= |Q*T1| <= |Q| + |T1| limbs with no leading zeros.
Xi = Xj - (Qj - Qi) - (T1j - T1i);
_reset(X, Xi, Xj);
mul(r, Q, Qi, Qj, T1, T1i, T1j, X, Xi, Xj);
if (X[Xi] === 0) ++Xi; // Remove leading zero if no carry
// t_{i+1} = t_{i-1} - q_i * t_i
// T0 is t_{i-1} and becomes t_{i+1}
T0i = T0j - (Xj - Xi + 1);
T0i = Math.max(0, T0i); // Next addition never overflows beyond bounds
_iadd(r, T0, T0i, T0j, X, Xi, Xj);
if (T0[T0i] === 0) ++T0i;
// Invariants
// ----------
//
// 1. No leading zeros in R0
assert(R0i >= 0 && R0j <= R0.length);
assert(R0j - R0i <= 0 || R0[R0i] !== 0);
// 2. No leading zeros in R1
assert(R1i >= 0 && R1j <= R1.length);
assert(R1j - R1i <= 0 || R1[R1i] !== 0);
// 3. |Q| = |R1| (why ???)
// 4. s_{i-1} = S1 < 0
// 5. s_i = S0 > 0
// 6. t_{i-1} = T1 > 0
// 7. t_i = T0 < 0
if (R0i === R0j) return [R1i, S0i, T0i, S1i, T1i, steps];
++steps;
// Q_i = (r_{i-1} - r_{i+1}) / r_i
// R1 is r_{i-1} and becomes r_{i+1}
// R0 is r_i
// Q is q_i
Qi = Qj - (R1j - R1i);
_reset(Q, Qi, Qj);
_idivmod(r, R1, R1i, R1j, R0, R0i, R0j, Q, Qi, Qj);
// Remove leading zeros from Q
// since Q = R1 / R0 we have |R1| - |R0| <= |Q| <= |R1| - |R0| + 1
Qi = Qj - (R1j - R1i - R0j + R0i + 1);
if (Q[Qi] === 0) ++Qi;
assert(Qi < Qj && Q[Qi] !== 0);
// Remove leading zeros from R1
// since R1 = R1 % R0 we have |R1| <= |R0|
R1i = _trim_positive(R1, R1j - (R0j - R0i), R1j);
// Q_i * s_i
// since Q and S0 have no leading zeros then
// Q * S0 has |Q| + |S0| - 1 <= |Q*S0| <= |Q| + |S0| limbs with no leading zeros.
Xi = Xj - (Qj - Qi) - (S0j - S0i);
_reset(X, Xi, Xj);
mul(r, Q, Qi, Qj, S0, S0i, S0j, X, Xi, Xj);
if (X[Xi] === 0) ++Xi; // Remove leading zero if no carry
// s_{i+1} = s_{i-1} - q_i * s_i
// S1 is s_{i-1} and becomes s_{i+1}
S1i = S1j - (Xj - Xi + 1);
S1i = Math.max(0, S1i); // Next addition never overflows beyond bounds
_iadd(r, S1, S1i, S1j, X, Xi, Xj);
if (S1[S1i] === 0) ++S1i;
// Q_i * t_i
// since Q and T0 have no leading zeros then
// Q * T0 has |Q| + |T0| - 1 <= |Q*T0| <= |Q| + |T0| limbs with no leading zeros.
Xi = Xj - (Qj - Qi) - (T0j - T0i);
_reset(X, Xi, Xj);
mul(r, Q, Qi, Qj, T0, T0i, T0j, X, Xi, Xj);
if (X[Xi] === 0) ++Xi; // Remove leading zero if no carry
// t_{i+1} = t_{i-1} - q_i * t_i
// T1 is t_{i-1} and becomes t_{i+1}
T1i = T1j - (Xj - Xi + 1);
T1i = Math.max(0, T1i); // Next addition never overflows beyond bounds
_iadd(r, T1, T1i, T1j, X, Xi, Xj);
if (T1[T1i] === 0) ++T1i;
}
}