Exponenciación rápida
En resumen: La exponenciación rápida calcula aⁿ con O(log n) multiplicaciones elevando al cuadrado y leyendo los dígitos binarios de n; tomar el resto tras cada paso da la exponenciación modular.
Problema
La clave de una conexión HTTPS
Cuando el navegador abre una página HTTPS, él y el servidor acuerdan una clave secreta compartida sin enviarla por la red. En la versión clásica (Diffie–Hellman) cada lado calcula una potencia modular: gᵃ mod p, donde a es un número secreto de decenas de dígitos. El mismo cálculo está en las firmas RSA y en el hash de un texto, que multiplica por potencias de una base.
La forma ingenua es un bucle que multiplica n veces:
def power_loop(a, n, m): result = 1 for _ in range(n): result = result * a % m return resultCon un exponente de 77 dígitos (por ejemplo 2²⁵⁵ − 19) son unas 6 · 10⁷⁶ multiplicaciones. A mil millones de multiplicaciones por segundo el bucle se ejecutaría durante 10⁶⁷ segundos, y el universo tiene unos 4 · 10¹⁷ segundos de edad. Sin embargo, el servidor hace miles de conexiones así por segundo.
La idea
El cuadrado de la mitad
3¹³ no necesita 13 multiplicaciones. 3¹³ = 3⁶ · 3⁶ · 3, y 3⁶ = 3³ · 3³. Cada paso divide el exponente por la mitad:
- aⁿ = (aⁿᐟ²)², si n es par;
- aⁿ = (aⁿᐟ²)² · a, si n es impar (n / 2 es división entera);
- a⁰ = 1: este es el caso base.
def power(a, n, m): if n == 0: return 1 half = power(a, n // 2, m) square = half * half % m if n % 2 == 1: return square * a % m return square
print(power(3, 13, 1_000_000_007))print(power(5, 1_000_000_000, 1_000_000_007))1594323142848001El resto mod m después de cada multiplicación mantiene los números pequeños: (x · y) mod m es lo mismo que ((x mod m) · (y mod m)) mod m. Aquí m = 1 000 000 007, un número primo y un límite habitual en los hashes y en los concursos de programación.
¿Cuánto cuesta? Cada llamada divide n entre 2, igual que la búsqueda binaria, así que la profundidad es de unos log₂ n: para n = 10⁹ son 30 marcos. La llamada no es de cola (elevar al cuadrado viene después), pero desde la primera lección sabes que 30 marcos no son nada para la pila. El trabajo es como máximo dos multiplicaciones por nivel: O(log n) en lugar de O(n).
La idea
Sin recursión: los dígitos binarios de n
Un bucle puede hacer lo mismo. 13 en binario es 1101: 13 = 8 + 4 + 1, así que 3¹³ = 3⁸ · 3⁴ · 3¹. El bucle lee los dígitos de n de derecha a izquierda (n % 2 es el último, y n /= 2 lo quita) y mantiene un a elevado al cuadrado: 3¹, 3², 3⁴, 3⁸… Con un dígito 1 el a actual entra en el resultado.
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return result
print(power(3, 13, 1_000_000_007))print(power(5, 1_000_000_000, 1_000_000_007))1594323142848001Las multiplicaciones son tantas como dígitos tiene n (una elevación al cuadrado por dígito) más los unos que contiene (una multiplicación del resultado cada uno). Esta es la versión que usan las bibliotecas: sin recursión, y sin depender del tamaño de la pila.
Míralo
Los cuadrados entran en el resultado
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return result
print(power(3, 13, 1_000_000_007)) # 1594323Paso a paso, en palabras
- result empieza en 1. n = 13, en binario 1101: cada dígito binario decide si a entra en el resultado.
- n = 13 es impar: result = 1 · 3 = 3. Después a pasa a ser 3² = 9 y n pasa a ser 6.
- n = 6 es par: result se queda en 3. Solo a se eleva al cuadrado: 9 · 9 = 81 = 3⁴, y n = 3.
- n = 3: result = 3 · 81 = 243 = 3⁵. a = 6561 = 3⁸, n = 1.
- n = 1: result = 243 · 6561 = 1594323 = 3¹³, porque 13 = 8 + 4 + 1. El bucle aún eleva a al cuadrado una vez más: 3¹⁶ es una multiplicación que nadie usa.
- n = 0 y el bucle se detiene: 7 multiplicaciones en lugar de 13. Con n = 10⁹ son 43 en lugar de mil millones.
Míralo
Contar las multiplicaciones
Las dos versiones con un contador, para exponentes de 10 a 10 000 000, y la rápida también para mil millones:
multiplications = 0
def power_loop(a, n, m): global multiplications result = 1 for _ in range(n): multiplications += 1 result = result * a % m return result
def power(a, n, m): global multiplications result = 1 a %= m while n > 0: if n % 2 == 1: multiplications += 1 result = result * a % m multiplications += 1 a = a * a % m n //= 2 return result
m = 1_000_000_007for n in [10, 1_000, 100_000, 10_000_000]: multiplications = 0 power_loop(5, n, m) slow = multiplications multiplications = 0 power(5, n, m) print(f"n = {n}: bucle {slow}, rápido {multiplications} multiplicaciones")multiplications = 0power(5, 1_000_000_000, m)print(f"n = 1000000000: rápido {multiplications} multiplicaciones")n = 10: bucle 10, rápido 6 multiplicacionesn = 1000: bucle 1000, rápido 16 multiplicacionesn = 100000: bucle 100000, rápido 23 multiplicacionesn = 10000000: bucle 10000000, rápido 32 multiplicacionesn = 1000000000: rápido 43 multiplicacionesEjecútalas también aquí. Bench cuenta con la misma regla (una multiplicación mod m es una operación), así que los números son los de la salida de arriba. No ejecutamos mil millones con el bucle ingenuo: Bench lo estima a partir de la curva.
- Bucle: n multiplicaciones
O(n) - Elevar al cuadrado
O(log n)
| Algoritmo | n = 10 | n = 1000 | n = 100.000 | n = 10.000.000 | n = 1.000.000.000 |
|---|---|---|---|---|---|
| Bucle: n multiplicaciones | 10 | 1000 | 100.000 | 10.000.000 | ≈ 1.000.000.000 |
| Elevar al cuadrado | 6 | 16 | 23 | 32 | 43 |
def power(a, n, m): result = 1 a %= m while n > 0: if n % 2 == 1: result = result * a % m a = a * a % m n //= 2 return resultEl bucle ingenuo crece como n, el rápido como el número de dígitos de n: una n diez veces mayor añade unos 3,3 dígitos binarios, es decir, de 3 a 7 multiplicaciones.
Pruébalo
15 o 16
¿Qué exponente necesita más multiplicaciones?
multiplications = 0
def power(a, n, m): global multiplications result = 1 a %= m while n > 0: if n % 2 == 1: multiplications += 1 result = result * a % m multiplications += 1 a = a * a % m n //= 2 return result
for n in [15, 16]: multiplications = 0 power(3, n, 1_000_000_007) print(f"{n}: {multiplications} multiplicaciones")15 = 1111: cuatro dígitos y cuatro unos, 8 multiplicaciones. 16 = 10000: cinco dígitos y un uno, 6 multiplicaciones. Un exponente mayor no significa más trabajo: lo que cuenta son los dígitos y los unos.
La idea
Cuando el producto no cabe
El resto mantiene los números por debajo de m, pero el producto de dos números así es de hasta m². Si m² no cabe en el tipo, el cálculo es incorrecto, y en silencio. Aquí está la misma potencia con un módulo mayor:
m = 10_000_000_019print("power:", power(5, 1_000_000_000, m))print("pow: ", pow(5, 1_000_000_000, m))power: 475867028pow: 475867028int crece según haga falta, así que el mismo código da la respuesta correcta también con un módulo grande. El precio es que cada multiplicación de un número grande es más lenta. El pow(a, n, m) incorporado hace la misma elevación al cuadrado, pero en C.La idea
Lo incorporado, y un exponente de 77 dígitos
En código real una potencia modular no se escribe a mano:
secret = 2**255 - 19p = 2**127 - 1print("Dígitos en el exponente: " + str(len(str(secret))))print("Multiplicaciones al elevar al cuadrado: " + str(secret.bit_length() + secret.bit_count()))print("5^secret mod p =", pow(5, secret, p))Dígitos en el exponente: 77Multiplicaciones al elevar al cuadrado: 5085^secret mod p = 97285997395055909576908201857106207893pow(a, n, m) con tres argumentos eleva al cuadrado y toma el resto por el camino. a ** n % m es correcto, pero primero construye todo aⁿ: con un exponente de 77 dígitos no terminaría nunca. 508 multiplicaciones de números de 39 dígitos en lugar de 6 · 10⁷⁶: por eso la clave se acuerda en milisegundos. Las bibliotecas necesitan aún menos, porque miran varios dígitos a la vez, pero la idea es la misma.
Pruébalo
Encuentra el error
Esta versión devuelve el resultado correcto. Pero para n = 1 000 000 hace 2.048.575 multiplicaciones en lugar de 27. ¿Qué línea tiene la culpa?
Elige la línea con el error.
La mitad se calcula dos veces. Cada nivel duplica las llamadas: 2, 4, 8… y abajo del todo hay unas n. La profundidad sigue siendo log₂ n, así que la pila no nota nada, pero el trabajo vuelve a ser O(n). La solución es poner la mitad en una variable y luego half * half.
Pruébalo
Completa
¿Cuándo entra el a actual en el resultado?
def power(a, n, m):
result = 1
a %= m
while n > 0:
if :
result = result * a % m
a = a * a % m
n //= 2
return resultCuando el último dígito binario de n es 1, es decir, cuando n es impar.
Ejercicio
Intereses al céntimo
Resuélvelo en tu editor.
A cuánto crece el depósito. Pones 1000 € al 1 % de interés mensual, que se suma al importe. Lee el número de meses e imprime el importe redondeado al céntimo. Calcúlalo de forma exacta, con enteros: 1000 € son 100 000 céntimos, y tras n meses el importe es 100 000 · 101ⁿ / 100ⁿ céntimos. Escribe tu propia exponenciación rápida con int, sin módulo, y redondea medio céntimo hacia arriba.
120Tras 120 meses: 3300,39 €Con 120 meses las multiplicaciones son pocas incluso con el bucle ingenuo. Aquí, sin embargo, los números crecen (101¹²⁰ tiene 241 dígitos), así que cada multiplicación siguiente cuesta más, y ahorrar multiplicaciones importa aún más.
Más tareas
- Interés diario durante 30 años. Lo mismo con un 0,01 % diario durante 10 950 días. Cuenta las multiplicaciones de las dos versiones.
- Los últimos dígitos. Imprime los últimos 6 dígitos de 2¹⁰⁰⁰⁰⁰⁰⁰⁰⁰ (módulo 1 000 000), sin construir el número entero.
Más precisamente
Más sobre potencias
Más precisamente, para curiosos
Números de Fibonacci en O(log n). La misma elevación al cuadrado funciona con una matriz de 2 × 2 en lugar de un número: el n-ésimo número de Fibonacci sale de la n-ésima potencia de [[1, 1], [1, 0]], con unas 2 log₂ n multiplicaciones de matrices.
Curvas elípticas. Hoy TLS suele acordar la clave con una curva elíptica (X25519). Allí, en lugar de la multiplicación hay «sumar puntos», y en lugar de elevar al cuadrado hay «duplicar», pero el algoritmo es el mismo: duplicar y sumar sobre los dígitos binarios del número secreto.
Tiempo constante. El bucle de la lección multiplica el resultado solo con un dígito 1. Quien mide el tiempo puede adivinar cuántos unos tiene el número secreto. Por eso las bibliotecas criptográficas hacen el mismo trabajo para cada dígito (la escalera de Montgomery), y por eso la criptografía no se escribe a mano, sino que se toma de una biblioteca.
- Python:
pow(a, -1, m)(desde 3.8) da el inverso de a mod m.a ** ncon enteros es exacto, pero crece sin límite.
Test