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.

domingo, 11 de marzo de 2012

Subida a ritmo de M.C.M (2)

Relación con los factoriales

Dijimos en la entrada anterior que la sucesión A(n) subía rápido, pero la de factoriales más. Si dividimos los factoriales entre los MCM que estamos estudiando nos da esta otra sucesión C(n), que basta verla para comprender las distintas “velocidades”:

1, 1, 1, 2, 2, 12, 12, 48, 144, 1440, 1440, 17280, 241920, 3626800…(http://oeis.org/A025527)

¿Qué números no aportan nada y dejan los valores iguales? Los primos, porque para pasar de (n-1)! a n! hay que multiplicar por n, pero esta operación es la misma que hay que realizar en la sucesión A(n) si n es primo.

Podemos emprender el mismo estudio que para A(n) y es dividir cada término por  el anterior.



Tienes esta sucesión D(n) en http://oeis.org/A048671

En esta tabla vemos que las potencias de primos pr hacen crecer los términos en pr-1 y el resto aporta su propio valor. Para justificarlo volvemos a considerar el paso de (n-1)! a n! y de MCM(1,2,3,…n-1) a MCM(1,2,3,..n):

  •  Vimos que en los primos en ambos casos se multiplicaba por el mismo número primo y por eso en ellos C(N)/C(N-1)=1=p0=p1-1, luego se cumple.


  •  En el caso de las potencias de primos el factorial se incrementa multiplicándose por pr y los MCM s incrementan en p, luego el cociente se incrementará en pr-1, como hemos afirmado.


  •  En los demás casos el factorial se multiplica por n y el MCM queda igual, luego C(n) quedará también multiplicada por n.


Pero bueno, ¿qué es todo esto? Pues sencillamente, que B(n)*D(n)=n Era de esperar. Una aporta lo que le falta a la otra para ser n. Ahí lo tienes:

Los múltiplos y divisores nunca dejan de asombrarnos.

martes, 6 de marzo de 2012

Subida a ritmo de M.C.M (1)

Si te paras unos segundos, ¿sabrías descubrir cómo se ha generado esta sucesión?

1, 2, 6, 12, 60, 60, 420, 840, 2520, 2520, 27720, 27720, 360360, 360360, 360360, 720720, 12252240, 12252240, 232792560, 232792560, 232792560…(http://oeis.org/A003418)

Se parece a la de factoriales, pero crece a menos ritmo.

¿Ya lo sabes? Se trata del M.C.M. de los primeros números naturales: A(n)=MCM(1,2,3,…n). Así el 420=M.C.M(1,2,3,4,5,6,7)

La puedes engendrar con hoja de cálculo, escribiendo los primeros números y abajo encuentras el M.C.M del número de arriba y el número de la izquierda. No damos más detalles.


1
2
3
4
5
6
7
8
9
10
11
12
13
14
1
2
6
12
60
60
420
840
2520
2520
27720
27720
360360
360360













Una bonita pregunta es qué aporta cada número al resultado final del MCM. Observa en la tabla que el valor A(5)=60 y A(6)=60 también. ¿Por qué el 6 no ha aportado nada al cálculo? Parece ser que sus factores primos estaban ya contabilizados. Entonces, ¿cuáles aportan? Para verlo más claro dividiremos A(n) entre A(n-1)

Si dividimos cada MCM  por el anterior nos resulta la sucesión B(n), si definimos B(1)=1

1, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1, 13, 1, 1, 2, 17, 1, 19, 1, 1, 1, 23,…(http://oeis.org/A014963)

Con hoja de cálculo se ve mejor:


Sólo aportan un factor mayor que 1 los números primos y sus potencias. Es claro que es porque sólo ellos suponen algo nuevo. El resto, como el 12, usa factores que ya han aportado el 3 y el 4. ¿Qué ocurre entonces? Que al llegar a cada potencia de primo se habrá acumulado este tantas veces como indique esa potencia.

Estudia el 8. Antes de él ha aparecido el 2 como factor de sí mismo y como factor de 4. Con el 2 que aporta el 8 ya tenemos tres, que es precisamente el exponente correspondiente al 8.

En esta sucesión se van acumulando los factores primos de forma que al llegar sus potencias las reproducen exactamente.

Esto tiene una consecuencia muy elegante:


Si tomas todos los divisores de un número y multiplicas los factores que aportan al MCM (sucesión B(n)) nos resultará de nuevo ese número.

Por ejemplo, en el caso de 24 tendríamos:

Divisores:                 1,  2,  3,   4,   6,   8,   12,   24
Valores de B(n):     1,  2,  3,   2,   1,   2,   1,     1

Es evidente que el producto de los valores de B(n) vuelve a dar 24.

¿Conoces la función de Mangoldt? Si has leído a nuestro amigo Rafael Parra te sonará (http://hojamat.es/parra/funesp.pdf)

Pues bien, nuestra función B(n) es la exponencial de dicha función. Si tomas logaritmos en B(n) obtendrás 0, log(2), log(3), log(2), log(5), 0,… que es la definición de la función de Mahgoldt (tomamos la imagen de http://mathworld.wolfram.com/MangoldtFunction.html)


Quiere decir que si tomamos logaritmos en la fórmula de arriba nos resultará esta otra:


que podrás encontrar en textos de Teoría de Números. No seguimos por ahí.