Multiplicación de Montgomery

Muchos algoritmos de teoría de números, como el test de primalidad o la factorización de enteros, y en criptografía, como RSA, requieren muchas operaciones módulo un número grande. Una multiplicación como $x y \bmod{n}$ es bastante lenta de calcular con los algoritmos típicos, porque hace falta una división para saber cuántas veces hay que restar $n$ del producto. Y la división es una operación realmente costosa, sobre todo con números grandes.

La multiplicación (modular) de Montgomery es un método que permite calcular esas multiplicaciones más rápido. En lugar de dividir el producto y restar $n$ varias veces, suma múltiplos de $n$ para anular los bits bajos y luego simplemente descarta los bits bajos.

Representación de Montgomery

Sin embargo, la multiplicación de Montgomery no es gratis. El algoritmo solo funciona en el espacio de Montgomery. Y hay que transformar nuestros números a ese espacio antes de poder empezar a multiplicar.

Para el espacio necesitamos un entero positivo $r \ge n$ coprimo con $n$, es decir, $\gcd(n, r) = 1$. En la práctica siempre elegimos $r$ como $2^m$ para un entero positivo $m$, porque entonces las multiplicaciones, divisiones y operaciones módulo $r$ se pueden implementar de forma eficiente con desplazamientos y otras operaciones de bits. $n$ será un número impar en prácticamente todas las aplicaciones, porque no es difícil factorizar un número par. Así que toda potencia de $2$ será coprima con $n$.

El representante $\bar{x}$ de un número $x$ en el espacio de Montgomery se define como:

$$\bar{x} := x \cdot r \bmod n$$

Nótese que la transformación es, de hecho, justamente una multiplicación de las que queremos optimizar. Así que sigue siendo una operación costosa. Sin embargo, solo hace falta transformar un número una vez al espacio. En cuanto estamos en el espacio de Montgomery, podemos realizar tantas operaciones como queramos de forma eficiente. Y al final transformamos el resultado final de vuelta. Así que, mientras se hagan muchas operaciones módulo $n$, esto no será un problema.

Dentro del espacio de Montgomery todavía se pueden realizar la mayoría de las operaciones como de costumbre. Se pueden sumar dos elementos ($x \cdot r + y \cdot r \equiv (x + y) \cdot r \bmod n$), restar, comprobar igualdad e incluso calcular el máximo común divisor de un número con $n$ (porque $\gcd(n, r) = 1$). Todo con los algoritmos habituales.

Sin embargo, esto no vale para la multiplicación.

Esperamos que el resultado sea:

$$\bar{x} * \bar{y} = \overline{x \cdot y} = (x \cdot y) \cdot r \bmod n.$$

Pero la multiplicación normal nos da:

$$\bar{x} \cdot \bar{y} = (x \cdot y) \cdot r \cdot r \bmod n.$$

Por lo tanto, la multiplicación en el espacio de Montgomery se define como:

$$\bar{x} * \bar{y} := \bar{x} \cdot \bar{y} \cdot r^{-1} \bmod n.$$

Reducción de Montgomery

La multiplicación de dos números en el espacio de Montgomery requiere un cálculo eficiente de $x \cdot r^{-1} \bmod n$. Esta operación se llama reducción de Montgomery, y también se conoce como el algoritmo REDC.

Como $\gcd(n, r) = 1$, sabemos que existen dos números $r^{-1}$ y $n^{\prime}$ con $0 < r^{-1}, n^{\prime} < n$ tales que

$$r \cdot r^{-1} + n \cdot n^{\prime} = 1.$$

Tanto $r^{-1}$ como $n^{\prime}$ se pueden calcular usando el algoritmo de Euclides extendido.

Usando esta identidad podemos escribir $x \cdot r^{-1}$ como:

$$\begin{array}{rl}{\displaystyle x\cdot r^{-1}}& {\displaystyle =x\cdot r\cdot r^{-1}\mathrm{/}r=x\cdot (-n\cdot n^{\mathrm{\prime }}+1)\mathrm{/}r}\\ {\displaystyle }& {\displaystyle =(-x\cdot n\cdot n^{\mathrm{\prime }}+x)\mathrm{/}r\equiv (-x\cdot n\cdot n^{\mathrm{\prime }}+l\cdot r\cdot n+x)\mathrm{/}r<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>n}\\ {\displaystyle }& {\displaystyle \equiv ((-x\cdot n^{\mathrm{\prime }}+l\cdot r)\cdot n+x)\mathrm{/}r<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>n}\end{array}$$$$\begin{aligned} x \cdot r^{-1} &amp;= x \cdot r \cdot r^{-1} / r = x \cdot (-n \cdot n^{\prime} + 1) / r \ &amp;= (-x \cdot n \cdot n^{\prime} + x) / r \equiv (-x \cdot n \cdot n^{\prime} + l \cdot r \cdot n + x) / r \bmod n\ &amp;\equiv ((-x \cdot n^{\prime} + l \cdot r) \cdot n + x) / r \bmod n \end{aligned}$$

Las equivalencias valen para cualquier entero arbitrario $l$. Esto significa que podemos sumar o restar un múltiplo arbitrario de $r$ a $x \cdot n^{\prime}$, o en otras palabras, podemos calcular $q := x \cdot n^{\prime}$ módulo $r$.

Esto nos da el siguiente algoritmo para calcular $x \cdot r^{-1} \bmod n$:

function reduce(x):
    q = (x mod r) * n' mod r
    a = (x - q * n) / r
    if a < 0:
        a += n
    return a

Como $x < n \cdot n < r \cdot n$ (incluso si $x$ es el producto de una multiplicación) y $q \cdot n < r \cdot n$, sabemos que $-n < (x - q \cdot n) / r < n$. Por lo tanto, la operación de módulo final se implementa con una sola comprobación y una suma.

