aureooms/js-modular-arithmetic-big-endian

import assert from 'assert';

import {
    _alloc as n_alloc,
    _zeros as n_zeros,
    _copy as n_copy,
    _mul as n_mul,
    _idivmod as n_idivmod,
    _sub as n_sub,
    _extended_euclidean_algorithm as n_extended_euclidean_algorithm,
} from '@arithmetic-operations-for/naturals-big-endian';

import _redc from './_redc.js';

/**
 *
 * N has no leading zeroes
 *
 * @param {Number} b the radix
 * @param {Array} N number array of length k.
 *
 *
 *
 *
 */
export default function _montgomery(b, N) {
    assert(N.length > 0);
    assert(N[0] !== 0);

    const k = N.length;

    const _2kp1 = 2 * k + 1;
    const _R = n_zeros(_2kp1);
    _R[k] = 1; // B^k

    const [GCD, GCDi, _S, _Si, _M, _1, _2, _3, _4, _5, steps] =
        n_extended_euclidean_algorithm(b, _R, k, _2kp1, N, 0, k);

    // Assert that GCD(R,N) is 1.
    if (GCD.length - GCDi !== 1 || GCD[GCDi] !== 1)
        throw new Error('Montgomery: GCD(R,N) is not 1.');

    const M = n_alloc(k); // M mod R on k words
    if (steps % 2 === 0) {
        // We use _R[0:k] because it is filled with zeros.
        n_sub(b, _R, 0, k, _M, _M.length - k, _M.length, M, 0, k); // _M.length-k is always 1 ?
    } else {
        n_copy(_M, _M.length - k, _M.length, M, 0); // _M.length-k is always 1 ?
    }

    // R^2 mod N
    const _R2 = n_zeros(_2kp1);
    _R2[0] = 1;
    n_idivmod(b, _R2, 0, _2kp1, N, 0, k, _R, 0, _2kp1); // Use mod only function once implemented
    const R2 = n_alloc(k); // R^2 mod N on k words
    n_copy(_R2, k + 1, _2kp1, R2, 0);

    // Avoid using division for the computation of the other constants.
    // From Wikipedia:
    // Performing these operations requires knowing at least N′ and R2 mod N.
    // When R is a power of a small positive integer b, N′ can be computed by
    // Hensel's lemma: The inverse of N modulo b is computed by a naive
    // algorithm (for instance, if b = 2 then the inverse is 1), and Hensel's
    // lemma is used repeatedly to find the inverse modulo higher and higher
    // powers of b, stopping when the inverse modulo R is known; N′ is the
    // negation of this inverse. The constants R mod N and R3 mod N can be
    // generated as REDC(R2 mod N) and as REDC((R2 mod N)(R2 mod N)). The
    // fundamental operation is to compute REDC of a product. When standalone
    // REDC is needed, it can be computed as REDC of a product with 1 mod N. The
    // only place where a direct reduction modulo N is necessary is in the
    // precomputation of R2 mod N.

    // R mod N = REDC(R^2 mod N)
    _redc(b, k, N, 0, k, M, 0, k, _R2, 0, _2kp1);
    const R = n_alloc(k); // R mod N on k words
    n_copy(_R2, k + 1, _2kp1, R, 0);

    // R^3 mod N = REDC((R^2 mod N)(R^2 mod N))
    const _R3 = n_zeros(_2kp1);
    n_mul(b, R2, 0, k, R2, 0, k, _R3, 1, _2kp1);
    _redc(b, k, N, 0, k, M, 0, k, _R3, 0, _2kp1);
    const R3 = n_alloc(k); // R^3 mod N on k words
    n_copy(_R3, k + 1, _2kp1, R3, 0);

    // Console.debug({b, N, k, M, R, R2, R3});
    return {k, M, R, R2, R3};
}