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

domingo, 18 de marzo de 2012

Funciones recursivas en las hojas de cálculo

Cuando yo programaba hace años en Pascal se nos vendía su posibilidad de usar la recursión, es decir, que una función se llamara a sí misma, en declaraciones del tipo

Factorial(n)=n*factorial(n-1)

Esta y otras características nos hicieron abandonar el Basic como un lenguaje más primitivo y que no admitía funciones recursivas ni por asomo. Pasados bastantes años dejé la confección de programas ejecutables y consiguientemente el Pascal. Ahora que mis trabajos, por voluntad propia, los restrinjo a las hojas de cálculo y a su Basic, no uso la recursión…hasta hoy.

Preparando una próxima entrada se me ocurrió usar funciones recursivas en Excel, OpenOffice y LibreOffice (en Google Docs no funcionan las macros en Basic) con la sorpresa de que sí funcionaban bastante bien.

Toda función recursiva contiene una llamada a sí misma, directa o indirectamente a través de otra función. Como esto nos puede llevar a un proceso sin fin, debe contener también un código de parada, que suele ser una definición en un caso concreto, como veremos en los ejemplos.

La recursividad no se resuelve hasta que no desemboca en ese caso de parada. Mientras tanto hay que guardar los datos pendientes situándolos en una pila. Por tanto, ahí está el único problema de usar la recursividad en las hojas, y es que se puede agotar la pila si se alargan mucho los cálculos, con el consiguiente mensaje de error. Un fallo de principiante es programar una función recursiva sin facilitar su salida. En ese caso el error será más grave aún: un cálculo sin fin.

Explicamos a continuación algunas funciones recursivas, comenzando con el factorial, que es la más popular y que nos servirá para explicar algunos detalles:

Public Function factorial(n)
Dim f
If n = 0 Then f = 1 Else f = factorial(n - 1) * n
factorial = f
End Function

Es fácil entender  el código: Para evitar confusiones, comenzamos almacenando el factorial en la variable f, para al final recoger su valor en factorial. El cálculo de f es el clásico de la función n!: si n es cero, definimos el factorial como 1 y en caso contrario multiplicamos por n el factorial de n-1.

Prueba esta función en cualquiera de las tres hojas propuestas más arriba.

Este primer ejemplo contiene las dos partes imprescindibles en una función de este tipo:

- Alguna asignación de un valor concreto, que servirá para detener la primera fase de lectura de datos y comenzar los cálculos hacia atrás. Aquí es la asignación 0!=1

- La definición recursiva propiamente dicha, que, como es conocido, consiste en exigir que n!=n*(n-1)!

Lo explicamos con un esquema:






















Al intentar calcular el factorial de 7, el programa se encuentra con una referencia al factorial de 6, guarda el 7 en la pila y se dedica a calcular el nuevo factorial. Como no puede, almacena el 6 encima del 7 (es una pila) y lo intenta con el 5, y así va de fracaso en fracaso (flecha descendente) hasta llegar al 0.

El valor 0 admite el cálculo, porque está definido como 1. Resuelto esto, es como si el programa se preguntara: ¿por dónde iba? Acude a la pila y ve un 1, con lo que ya puede calcular 1!=1*0!=1 y así sigue (flecha ascendente) buscando datos en la pila y resolviendo los cálculos según la definición recursiva.

Es evidente que para números grandes la pila se puede agotar por falta de memoria asignada.

Con esta función recursiva (la más inútil que me he inventado nunca) puedes tener una idea de la amplitud de la pila de tu hoja de cálculo.

Public Function identidad(n)
Dim i
If n = 1 Then i = 1 Else i = identidad(n - 1) + 1
identidad = i
End Function


La he elegido para que no influya la magnitud de los números, sino la cantidad de ellos que permanecen en la pila. Con ella he llegado la valor n=3270 como el último que no me da error. En los siguientes no consigo realizar el cálculo.

¿Qué margen tendrá tu hoja de cálculo? Prueba a ver.


Ejemplos varios

Si deseas el enésimo número triangular, sólo tienes que usar este código:

Public Function triang(a)
Dim p
If a = 1 Then p = 1 Else p = triang(a - 1) + a
triang = p
End Function

También se entiende bien: en los números triangulares vamos añadiendo en cada paso una base del triángulo nueva con a elementos. Prueba esta función y si quieres compara los resultados con la clásica fórmula Tn=n(n+1)/2.

¿Puedes analizar esta función?

Public Function cuad(a)
Dim p
If a = 1 Then p = 1 Else p = cuad(a - 1) + 2 * a - 1
cuad = p
End Function

¿Por qué produce como resultado el cuadrado de a? Este es un bonito ejemplo de elevar un número al cuadrado sin multiplicar en ningún momento.

Y ya que estamos con números poligonales, podríamos generarlos todos con una función recursiva única que dependiera de a y también del número de lados del polígono. ¿Te atreves con ella?

¿Y qué opinas de esta, con dos variables?¿Qué resultado produce?

Public Function combi(m, n)
Dim c
If n = 0 Then c = 1 Else c = combi(m, n - 1) * (m - n + 1) / n
combi = c
End Function

Ahora un ejemplo más serio:

En una entrada anterior (http://hojaynumeros.blogspot.com/2012/02/suma-de-los-elementos-de-todos-los.html) descubrimos que la suma de todos los elementos de los subconjuntos de un conjunto de n elementos venía dada por la fórmula de recurrencia

    Sn=2Sn-1+ n*2n-1

y que da lugar a esta sucesión de valores en función de n

0, 1, 6, 24, 80, 240, 672, 1792, 4608, 11520, 28160, 67584, 159744, 372736, 860160, 1966080,…  (http://oeis.org/A001788)

Si definimos una función según esta recurrencia podremos reproducir esta lista en nuestra hoja de cálculo. Podría ser esta:

Public Function sumaelem(n)
Dim s
If n = 1 Then s = 1 Else s = 2 * sumaelem(n - 1) + n * 2 ^ (n - 1)
sumaelem = s
End Function

Con ella hemos construido esta tabla que coincide con la de OEIS



Un ejemplo elegante

Define esta función de texto

Public Function simetrico$(a$)
Dim s$
If a$ = "" Then s$ = "" Else s$ = simetrico(Right$(a$, Len(a$) - 1)) + Mid$(a$, 1, 1)
simetrico = s$
End Function

Escribe una palabra en una celda y aplícale esta función desde otra celda ¿Cuál es el resultado?

Como ves, todo esto es bastante divertido, pero no muy útil a causa del agotamiento del espacio de memoria asignado a la pila de datos.

Y ahora tú. ¿Cómo hallarías, mediante una función recursiva, el término general en estas sucesiones?

Progresiones aritméticas y geométricas.
La sucesión de Fibonacci (¡cómo no!)
La enésima potencia de un número dado.