src/Montgomery.js
import {
_alloc as n_alloc,
_zeros as n_zeros,
_reset as n_reset,
_copy as n_copy,
jz as njz,
_extended_euclidean_algorithm as n_extended_euclidean_algorithm,
_trim_positive as n_trim_positive,
_sub as n_sub,
convert as nconvert,
} from '@arithmetic-operations-for/naturals-big-endian';
import _iadd from './_iadd.js';
import _isub from './_isub.js';
import _montgomery from './_montgomery.js';
import _mul from './_mul.js';
import _redc from './_redc.js';
import modN from './modN.js';
import modR from './modR.js';
export default class Montgomery {
constructor(b, N) {
const {k, M, R, R2, R3} = _montgomery(b, N);
this.b = b;
this.N = N;
this.k = k;
this.M = M;
this.R = R;
this.R2 = R2;
this.R3 = R3;
// Use shared/pooled memory ?
}
one() {
return this.R;
}
zero() {
return n_zeros(this.k);
}
from(x) {
// Conversion into Montgomery form is done by computing .
// aR mod N = REDC((a mod N)(R^2 mod N))
const _2kp1 = 2 * this.k + 1;
const red = n_zeros(_2kp1); // TODO Use UintXArray ?
const amodN = modN(this.b, this.N, x);
_mul(this.b, this.N, this.M, this.R2, amodN, red);
// TODO many unnecessary copies/alloc can be avoided by
// allowing array offsets in methods.
return modR(this.k, red);
}
out(aRmodN) {
// Conversion out of Montgomery form is done by computing.
// a mod N = REDC(aR mod N)
const _2kp1 = 2 * this.k + 1;
const _red = n_zeros(_2kp1); // TODO Use UintXArray ?
n_copy(aRmodN, 0, this.k, _red, _2kp1 - this.k);
_redc(this.b, this.k, this.N, 0, this.k, this.M, 0, this.k, _red, 0, _2kp1);
const i = n_trim_positive(_red, this.k + 1, _2kp1);
const red = n_alloc(_2kp1 - i); // TODO Use UintXArray ?
n_copy(_red, i, _2kp1, red, 0);
return red;
}
mul(aRmodN, bRmodN) {
const _2kp1 = 2 * this.k + 1;
const abRmodN = n_zeros(_2kp1);
_mul(this.b, this.N, this.M, aRmodN, bRmodN, abRmodN);
return modR(this.k, abRmodN);
}
add(aRmodN, bRmodN) {
const aRpbRmodN = n_alloc(this.k);
n_copy(aRmodN, 0, this.k, aRpbRmodN, 0);
_iadd(this.b, this.N, aRpbRmodN, bRmodN);
return aRpbRmodN;
}
sub(aRmodN, bRmodN) {
const aRpbRmodN = n_alloc(this.k);
n_copy(aRmodN, 0, this.k, aRpbRmodN, 0);
_isub(this.b, this.N, aRpbRmodN, bRmodN);
return aRpbRmodN;
}
inv(aRmodN) {
// The modular inverse
// Compute (aR mod N)^-1 using Euclidean algo
const ai = n_trim_positive(aRmodN, 0, this.k);
let [GCD, GCDi, _S, _Si, aRmodNi, _1, _2, _3, _4, _5, steps] =
n_extended_euclidean_algorithm(
this.b,
this.N,
0,
this.k,
aRmodN,
ai,
this.k,
);
// Assert that GCD(N,aRmodN) is 1.
if (GCD.length - GCDi !== 1 || GCD[GCDi] !== 1)
throw new Error('aRmodN has no inverse modulo N');
const _2kp1 = 2 * this.k + 1;
const red = n_zeros(_2kp1); // TODO Use UintXArray ?
if (steps % 2 === 1) {
// We compute N - aRmodNi
const temporary = n_zeros(this.k);
n_sub(
this.b,
this.N,
0,
this.k,
aRmodNi,
0,
this.k,
temporary,
0,
this.k,
);
aRmodNi = temporary;
}
// A^-1 R mod N = REDC((aR mod N)^-1(R^3 mod N)).
_mul(this.b, this.N, this.M, this.R3, aRmodNi, red);
return modR(this.k, red);
}
pown(aRmodN, x) {
// Modular
// exponentiation can be done using exponentiation by squaring by initializing the
// initial product to the Montgomery representation of 1, that is, to R mod N, and
// by replacing the multiply and square steps by Montgomery multiplies.
const nonneg = x >= 0;
if (!nonneg) x = -x;
if (x === 0) return this.R;
if (x === 1) return nonneg ? aRmodN : this.inv(aRmodN);
const xbits = [];
do {
xbits.push(x & 1); // eslint-disable-line no-bitwise
x >>= 1; // eslint-disable-line no-bitwise
} while (x !== 1);
return this._powb(aRmodN, xbits, nonneg);
}
_powb(aRmodN, xbits, nonneg) {
// The binary expansion of the exponent is 1 concatenanted with xbits
// reversed. Must have xbits.length >= 1.
const aRmodNpown = n_alloc(this.k);
n_copy(aRmodN, 0, this.k, aRmodNpown, 0);
const _2kp1 = 2 * this.k + 1;
const temporary = n_alloc(_2kp1);
do {
n_reset(temporary, 0, _2kp1);
_mul(this.b, this.N, this.M, aRmodNpown, aRmodNpown, temporary);
n_copy(temporary, _2kp1 - this.k, _2kp1, aRmodNpown, 0);
if (xbits.pop() === 1) {
n_reset(temporary, 0, _2kp1);
_mul(this.b, this.N, this.M, aRmodNpown, aRmodN, temporary);
n_copy(temporary, _2kp1 - this.k, _2kp1, aRmodNpown, 0);
}
} while (xbits.length > 0);
return nonneg ? aRmodNpown : this.inv(aRmodNpown);
}
pow(aRmodN, b, nonneg = true) {
if (njz(b, 0, b.length - 1)) {
// B consists of a single limb
return this.pown(aRmodN, nonneg ? b.at(-1) : -b.at(-1));
}
const xbits = nconvert(this.b, 2, b, 0, b.length);
xbits.reverse();
xbits.pop();
return this._powb(aRmodN, xbits, nonneg);
}
}