Implementación de algoritmos de teoría de números/Exponenciación modular

Una exponenciación modular calcula el residuo cuando un número entero positivo b (la base) se eleva a la e-ésima potencia (el exponente), b^e, y es dividido por el entero positivo m, llamado módulo. En notación matemática, dada la base b, el exponente e, y el módulo m, la exponenciación modular c se escribe:

$c\equiv b^{e}{\pmod {m}}$

Es particularmente útil en ciencias de la computación, especialmente en el campo de la criptografía.

La exponenciación modular se puede realizar con exponente negativo e encontrando el inverso multiplicativo modular d de b módulo m usando el algoritmo extendido de Euclides. Esto es:

$c\equiv b^{e}\equiv d^{\left|e\right|}{\pmod {m}}$ donde $e<0$ y $b\cdot d\equiv 1{\pmod {m}}$

Problemas de exponenciación modular similares al descrito arriba son considerados fáciles de resolver, incluso cuando los números que se manejan son enormes.

Por otro lado, el cálculo del logaritmo discreto — es decir, la tarea de encontrar el exponente e si es dado un b, c, y m — es un problema de los considerados difíciles. Este comportamiento de función unidireccional hace a la exponenciación modular un candidato para su uso en algoritmos criptográficos.

Método de memoria eficiente

Este requiere $O(e)$ multiplicaciones para completarse. En pseudocódigo, este método puede realizarse de la siguiente manera:

function modular_pow(base, exponent, modulus) is
    if modulus = 1 then return 0 then
    c := 1
    for e_prime = 0 to exponent-1 do
        c := (c * base) mod modulus
    return c

Método binario derecha a izquierda

Este método reduce drásticamente el número de operaciones para realizar la exponenciación modular, mientras se mantiene la misma cantidad de memoria que en el método anterior. Es una combinación del método anterior y un principio más general llamado exponenciación por cuadrados (más conocido como exponenciación binaria).

Primero, se requiere que el exponente e sea convertido a notación binaria. Esto es, e puede ser escrito como:

$e=\sum _{i=0}^{n-1}a_{i}2^{i}$

En tal notación, la longitud de e es n bits. a_i puede tomar el valor de 0 o 1 para cualquier i tal que 0 ≤ i < n. Por definición, a_{n − 1} = 1.

El valor b^e puede ser escrito entonces como:

$b^{e}=b^{\left(\sum _{i=0}^{n-1}a_{i}2^{i}\right)}=\prod _{i=0}^{n-1}\left(b^{2^{i}}\right)^{a_{i}}$

La solución c es por lo tanto:

$c\equiv \prod _{i=0}^{n-1}\left(b^{2^{i}}\right)^{a_{i}}\ ({\mbox{mod}}\ m)$

Pseudocódigo

Un ejemplo en pseudocódigo basado de Applied Cryptography por Bruce Schneier. Las entradas base, exponent, y modulus corresponden a b, e, y m en las ecuaciones dadas anteriormente.

function modular_pow(base, exponent, modulus) is
    if modulus = 1 then
        return 0
    Assert :: (modulus - 1) * (modulus - 1) does not overflow base
    result := 1
    base := base mod modulus
    while exponent > 0 do
        if (exponent mod 2 == 1) then
            result := (result * base) mod modulus
        exponent := exponent >> 1
        base := (base * base) mod modulus
    return result

Implementaciones

C

Algoritmo de derecha a izquierda

long powmod(long a, long e, long n){
  long accum = 1, apow, x;
  apow = a; x = e;
  while(x != 0){
	if (x & 0x01){ accum = (accum*apow) % n;};
	x >>= 1;
	apow = (apow * apow) % n;
  };
  return accum;
}

Algoritmo de derecha a izquierda con GMP

La biblioteca GMP tiene una función mpz_powm(result,a,e,n), así que no es necesaria esta implementación, pero se muestra aquí como método explicativo.

// High to Low powmod algorithm
int powmod(mpz_t result, mpz_t a,  mpz_t e,  mpz_t n) {
  // Use result as accum (temp variable)
  if (mpz_cmp_si(e,0) == 0) { // If exponent is zero
	mpz_set_ui(result, 1); // Set result to 1
	return 1;
  };
  mpz_set(result, a); // Set value of accum to a
 
  int bitptr = mpz_sizeinbase(e,2) - 1; // Find top bit in exponent
 
  for(bitptr--; bitptr >= 0; bitptr--) {
	mpz_mul(result,result,result); // result <-- result^2
	mpz_fdiv_r(result,result,n);   // result <-- result (mod n)
	if(mpz_tstbit(e,bitptr)) { // Is bit e[bitptr] == 1?
  	   mpz_mul(result,result,a);  	// result <-- result*a
  	   mpz_fdiv_r(result,result,n);   // result <-- result (mod n)
	};
  };
  return 1;
}

Java

Algoritmo de derecha a izquierda

public static long powmod(long a, long e, long n){
   long accum = 1;
   long x = e;
   long apow = a;
   while (x != 0){
      if ((x & 0x01) == 0x01){
         accum = (accum * apow) % n;
      };
      x >>= 1;
      apow = (apow * apow) % n;
   };
   return accum;
}

Algoritmo de derecha a izquierda con BigInteger

La clase BigInteger de Java tiene un método modPow(e, n), así que no es necesaria esta implementación, pero se muestra aquí como método explicativo.

// High to low powmod algorithm
public static BigInteger POWMOD(BigInteger a, BigInteger e, BigInteger n) {
    BigInteger accum = a;
    if (e.equals(BigInteger.ZERO)) {return BigInteger.ONE;};
    int bitptr = e.bitLength()-1;
    for(bitptr--; bitptr >= 0; bitptr--){
        accum = accum.multiply(accum).mod(n);
        if(e.testBit(bitptr)){
            accum = accum.multiply(a).mod(n);
        };
    };
    return accum;
}

Python

Algoritmo de derecha a izquierda

Python tiene una función pow(a, e, n), así que no es necesaria esta implementación, pero se muestra aquí como método explicativo.

def powmod(a,e,n):
    accum = 1; i = 0; apow2 = a
    while ((e>>i)>0):
        if ((e>>i) & 1):
            accum = (accum*apow2) % n
        apow2 = (apow2*apow2) % n
        i += 1
    return accum