Exponenciación binaria por factorización

Consideremos el problema de calcular $ax^y \pmod{2^d}$, dados enteros $a$, $x$, $y$ y $d \geq 3$, donde $x$ es impar.

El algoritmo de abajo permite resolver este problema con $O(d)$ sumas y operaciones binarias y una sola multiplicación por $y$.

Por la estructura del grupo multiplicativo módulo $2^d$, cualquier número $x$ tal que $x \equiv 1 \pmod 4$ se puede representar como

$$x \equiv b^{L(x)} \pmod{2^d},$$

donde $b \equiv 5 \pmod 8$. Sin pérdida de generalidad asumimos que $x \equiv 1 \pmod 4$, ya que podemos reducir $x \equiv 3 \pmod 4$ a $x \equiv 1 \pmod 4$ sustituyendo $x \mapsto -x$ y $a \mapsto (-1)^{y} a$. En esta notación, $ax^y$ se representa como

$$a x^y \equiv a b^{yL(x)} \pmod{2^d}.$$

La idea central del algoritmo es simplificar el cálculo de $L(x)$ y $b^{y L(x)}$ usando el hecho de que trabajamos módulo $2^d$. Por razones que quedarán claras más adelante, trabajaremos con $4L(x)$ en lugar de $L(x)$, pero tomado módulo $2^d$ en vez de $2^{d-2}$.

En este artículo cubriremos la implementación para enteros de $32$ bits. Sean

• mbin_log_32(r, x) una función que calcula $r+4L(x) \pmod{2^d}$;
• mbin_exp_32(r, x) una función que calcula $r b^{\frac{x}{4}} \pmod{2^d}$;
• mbin_power_odd_32(a, x, y) una función que calcula $ax^y \pmod{2^d}$.

Entonces mbin_power_odd_32 se implementa así:

uint32_t mbin_power_odd_32(uint32_t rem, uint32_t base, uint32_t exp) {
    if (base & 2) {
        /* el divisor se considera negativo */
        base = -base;
        /* comprobar si el resultado debe ser negativo */
        if (exp & 1) {
            rem = -rem;
        }
    }
    return (mbin_exp_32(rem, mbin_log_32(0, base) * exp));
}

Cálculo de 4L(x) a partir de x

Sea $x$ un número impar tal que $x \equiv 1 \pmod 4$. Se puede representar como

$$x \equiv (2^{a_1}+1)\dots(2^{a_k}+1) \pmod{2^d},$$

donde $1 < a_1 < \dots < a_k < d$. Aquí $L(\cdot)$ está bien definido para cada factor, porque son iguales a $1$ módulo $4$. Por lo tanto,

$$4L(x) \equiv 4L(2^{a_1}+1)+\dots+4L(2^{a_k}+1) \pmod{2^{d}}.$$

Así, si precomputamos $t_k = 4L(2^n+1)$ para todo $1 < k < d$, podremos calcular $4L(x)$ para cualquier número $x$.

Para enteros de 32 bits podemos usar la siguiente tabla:

const uint32_t mbin_log_32_table[32] = {
    0x00000000, 0x00000000, 0xd3cfd984, 0x9ee62e18,
    0xe83d9070, 0xb59e81e0, 0xa17407c0, 0xce601f80,
    0xf4807f00, 0xe701fe00, 0xbe07fc00, 0xfc1ff800,
    0xf87ff000, 0xf1ffe000, 0xe7ffc000, 0xdfff8000,
    0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000,
    0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
    0xff000000, 0xfe000000, 0xfc000000, 0xf8000000,
    0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
};

En la práctica se usa un enfoque un poco distinto del descrito arriba. En lugar de encontrar la factorización de $x$, multiplicaremos $x$ sucesivamente por $2^n+1$ hasta convertirlo en $1$ módulo $2^d$. De este modo encontraremos la representación de $x^{-1}$, es decir

$$x (2^{a_1}+1)\dots(2^{a_k}+1) \equiv 1 \pmod {2^d}.$$