Como se ve, podemos realizar la reducción de Montgomery sin ninguna operación de módulo pesada. Si elegimos $r$ como una potencia de $2$, las operaciones de módulo y las divisiones del algoritmo se pueden calcular con máscaras de bits y desplazamientos.

Una segunda aplicación de la reducción de Montgomery es transferir un número de vuelta del espacio de Montgomery al espacio normal.

Truco del inverso rápido

Para calcular el inverso $n^{\prime} := n^{-1} \bmod r$ de forma eficiente, podemos usar el siguiente truco (inspirado en el método de Newton):

$$a \cdot x \equiv 1 \bmod 2^k \Longrightarrow a \cdot x \cdot (2 - a \cdot x) \equiv 1 \bmod 2^{2k}$$

Esto se puede demostrar fácilmente. Si tenemos $a \cdot x = 1 + m \cdot 2^k$, entonces:

$$\begin{array}{rl}{\displaystyle a\cdot x\cdot (2-a\cdot x)}& {\displaystyle =2\cdot a\cdot x-(a\cdot x{)}^2}\\ {\displaystyle }& {\displaystyle =2\cdot (1+m\cdot 2^k)-(1+m\cdot 2^k{)}^2}\\ {\displaystyle }& {\displaystyle =2+2\cdot m\cdot 2^k-1-2\cdot m\cdot 2^k-m^2\cdot 2^{2k}}\\ {\displaystyle }& {\displaystyle =1-m^2\cdot 2^{2k}}\\ {\displaystyle }& {\displaystyle \equiv 1<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>2^{2k}\mathrm{.}}\end{array}$$$$\begin{aligned} a \cdot x \cdot (2 - a \cdot x) &amp;= 2 \cdot a \cdot x - (a \cdot x)^2 \ &amp;= 2 \cdot (1 + m \cdot 2^k) - (1 + m \cdot 2^k)^2 \ &amp;= 2 + 2 \cdot m \cdot 2^k - 1 - 2 \cdot m \cdot 2^k - m^2 \cdot 2^{2k} \ &amp;= 1 - m^2 \cdot 2^{2k} \ &amp;\equiv 1 \bmod 2^{2k}. \end{aligned}$$

Esto significa que podemos empezar con $x = 1$ como inverso de $a$ módulo $2^1$, aplicar el truco unas cuantas veces y en cada iteración duplicamos la cantidad de bits correctos de $x$.

Implementación

Con el compilador GCC todavía podemos calcular $x \cdot y \bmod n$ de forma eficiente cuando los tres números son enteros de 64 bits, porque el compilador soporta enteros de 128 bits con los tipos __int128 y __uint128.

long long result = (__int128)x * y % n;

Sin embargo, no hay un tipo para enteros de 256 bits. Por eso aquí mostramos una implementación para una multiplicación de 128 bits.

using u64 = uint64_t;
using u128 = __uint128_t;
using i128 = __int128_t;
 
struct u256 {
    u128 high, low;
 
    static u256 mult(u128 x, u128 y) {
        u64 a = x >> 64, b = x;
        u64 c = y >> 64, d = y;
        // (a*2^64 + b) * (c*2^64 + d) =
        // (a*c) * 2^128 + (a*d + b*c)*2^64 + (b*d)
        u128 ac = (u128)a * c;
        u128 ad = (u128)a * d;
        u128 bc = (u128)b * c;
        u128 bd = (u128)b * d;
        u128 carry = (u128)(u64)ad + (u128)(u64)bc + (bd >> 64u);
        u128 high = ac + (ad >> 64u) + (bc >> 64u) + (carry >> 64u);
        u128 low = (ad << 64u) + (bc << 64u) + bd;
        return {high, low};
    }
};
 
struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;
    }
 
    u128 init(u128 x) {
        x %= mod;
        for (int i = 0; i < 128; i++) {
            x <<= 1;
            if (x >= mod)
                x -= mod;
        }
        return x;
    }
 
    u128 reduce(u256 x) {
        u128 q = x.low * inv;
        i128 a = x.high - u256::mult(q, mod).high;
        if (a < 0)
            a += mod;
        return a;
    }
 
    u128 mult(u128 a, u128 b) {
        return reduce(u256::mult(a, b));
    }
 
    u128 mod, inv;
};

Transformación rápida

El método actual de transformar un número al espacio de Montgomery es bastante lento. Hay formas más rápidas.

Se puede observar la siguiente relación:

$$\bar{x} := x \cdot r \bmod n = x \cdot r^2 / r = x * r^2$$

Transformar un número al espacio es simplemente una multiplicación dentro del espacio del número por $r^2$. Por lo tanto, podemos precomputar $r^2 \bmod n$ y realizar una multiplicación en lugar de desplazar el número 128 veces.

En el código siguiente inicializamos r2 con -n % n, que es equivalente a $r - n \equiv r \bmod n$, y lo desplazamos 4 veces para obtener $r \cdot 2^4 \bmod n$. Este número se puede interpretar como $2^4$ en el espacio de Montgomery. Si lo elevamos al cuadrado $5$ veces, obtenemos $(2^4)^{2^5} = (2^4)^{32} = 2^{128} = r$ en el espacio de Montgomery, que es exactamente $r^2 \bmod n$.

struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1), r2(-n % n) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;

        for (int i = 0; i < 4; i++) {
            r2 <<= 1;
            if (r2 >= mod)
                r2 -= mod;
        }
        for (int i = 0; i < 5; i++)
            r2 = mul(r2, r2);
    }

    u128 init(u128 x) {
        return mult(x, r2);
    }

    u128 mod, inv, r2;
};