viernes, 23 de marzo de 2012

De la multiplicación rusa a la exponenciación modular

(Con esta entrada participamos en el 3.14 Carnaval de Matemáticas, cuyo anfitrión es el blog Hablando de la Ciencia)


En algunos problemas, como en el criterio de primalidad de Fermat, que veremos en una próxima entrada, debemos elevar un elemento de Zm a un exponente grande. Son frecuentes los problemas del tipo “¿en qué cifra termina 263721?” o “¿es cierta la congruencia 234125º1 (mod 23)?


Si el exponente es grande pueden desembocar en cálculos muy complicados, por lo que se acude a la exponenciación por duplicación. Estas técnicas que se basan en duplicar son muy antiguas. Ya conocemos el método usado en Egipto (ver http://hojamat.es/parra/mat_antig.pdf) y posteriormente la llamada multiplicación a la rusa.

Tenemos implementada, como una curiosidad, esta multiplicación para hojas de cálculo (ver http://hojamat.es/sindecimales/aritmetica/herramientas/herrarit.htm#peque) y también tienes una exposición teórica en http://tiopetrus.blogia.com/2005/042501-multiplicacion-a-la-rusa-1-.php

Aquí nos va a interesar la parte común de los algoritmos de duplicación.

Recordemos el algoritmo de la multiplicación rusa:
















Si multiplicamos, por ejemplo, 87 por 456, vamos dividiendo 87 entre dos de forma entera, sin decimales, hasta llegar al 1. Simultáneamente duplicamos el otro factor 456 cada vez que dividamos el otro. Después sumamos los múltiplos de este número que se correspondan con los cocientes impares del otro, en el ejemplo 456+912+1824+7296+29184=39672 que coincide con el producto de 87*456.

Esto funciona porque los cocientes impares producen un 1 en la representación binaria de 87 y los pares un cero, por lo que tiene sentido sumar sólo los primeros, y como estamos duplicando en cada proceso, lo que hemos conseguido es lo siguiente:

87*456=(1+2+4+16+64)*456=456+912+1824+7296+29184=39672

En forma de algoritmo podría expresarse así:

Public Function rusa(a,b)
Dim  s
s = 0      ‘Se inicia la suma a cero
While a > 0  ‘Mientras a no llegue a cero, se divide entre 2
If (a / 2) <> Int(a / 2) Then s = s + b  ‘Si es impar se suma
b = 2 * b   ‘Se duplica b
a = Int(a / 2)  ‘ Se divide a
Wend
rusa = s  ‘ La función recoge el valor de s
End Function

Si copias este listado, lo puedes trasladar al módulo Basic de una hoja de cálculo para comprobarlo. Los parámetros a y b son los factores.

Podíamos intentar un proceso similar con la potenciación. Por ejemplo, para calcular 713 podríamos proceder así:
















En lugar de duplicar, elevamos al cuadrado, y al final multiplicamos en vez de sumar.

Para el conjunto Z la potenciación desemboca pronto en números grandes, pero no así en Zm, pues los resultados siempre tendrán como cota m y este método puede ser muy útil. Incluso se puede intentar mentalmente. Por ejemplo, calcular 763 (mod 5): 7º2 (mod 5); 72º4 (mod 5); 74º1 (mod 5); 78º1 (mod 5) y ya todos valen 1, luego 763=732+16+8+4+2+1º2*4*1*1*1*1º8º3, luego 763º3 (mod 5)

El algoritmo de la multiplicación rusa se puede adaptar fácilmente a esta exponenciación

Public Function expo(a, e, m) ‘Los parámetros son base, exponente y módulo
Dim p
p = 1 ‘inicia el producto final
While e > 0
If (e / 2) <> Int(e / 2) Then p = p * a: p = p - m * Int(p / m) ‘si es impar multiplica y reduce a resto módulo m
a = a ^ 2: a = a - Int(a / m) * m ‘se eleva al cuadrado reduciendo a módulo m
e = Int(e / 2) ‘se divide el exponente
Wend
expo = p
End Function

Esta sería la exponenciación modular o binaria, que resulta imprescindible en cálculos con grandes números, porque sólo utiliza un número de multiplicaciones del orden del logaritmo de n, y no de n como el algoritmo clásico.

No resisto incluir la versión recursiva que aparece en Wikipedia (para Z)


Si has leído nuestra anterior entrada, sabrás que, con limitaciones, el Basic de las hojas admite recursividad. En efecto, hemos probado esta definición para módulo m y funciona:

Public Function expo2(a, e, m)
Dim ep
If e = 1 Then
ep = a  ‘Definición de parada de la recursión
ElseIf e = Int(e / 2) * 2 Then ‘Caso par
ep = (expo2(a, e / 2, m)) ^ 2: ep = ep - m * Int(ep / m) ‘Elevación al cuadrado
Else
ep = a * expo2(a, e - 1, m): ep = ep - m * Int(ep / m) ‘Caso impar
End If
expo2 = ep
End Function