Para ello iteramos sobre $n$ tales que $1 < n < d$. Si el $x$ actual tiene el $n$-ésimo bit activado, multiplicamos $x$ por $2^n+1$, lo cual en C++ se hace convenientemente como x = x + (x << n). Esto no cambia los bits menores que $n`, pero pone a cero el$n$-ésimo bit, porque$x$ es impar.

Con todo esto en mente, la función mbin_log_32(r, x) se implementa así:

uint32_t mbin_log_32(uint32_t r, uint32_t x) {
    uint8_t n;
 
    for (n = 2; n < 32; n++) {
        if (x & (1 << n)) {
            x = x + (x << n);
            r -= mbin_log_32_table[n];
        }
    }
 
    return r;
}

Nótese que $4L(x) = -4L(x^{-1})$, así que en lugar de sumar $4L(2^n+1)$, lo restamos de $r$, que inicialmente vale $0$.

Cálculo de x a partir de 4L(x)

Nótese que para $k \geq 1$ se cumple

$$(a 2^{k}+1)^2 = a^2 2^{2k} +a 2^{k+1}+1 = b2^{k+1}+1,$$

de donde (elevando al cuadrado repetidamente) se deduce que

$$(2^a+1)^{2^b} \equiv 1 \pmod{2^{a+b}}.$$

Aplicando este resultado a $a=2^n+1$ y $b=d-k$ deducimos que el orden multiplicativo de $2^n+1$ es un divisor de $2^{d-n}$.

Esto, a su vez, significa que $L(2^n+1)$ debe ser divisible por $2^{n}$, porque el orden de $b$ es $2^{d-2}$ y el orden de $b^y$ es $2^{d-2-v}$, donde $2^v$ es la mayor potencia de $2$ que divide a $y$, así que necesitamos

$$2^{d-k} \equiv 0 \pmod{2^{d-2-v}},$$

por lo tanto $v$ debe ser mayor o igual que $k-2$. Esto es un poco feo y para mitigarlo dijimos al principio que multiplicamos $L(x)$ por $4$. Ahora, si conocemos $4L(x)$, podemos descomponerlo de forma única como suma de $4L(2^n+1)$ revisando sucesivamente los bits de $4L(x)$. Si el $n$-ésimo bit está en $1$, multiplicamos el resultado por $2^n+1$ y reducimos el $4L(x)$ actual en $4L(2^n+1)$.

Así, mbin_exp_32 se implementa de la siguiente forma:

uint32_t mbin_exp_32(uint32_t r, uint32_t x) {
    uint8_t n;
 
    for (n = 2; n < 32; n++) {
        if (x & (1 << n)) {
            r = r + (r << n);
            x -= mbin_log_32_table[n];
        }
    }
 
    return r;
}

Optimizaciones adicionales

Es posible reducir a la mitad el número de iteraciones si se observa que $4L(2^{d-1}+1)=2^{d-1}$ y que para $2k \geq d$ se cumple

$$(2^n+1)^2 \equiv 2^{2n} + 2^{n+1}+1 \equiv 2^{n+1}+1 \pmod{2^d},$$

lo que permite deducir que $4L(2^n+1)=2^n$ para $2n \geq d$. Así, se puede simplificar el algoritmo recorriendo solo hasta $\frac{d}{2}$ y luego usar el hecho anterior para calcular la parte restante con operaciones de bits:

uint32_t mbin_log_32(uint32_t r, uint32_t x) {
    uint8_t n;
 
    for (n = 2; n != 16; n++) {
        if (x & (1 << n)) {
            x = x + (x << n);
            r -= mbin_log_32_table[n];
        }
    }
 
    r -= (x & 0xFFFF0000);
 
    return r;
}
 
uint32_t mbin_exp_32(uint32_t r, uint32_t x) {
    uint8_t n;
 
    for (n = 2; n != 16; n++) {
        if (x & (1 << n)) {
            r = r + (r << n);
            x -= mbin_log_32_table[n];
        }
    }
 
    r *= 1 - (x & 0xFFFF0000);
 
    return r;
}

Cálculo de la tabla de logaritmos

Para calcular la tabla de logaritmos, se puede modificar el algoritmo de Pohlig–Hellman  para el caso en que el módulo es una potencia de $2$.

Nuestra tarea principal aquí es calcular $x$ tal que $g^x \equiv y \pmod{2^d}$, donde $g=5$ e $y$ es un número de la forma $2^n+1$.

Elevando ambos lados al cuadrado $k$ veces llegamos a

$$g^{2^k x} \equiv y^{2^k} \pmod{2^d}.$$

Nótese que el orden de $g$ no es mayor que $2^{d}$ (de hecho, que $2^{d-2}$, pero nos quedaremos con $2^d$ por conveniencia); por lo tanto, usando $k=d-1$ tendremos $g^1$ o $g^0$ en el lado izquierdo, lo que nos permite determinar el bit más pequeño de $x$ comparando $y^{2^k}$ con $g$. Ahora supongamos que $x=x_0 + 2^k x_1$, donde $x_0$ es una parte conocida y $x_1$ todavía no se conoce. Entonces

$$g^{x_0+2^k x_1} \equiv y \pmod{2^d}.$$

Multiplicando ambos lados por $g^{-x_0}$, obtenemos

$$g^{2^k x_1} \equiv (g^{-x_0} y) \pmod{2^d}.$$

Ahora, elevando ambos lados al cuadrado $d-k-1$ veces podemos obtener el siguiente bit de $x$, y eventualmente recuperar todos sus bits.

Referencias

• M30, Hans Petter Selasky, 2009