Ada Lovelace: Pionera de la programación

Esta entrada del blog de EduPython fue elaborada por mi buena amiga y ex-alumna Samantha Montserrat Ponce Aparicio, quien se tomó el tiempo para leer y escribir sus impresiones del libro titulado “Ada, the Enchantress of Numbers: Poetical Science” de Betty Alexandra Toole, publicado por Critical Connection en octubre del 2010.

Samantha es una de las mujeres programadoras más talentosas que he conocido. Actualmente estudia la carrera de Ingeniero en Sistemas Computacionales en el Tecnológico de Monterrey, Campus Estado de México, y está por graduarse en unos pocos meses.

Es un gran gusto y honor tener a Samantha como la primera autora invitada del blog de EduPython. 

Ariel Ortiz.
Agosto, 2016.

---

Antes de leer sobre la vida de Ada Lovelace, lo único que sabía es que era considerada la primera programadora. Nunca se me ocurrió averiguar qué programó, en qué computadora, cómo era su vida y si realmente había sido la primera persona en programar.

Los primeros años de vida de Ada

La vida de Ada Lovelace fue verdaderamente complicada. Fue hija de Lord Byron, sí, el poeta romántico, y Lady Byron quienes se separaron unos meses después de su nacimiento por las diferencias de pensamiento que tenían, ella pensaba en números y él en versos.

Ada de niña, retrato en exhibición
en Somerville College, Oxford.

Augusta Ada Byron nació en Londres el 10 de diciembre de 1815. Después de la separación, Lady Byron se quedó con la custodia de Ada y decidió educarla promoviendo en ella el gusto por la ciencia y las matemáticas, con la esperanza de que no fuera a seguir los malos pasos de su padre en el mundo del arte y la poesía. Ada aprendió a concretar ideas en forma de bloques en lugar de simplemente seguir la teoría, a pesar de que en ese tiempo solamente se utilizaban conceptos intangibles, incluso para enseñar a los niños.

El poeta George Gordon Byron, padre de Ada.

Uno de los datos tristes sobre Ada es que nunca conoció a su padre. Él murió en abril de 1824 y ella solamente lo vio en un retrato poco antes de cumplir 20 años. Por eso resintió a su madre durante toda su vida.

Al no tener hermanos, decidió enfocarse en observar a su gato Mrs. Puff, porque todos sabemos que la vida con gatos es mejor. A los 12 años decidió que quería volar, así que después de investigar e imaginar, escribió todo lo que encontró en un libro llamado Flyology, pero como su madre se burló de ella, abandonó el proyecto, a pesar de que sus teorías eran bastante acertadas. William Henson utilizó las mismas ideas para el diseño de su Aerial Steam Carriage, también llamado Ariel, como un profesor que conozco ☺.

Sus tutores se quejaban porque ella cuestionaba todo lo que enseñaban. Siempre buscó aplicar conceptos matemáticos en juegos, diagramas y metáforas, pues esa era la forma en la que ella aprendía mejor.

Después de haber estado enferma y en cama por tres años, decidió que además de estudiar matemáticas, quería aprender música y equitación. Cuando leí sobre esto me sentí identificada con ella, porque toda la vida me han gustado las matemáticas y las computadoras, pero también me gusta cantar, bailar y probar artes nuevas, por ejemplo, hace pocos meses empecé a tomar clases de danza aérea en telas.

La máquina analítica

Cuando Ada tenía casi 18 años conoció a Charles Babbage, y la vida de ambos cambió. Charles Babbage es conocido por haber inventado la primera computadora mecánica, que se programaba por medio de tarjetas perforadas. Su único y gran problema es que el gobierno no le daba los fondos necesarios para sus proyectos.

Charles Babbage, matemático,
filósofo, ingeniero e inventor.

Es ahí cuando Ada se ofreció a hacer notas sobre su invento, la máquina analítica, para conseguir fondos para toda la investigación de Babbage. En esas notas, Ada incluía diagramas, tablas y párrafos enteros sobre el poder de la máquina analítica. Por primera vez se utilizaron conceptos tales como: ciclo, índice y paralelismo, los cuales son básicos en la programación hoy en día.

Máquina analítica de Charles Babbage.
Ejemplar en exhibición en el museo de ciencia de
Londres. Solo una parte de esta máquina
se construyó antes de su muerte en 1871.

Como era muy importante destacar sobre invenciones anteriores, incluyendo la máquina diferencial diseñada también por Babbage, Ada se dio a la tarea de buscar un problema que solo pudiera resolverse con la máquina analítica. El problema que finalmente escogió y programó fue el de los números de Bernoulli, por esto se dice que ella fue la primera programadora. Lo cierto es que Babbage fue el primer programador, ya que tenía algunos ejemplos para poder probar su máquina, pero Ada fue quien documentó un programa computacional completo por primera vez.

El legado de Ada

Lo que más me encanta de Ada es su forma de pensar. Siempre cuestionaba todo y no se quedaba callada si no entendía algún concepto o no le parecía correcto. Quería aprender los principios y saber porqué las cosas eran como decían ser. Su pensamiento artístico le ayudaba a ver las matemáticas desde otro punto de vista y eso fue su gran colaboración con el trabajo de Babbage. Estoy segura de que sin las notas de Ada, la máquina analítica de Babbage hubiera pasado desapercibida, ya que ni el mismo Babbage sabía el poder de su invento hasta que Ada comenzó a ayudarlo.

Ada era, a mi parecer, bastante popular, o por lo menos se juntaba con personas muy reconocidas. Su tutor por un tiempo fue Augustus De Morgan, quien hizo las leyes de De Morgan, su mejor amigo fue Charles Babbage, sus amigos fueron Sir David Brewster, el inventor del caleidoscopio, Andrew Crosse, conocido por la electrocristalización, Charles Dickens, Michael Faraday y Florence Nightingale.

Se casó con William King en 1835 y en 1838 se volvieron conde y condesa de Lovelace, de ahí su nombre. Durante toda su vida se enfermó bastante y terminó muriendo en 1852 a la temprana edad de 37 años.

Daguerrotipo de Augusta Ada King,
Condesa de Lovelace, 1844.

Me parece importante que las personas en el mundo de la programación conozcan sobre Ada Lovelace, en especial las mujeres, para que se den cuenta que una mentalidad diferente puede lograr grandes cosas. Ada, a pesar de haber tenido una vida corta, hizo una contribución tan grande que incluso existe un importante lenguaje de programación que lleva su nombre.

Es sorprendente lo que se puede lograr con una forma diferente de pensar, puede que con eso logremos cosas aún más grandes, superando incluso a la invención de la computadora.

Apéndice: Los números de Bernoulli

El siguiente diagrama podría ser considerado el primer programa computacional publicado. Describe la manera de calcular los números de Bernoulli utilizando la máquina analítica de Charles Babbage:

Diagrama para calcular los números de Bernoulli.

Los números de Bernoulli se utilizan en algunas series de expansión para diferentes funciones (trigonométricas, hiperbólicas, gamma, etc.) y son extremadamente importantes en el análisis y la teoría de números.

De manera simplificada, para evitar meternos en demasiados detalles, el k-ésimo número de Bernoulli \(B_k\) se define usando la siguiente relación de recurrencia: \begin{equation} B_k = \left\{\begin{matrix} 1 & \textrm{si} \;\; k = 0\\ - \sum\limits_{i = 0}^{k - 1} \binom{k}{i} \frac{B_i}{k + 1 - i} & \textrm{si} \;\; k > 0 \end{matrix}\right. \end{equation} La fórmula anterior utiliza el coeficiente binomial, el cual se define de esta manera: \begin{equation*} \binom{k}{i} = \frac{k!}{i! \;(k-i)!} \end{equation*} El siguiente código es la traducción directa a Python 3 de estas dos fórmulas:

from math import factorial

def binomial(k, i):
    """Regresa el coeficiente binomial C(k, i),
    es decir, el número de formas distintas de
    seleccionar i elementos a partir de un
    conjunto de tamaño k.
    """
    return (factorial(k) //
            (factorial(i) * factorial(k - i)))

def bernoulli(k):
    """Regresa el k-ésimo elemento de la
    secuencia de números de Beroulli.
    """
    if k == 0:
        return 1
    else:
        return -sum([binomial(k, i) * bernoulli(i)
                     / (k + 1 - i)
                     for i in range(k)])

Con esto podemos usar el siguiente ciclo for para imprimir los números de Bernoulli B₀ al B₂₀:

for i in range(21):
    print('B({:2}) = {:10.5f}'.format(i, bernoulli(i)))

El uso del método format() permite que la salida aparezca adecuadamente alineada y con exactamente cinco dígitos después del punto decimal. La salida es la siguiente:

B( 0) =    1.00000
B( 1) =   -0.50000
B( 2) =    0.16667
B( 3) =   -0.00000
B( 4) =   -0.03333
B( 5) =   -0.00000
B( 6) =    0.02381
B( 7) =    0.00000
B( 8) =   -0.03333
B( 9) =   -0.00000
B(10) =    0.07576
B(11) =   -0.00000
B(12) =   -0.25311
B(13) =    0.00000
B(14) =    1.16667
B(15) =   -0.00000
B(16) =   -7.09216
B(17) =    0.00000
B(18) =   54.97118
B(19) =   -0.00000
B(20) = -529.12424

Vale la pena notar que B_k = 0 para todos los valores de k impares, con la excepción de k = 1.

Es conveniente comentar también que nuestra implementación es extremadamente ineficiente. Entre más grande sea el valor de k el programa se vuelve perceptiblemente más lento. Alternativamente, el algoritmo de Akiyama–Tanigawa permite calcular los números de Bernoulli de manera optimizada. El sitio de RosettaCode.org muestra cómo implementar dicho algoritmo usando Python y otra treintena de lenguajes de programación.

Combinaciones, permutaciones y otras diversiones

Imaginemos que estamos en una heladería y queremos un Tres Marías. Este delicioso postre se prepara con tres bolas de helado de diferentes sabores. La heladería cuenta con helado de los siguientes sabores: cereza, chocolate, fresa, nuez y vainilla. Con esta información, ¿qué elecciones de sabores tenemos para crear nuestro Tres Marías?      

La respuesta son las siguientes diez opciones:

• Cereza, chocolate y fresa.
• Cereza, chocolate y nuez.
• Cereza, chocolate y vainilla.
• Cereza, fresa y nuez.
• Cereza, fresa y vainilla.
• Cereza, nuez y vainilla.
• Chocolate, fresa y nuez.
• Chocolate, fresa y vainilla.
• Chocolate, nuez y vainilla.
• Fresa, nuez y vainilla.

Analicemos otro problema. Tenemos las siguientes tres letras: i, o, r. ¿Qué palabras podemos formar usando solamente esas tres letras?  Supongamos que no nos importa que algunas de las palabras formadas no existan en el idioma español. La respuesta son seis palabras:

• ior
• iro
• oir
• ori
• rio
• roi

El problema de los helados es un problema de combinaciones: deseamos determinar las formas de agrupar los elementos de un conjunto en donde no importa el orden en que se colocan dichos elementos. Un Tres Marías de chocolate, fresa y vainilla es igual a uno de fresa, vainilla y chocolate. Por otro lado, el problema de las letras es un problema de permutaciones: a diferencia de las combinaciones, aquí el orden sí nos interesa. No es lo mismo “rio” que “oir”.

Tradicionalmente, muchos libros de probabilidad y estadística comienzan con una descripción de cómo calcular el número de combinaciones o permutaciones de un conjunto. Sin embargo, si lo que nos interesa es obtener un listado con dichas combinaciones o permutaciones es muy probable que estos libros no expliquen la manera de hacerlo. Afortunadamente, con Python y un poco de ingenio podemos resolver este problema.

NOTA: Todo el código que aquí se presenta fue probado con Python 3.5.

Recordando al conjunto potencia

En nuestra entrada anterior discutimos cómo programar una función que calcula el conjunto potencia. Repetiremos aquí el código correspondiente por conveniencia, ya que estaremos haciendo uso de él más adelante.

def potencia(c):
    """Calcula y devuelve el conjunto potencia del 
       conjunto c.
    """
    if len(c) == 0:
        return [[]]
    r = potencia(c[:-1])
    return r + [s + [c[-1]] for s in r]

def imprime_ordenado(c):
    """Imprime en la salida estándar todos los
       subconjuntos del conjunto c (una lista de
       listas) ordenados primero por tamaño y
       luego lexicográficamente. Cada subconjunto
       se imprime en su propia línea. Los
       elementos de los subconjuntos deben ser
       comparables entre sí, de otra forma puede
       ocurrir un TypeError.
    """
    for e in sorted(c, key=lambda s: (len(s), s)):
        print(e)

Combinaciones

Se llama combinaciones de m elementos tomando n elementos a la vez (donde m ≥ n) a todas las agrupaciones posibles de tamaño n que pueden hacerse con los m elementos. Recordemos que no importa el orden de los elementos. Para esta discusión también supondremos que los elementos no se repiten.

Regresemos al problema original de las tres bolas de helado. Lo que queremos es encontrar todas las combinaciones que se pueden formar a partir de 5 sabores de helado tomando 3 sabores a la vez.

Primero veamos qué obtenemos cuando calculamos el conjunto potencia con los cinco sabores de helado:

>>> imprime_ordenado(
...     potencia(['cereza', 'chocolate', 'fresa', 
...               'nuez', 'vainilla']))
[]
['cereza']
['chocolate']
['fresa']
['nuez']
['vainilla']
['cereza', 'chocolate']
['cereza', 'fresa']
['cereza', 'nuez']
['cereza', 'vainilla']
['chocolate', 'fresa']
['chocolate', 'nuez']
['chocolate', 'vainilla']
['fresa', 'nuez']
['fresa', 'vainilla']
['nuez', 'vainilla']
['cereza', 'chocolate', 'fresa']
['cereza', 'chocolate', 'nuez']
['cereza', 'chocolate', 'vainilla']
['cereza', 'fresa', 'nuez']
['cereza', 'fresa', 'vainilla']
['cereza', 'nuez', 'vainilla']
['chocolate', 'fresa', 'nuez']
['chocolate', 'fresa', 'vainilla']
['chocolate', 'nuez', 'vainilla']
['fresa', 'nuez', 'vainilla']
['cereza', 'chocolate', 'fresa', 'nuez']
['cereza', 'chocolate', 'fresa', 'vainilla']
['cereza', 'chocolate', 'nuez', 'vainilla']
['cereza', 'fresa', 'nuez', 'vainilla']
['chocolate', 'fresa', 'nuez', 'vainilla']
['cereza', 'chocolate', 'fresa', 'nuez', 'vainilla']

El resultado consta de 32 subconjuntos (2⁵=32). Si observamos cuidadosamente, hay exactamente 10 subconjuntos que son de cardinalidad 3. Esos 10 subconjuntos son precisamente las 10 diez combinaciones que se pueden formar a partir de un total de 5 elementos tomando 3 elementos a la vez. Sabiendo esto, podemos definir en una sola línea la función en Python que sirve para calcular combinaciones:

def combinaciones(c, n):
    """Calcula y devuelve una lista con todas las
       combinaciones posibles que se pueden hacer
       con los elementos contenidos en c tomando n
       elementos a la vez.
    """
    return [s for s in potencia(c) if len(s) == n]

Como podemos observar, aquí utilizamos una lista por comprensión, la cual se puede leer así: para cada subconjunto s que pertenece al resultado del conjunto potencia de c, conservar solo aquellos valores de s en donde su cardinalidad sea igual a n. Probando la función podemos ver que sí produce los resultados esperados:

>>> imprime_ordenado(
        combinaciones(['cereza', 'chocolate', 'fresa',
                       'nuez', 'vainilla'], 3))
['cereza', 'chocolate', 'fresa']
['cereza', 'chocolate', 'nuez']
['cereza', 'chocolate', 'vainilla']
['cereza', 'fresa', 'nuez']
['cereza', 'fresa', 'vainilla']
['cereza', 'nuez', 'vainilla']
['chocolate', 'fresa', 'nuez']
['chocolate', 'fresa', 'vainilla']
['chocolate', 'nuez', 'vainilla']
['fresa', 'nuez', 'vainilla']

El número de posibles combinaciones se denota \({}_m C_n \) y se calcula así: $$ {}_m C_n = \frac{m!}{n! \; (m - n)!} $$ En donde m es la cardinalidad del conjunto inicial y n es la cardinalidad de los subconjuntos que deseamos formar. El signo de admiración es la operación de factorial. Verificando nuestro ejemplo, donde m = 5 y n = 3, el resultado obtenido es el que hemos estado manejando: 5!/(3!*2!) = 120/(6*2) = 120/12 = 10.

La función en Python que calcula el número de combinaciones queda así:

from math import factorial

def numero_combinaciones(m, n):
    """Calcula y devuelve el número de combinaciones
       posibles que se pueden hacer con m elementos
       tomando n elementos a la vez.
    """
    return factorial(m) // (factorial(n) * factorial(m - n))

Podemos asegurarnos que las dos funciones que definimos, combinaciones() y numero_combinaciones(), son consistentes entre sí:

>>> len(combinaciones(range(20), 1))
20
>>> len(combinaciones(range(20), 1)) \
    == numero_combinaciones(20, 1)
True
>>> len(combinaciones(range(20), 10))
184756
>>> len(combinaciones(range(20), 10)) \
    == numero_combinaciones(20, 10)
True
>>> len(combinaciones(range(20), 20))
1
>>> len(combinaciones(range(20), 20)) \
    == numero_combinaciones(20, 20)
True

NOTA: El carácter de diagonal invertida (\) usado arriba indica que la instrucción continúa en la línea de abajo.

Permutaciones

Una permutación es la variación de la disposición u orden de los elementos de un conjunto.  Usaremos recursión para diseñar un algoritmo que permita permutar una lista. Hay que definir, entonces, dos cosas: el caso base y la llamada recursiva.

Caso base: El resultado de permutar un conjunto vacío es un conjunto que contiene al conjunto vacío.

Llamada recursiva: Nos interesa resolver una versión más simple del problema. Puede ser algo así: permutar el mismo conjunto de entrada pero quitándole un elemento, por ejemplo el primero. Analicemos la siguiente figura para comprender lo que deseamos lograr:

La llamada original es: permuta({x, y, z}). Vamos a hacer una llamada recursiva con el mismo conjunto pero sin su primer elemento. La llamada sería: permuta({y, z}). Por definición (o sea, por acto de fe) la llamada recursiva debe devolvernos lo siguiente: {{y, z}, {z, y}}. Ahora respondamos a la pregunta: ¿Qué se le debe hacer al resultado de la llamada recursiva permuta({y, z}) para convertirlo en el resultado esperado de la llamada original permuta({x, y, z})? El uso de los colores en la figura de arriba nos permiten responder a esta pregunta.

Tomemos el primer subconjunto del resultado de la llamada recursiva: {y, z}. A partir de éste, debemos insertar x (el elemento eliminado al momento de hacer la llamada recursiva) en cada posición posible de dicho subconjunto. Dado que hay dos elementos en el subconjunto, existen tres posibles posiciones de inserción: al inicio, en medio y al final. Los tres nuevos subconjuntos serían: {x, y, z}, {y, x, z} y {y, z, x}. Lo anterior hay que repetirlo para el segundo subconjunto de la llamada recursiva, en esta caso: {z, y}. Ahora obtendríamos: {x, z, y}, {z, x, y} y {z, y, x}. El proceso se repetiría de esta misma forma en caso de que el resultado de la llamada recursiva tuviera más subconjuntos. Todos los subconjuntos que generamos aquí conforman el resultado esperado de la llamada original.

Dividamos el problema en varias funciones de Python para simplificar su resolución.

La siguiente función crea una nueva lista a partir de una lista existente lst pero insertando el elemento x en la posición del índice i:

def inserta(x, lst, i):
    """Devuelve una nueva lista resultado de insertar
       x dentro de lst en la posición i.
    """
    return lst[:i] + [x] + lst[i:]

La expresión lst[:i] es una rebanada de la lista lst comenzando desde el inicio y hasta el elemento anterior al índice i. De manera similar, la expresión lst[i:] es una rebanada de la lista lst comenzando en el elemento del índice i y hasta el final de la lista. En medio de esas dos rebanadas creamos una nueva lista conteniendo a x. Por último, concatenamos todo en orden para crear la lista resultante.

NOTA:  Las listas de Python tienen un método insert() que podría parecerse mucho a nuestra función. Sin embargo nuestra versión no modifica la lista original, mientras que la de Python sí. Las funciones que definimos más adelante suponen que las inserciones no alteran la lista original.

Ejemplos de uso:

>>> inserta(42, [1, 2, 3], 0)
[42, 1, 2, 3]
>>> inserta(42, [1, 2, 3], 3)
[1, 2, 3, 42]
>>> inserta(42, [1, 2, 3], 1)
[1, 42, 2, 3]

La siguiente función nos genera una lista de listas con todas las posibles inserciones de un elemento:

def inserta_multiple(x, lst):
    """Devuelve una lista con el resultado de
       insertar x en todas las posiciones de lst.  
    """
    return [inserta(x, lst, i) for i in range(len(lst) + 1)]

La función range() produce los índices de todas las posiciones en las que deseamos hacer una inserción (desde 0 hasta el número de elementos de la lista original) y la lista por comprensión se encarga de llamar nuestra función inserta() para hacer el resto del trabajo.

La función inserta_multiple() puede usarse así:

>>> inserta_multiple(42, [])
[[42]]
>>> inserta_multiple(42, [1])
[[42, 1], [1, 42]]
>>> inserta_multiple(42, [1, 2])
[[42, 1, 2], [1, 42, 2], [1, 2, 42]]
>>> inserta_multiple(42, [1, 2, 3])
[[42, 1, 2, 3], [1, 42, 2, 3], [1, 2, 42, 3], [1, 2, 3, 42]]

Ahora sí, ya estamos en condiciones de implementar la función permuta():

def permuta(c):
    """Calcula y devuelve una lista con todas las
       permutaciones posibles que se pueden hacer
       con los elementos contenidos en c.
    """
    if len(c) == 0:
        return [[]]
    return sum([inserta_multiple(c[0], s)
                for s in permuta(c[1:])],
               [])

La expresión permuta(c[1:]) genera recursivamente todas las permutaciones de nuestro conjunto original c pero sin su primer elemento. Por cada subconjunto s devuelto por permuta(), la lista por comprensión realiza el inserta_multiple() con el primer elemento de c como argumento. Hay que tener en cuenta que esta lista por comprensión devuelve una lista de listas de listas. Esto quiere decir que hay un nivel en exceso de listas anidadas. Usamos la función sum() para concatenar las listas que están anidadas en la lista más externa y con ello nos queda el resultado que deseamos. El segundo parámetro de sum() es el valor con el que inicia la sumatoria, que por omisión es el número cero pues sum() se usa principalmente para sumar números. Sin embargo, si en su lugar se le manda una lista vacía como argumento la función sum() realiza una concatenación de listas, que es justamente lo que queremos.

Probemos permuta() con el segundo problema que propusimos al inicio de esta entrada:

¿Qué palabras podemos formar usando las letras i, o, r?

>>> imprime_ordenado(permuta(['i', 'o', 'r']))
['i', 'o', 'r']
['i', 'r', 'o']
['o', 'i', 'r']
['o', 'r', 'i']
['r', 'i', 'o']
['r', 'o', 'i']

Probando con otras entradas:

>>> permuta([])
[[]]
>>> permuta([1])
[[1]]
>>> permuta([1, 2])
[[1, 2], [2, 1]]
>>> imprime_ordenado(permuta([1, 2, 3]))
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]
>>> imprime_ordenado(permuta([1, 2, 3, 4]))
[1, 2, 3, 4]
[1, 2, 4, 3]
[1, 3, 2, 4]
[1, 3, 4, 2]
[1, 4, 2, 3]
[1, 4, 3, 2]
[2, 1, 3, 4]
[2, 1, 4, 3]
[2, 3, 1, 4]
[2, 3, 4, 1]
[2, 4, 1, 3]
[2, 4, 3, 1]
[3, 1, 2, 4]
[3, 1, 4, 2]
[3, 2, 1, 4]
[3, 2, 4, 1]
[3, 4, 1, 2]
[3, 4, 2, 1]
[4, 1, 2, 3]
[4, 1, 3, 2]
[4, 2, 1, 3]
[4, 2, 3, 1]
[4, 3, 1, 2]
[4, 3, 2, 1]

El número de permutaciones que se pueden obtener para un conjunto de cardinalidad n es: n!. Esto lo podemos verificar con el siguiente código:

>>> from math import factorial
>>> len(permuta(range(3)))
6
>>> len(permuta(range(3))) == factorial(3)
True
>>> len(permuta([]))
1
>>> len(permuta([])) == factorial(0)
True
>>> len(permuta(range(1)))
1
>>> len(permuta(range(1))) == factorial(1)
True
>>> len(permuta(range(5)))
120
>>> len(permuta(range(5))) == factorial(5)
True
>>> len(permuta(range(8)))
40320
>>> len(permuta(range(8))) == factorial(8)
True

En ocasiones puede resultar deseable calcular las permutaciones de un conjunto pero tomando solamente n elementos a la vez. Para ello podemos recurrir a la función combinaciones() junto con permuta():

def permutaciones(c, n):
    """Calcula y devuelve una lista con todas las
       permutaciones posibles que se pueden hacer
       con los elementos contenidos en c tomando n
       elementos a la vez.
    """
    return sum([permuta(s)
                for s in combinaciones(c, n)],
               [])

En este código usamos una lista por comprensión para calcular todas las combinaciones de c tomando n elementos a la vez, y luego permutar cada subconjunto obtenido. Finalmente concatenamos con sum() todas las listas con las permutaciones resultantes.

Ejemplo de uso:

>>> permutaciones([1, 2, 3], 0)
[[]]
>>> permutaciones([1, 2, 3], 1)
[[1], [2], [3]]
>>> permutaciones([1, 2, 3], 2)
[[1, 2], [2, 1], [1, 3], [3, 1], [2, 3], [3, 2]]
>>> imprime_ordenado(permutaciones([1, 2, 3, 4], 3))
[1, 2, 3]
[1, 2, 4]
[1, 3, 2]
[1, 3, 4]
[1, 4, 2]
[1, 4, 3]
[2, 1, 3]
[2, 1, 4]
[2, 3, 1]
[2, 3, 4]
[2, 4, 1]
[2, 4, 3]
[3, 1, 2]
[3, 1, 4]
[3, 2, 1]
[3, 2, 4]
[3, 4, 1]
[3, 4, 2]
[4, 1, 2]
[4, 1, 3]
[4, 2, 1]
[4, 2, 3]
[4, 3, 1]
[4, 3, 2]

Si m es la cardinalidad del conjunto original y vamos a formar grupos tomando n elementos a la vez, entonces el número de permutaciones posibles se denota \({}_m P_n\), y su fórmula es: $$ {}_m P_n = \frac{m!}{(m - n)!} $$ La función en Python que realiza este cálculo es:

def numero_permutaciones(m, n):
    """Calcula y devuelve el número de permutaciones
       posibles que se pueden hacer con m elementos
       tomando n elementos a la vez.
    """
    return factorial(m) // factorial(m - n)

Ahora verifiquemos si numero_permutaciones() es consistente con lo que permutaciones() nos devuelve:

>>> len(permutaciones(range(10), 0))
1
>>> len(permutaciones(range(10), 0)) \
    == numero_permutaciones(10, 0)
True
>>> len(permutaciones(range(20), 1))
20
>>> len(permutaciones(range(20), 1)) \
    == numero_permutaciones(20, 1)
True
>>> len(permutaciones(range(20), 4))
116280
>>> len(permutaciones(range(20), 4)) \
    == numero_permutaciones(20, 4)
True

Conclusión

El conjunto potencia y las listas por comprensión son nuestros mejores amigos para divertirnos diseñando algoritmos que generan combinaciones y permutaciones.

☺

Finalmente, vale la pena comentar que el módulo itertools, disponible tanto en Python 2 como en Python 3, proporciona toda la funcionalidad que describimos aquí (y mucha más) en forma de generadores combinatorios. Sugiero utilizar dicho módulo si se tiene interés de obtener combinaciones y/o permutaciones sin tener que conocer los detalles de su implementación:

>>> from itertools import combinations, permutations
>>> list(combinations([1, 2, 3, 4], 3))
[(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)]
>>> list(permutations([1, 2, 3, 4], 3))
[(1, 2, 3), (1, 2, 4), (1, 3, 2), (1, 3, 4), (1, 4, 2), 
(1, 4, 3), (2, 1, 3), (2, 1, 4), (2, 3, 1), (2, 3, 4), 
(2, 4, 1), (2, 4, 3), (3, 1, 2), (3, 1, 4), (3, 2, 1), 
(3, 2, 4), (3, 4, 1), (3, 4, 2), (4, 1, 2), (4, 1, 3), 
(4, 2, 1), (4, 2, 3), (4, 3, 1), (4, 3, 2)]

Potenciando conjuntos

En nuestra entrada anterior conocimos lo que son las listas por comprensión y pudimos ver algunos ejemplos sencillos de cómo utilizarlas en Python. En esta entrada veremos un ejemplo más sofisticado de este tema: el conjunto potencia.

Portada de la sexta edición del libro
“Fundamentos de Matemáticas”
escrito por Juan Silva y Adriana Lazo.
Editorial Limusa, 2003.

Definición de conjunto potencia tal como
aparece en el texto de Silva y Lazo.

El conjunto potencia es un concepto matemático que recuerdo haber aprendido cuando estudiaba el tema de lógica y conjuntos en la preparatoria, por ahí de a mediados de la década de los años ochenta. Si no me falla la memoria, fue en el libro de “Fundamentos de Matemáticas” de Silva y Lazo donde lo vi por primera vez. Sin embargo no fue sino hasta algunos años después que me detuve a pensar cómo es que se podría programar el conjunto potencia en un lenguaje de programación. Esto fue a raíz de un taller sobre el lenguaje de programación Scheme (un dialecto del lenguaje Lisp) al que asistí en la Conferencia Latinoamericana de Informática (CLEI) de 1994 en México, impartido por el Dr. Daniel P. Friedman.

Friedman es profesor e investigador de ciencia de la computación en la Universidad de Indiana. Su principal área de interés es la teoría y aplicación de los lenguajes de programación. Es coautor de varios libros, de los cuales destacan principalmente dos: “The Little Schemer” que escribió junto con Matthias Felleisen, y “Essentials of Programming Languages” que escribió con Mitchell Wand y Christopher Haynes.

El profesor Daniel Paul Friedman.

Portada de la cuarta edición del libro
“The Little Schemer” escrito por
Daniel Friedman y Matthias Felleisen,
Editorial The MIT Press, 1995.

En el taller al que hago alusión, Friedman explicó la función map de Scheme y mostró cómo utilizarla para implementar en seis líneas de código una función recursiva que calculaba el conjunto potencia. La sencillez y elegancia de dicha implementación me deslumbró. En casi una década de programar cotidianamente en lenguajes como Basic, Pascal y C++, nunca había visto algo que se le pareciera.

El código de Scheme que presentó Friendman era algo así:

(define potencia 
  (λ (c)
    (if (null? c)
        '(())
        (let ((r (potencia (cdr c))))
          (append r (map (λ (s) (cons (car c) s)) r))))))

El código anterior puede resultar bastante críptico para quien nunca ha programado en Lisp, sobre todo por el uso aparentemente excesivo de paréntesis (Lisp es el acrónimo de: List Processing [procesamiento de listas], aunque algunas personas dicen que debería significar algo así como: Lost In a Sea of Parenthesis [Perdido en un mar de paréntesis]). Sin embargo esas seis líneas están repletas de conceptos poderosos que son aplicables en otros lenguajes de programación (incluyendo Python).

Lo más importante a recalcar aquí es que la función map de Scheme puede ser emulada directamente en Python a través de listas por comprensión (de hecho, Python también cuenta con la función map(), pero en general son más fáciles se usar las listas por comprensión). Así pues, presento a continuación, usando Python, el razonamiento lógico detrás del código que en su momento me impresionó tan positivamente.

NOTA: Todo el código que se presenta posteriormente fue elaborado para Python 3.5.

Definición matemática

Dado un conjunto C, el conjunto potencia de C es otro conjunto formado exclusivamente por todos los subconjuntos de C.

Por ejemplo, si C = {x, y, z}, el conjunto potencia de C, representado como ℘(C), es: {{}, {x}, {y}, {z}, {x, y}, {x, z}, {y, z}, {x, y, z}}. Vale la pena notar que el conjunto resultante incluye tanto al conjunto vacío {} como al conjunto original C.

Como segundo ejemplo, supongamos que tenemos cuatro monedas con los siguientes valores: 1 peso, 2 pesos, 5 pesos y 10 pesos. ¿De qué maneras distintas podemos combinarlas y qué valor monetario obtenemos en cada caso? Resolviendo con el conjunto potencia:

• {} = 0 pesos
• {1} = 1 peso
• {2} = 2 pesos
• {1, 2} = 3 pesos
• {5} = 5 pesos
• {1, 5} = 6 pesos
• {2, 5} = 7 pesos
• {1, 2, 5} = 8 pesos
• {10} = 10 pesos
• {1, 10} = 11 pesos
• {2, 10} = 12 pesos
• {1, 2, 10} = 13 pesos
• {5, 10} = 15 pesos
• {1, 5, 10} = 16 pesos
• {2, 5, 10} = 17 pesos
• {1, 2, 5, 10} = 18 pesos

Implementación en Python

Python cuenta con un tipo de dato Set, sin embargo para lo que queremos hacer aquí resulta más adecuado usar listas convencionales para representar conjuntos. Así pues, la lista vacía [] representa un conjunto vacío y la lista [1, 2, 3] representa un conjunto de tres elementos.

Definamos un algoritmo recursivo para calcular el conjunto potencia. Como en cualquier procedimiento recursivo, debemos definir dos cosas:

1. El caso base.
2. La llamada recursiva. 

El caso base es bastante sencillo: el conjunto potencia de un conjunto vacío es un conjunto que contiene al conjunto vacío. Es decir, ℘({}) = {{}}. El código de Python correspondiente sería el siguiente (donde c es el parámetro de la función que contendrá el conjunto de entrada):

if len(c) == 0:
    return [[]]

Hay que notar que el valor devuelto no es una lista vacía sino una lista que contiene a la lista vacía.

La llamada recursiva es más elaborada. Para empezar, es importante que dicha llamada corresponda a una versión más simple del problema siendo resuelto. Por ejemplo, si estamos resolviendo ℘({x, y, z}), entonces la llamada recursiva puede involucrar calcular ℘({x, y}), en donde el nuevo problema es más simple dado que consiste en calcular el conjunto potencia para un conjunto con un elemento menos. Al quitar un elemento nos estamos acercando más al caso base, lo cual es indispensable pues debemos tener la certeza de que la función termina en algún momento. El código de Python sería:

r = potencia(c[:-1])

La expresión c[:-1] obtiene una rebanada (slice) de la lista c comenzando implícitamente en el índice 0 (el índice del primer de elemento la lista) y hasta un elemento antes del índice -1 (el índice del último elemento de la lista). Esto efectivamente devuelve una copia de la lista c pero sin su último elemento. El resultado de la llamada recursiva queda en la variable r, la cual usaremos más adelante.

Veamos el siguiente esquema para determinar lo que sigue:

Esta figura muestra que la llamada recursiva ℘({x, y}) devuelve el conjunto {{}, {x}, {y}, {x, y}}. Por definición eso es lo que debe devolver la función conjunto potencia cuando recibe {x, y} como entrada. Aquí debemos suponer (algunos dirían: “como acto de fe”) que la llamada recursiva efectivamente regresa el valor correcto. Lo que resta ahora es responder a la siguiente pregunta: ¿Qué debo hacerle al resultado de la llamada recursiva ℘({x, y}) para convertirlo en el resultado esperado de la llamada original ℘({x, y, z})?

Hay tres cosas que podemos notar en la figura:

1. El resultado de ℘({x, y, z}) tiene el doble de elementos que el resultado de  ℘({x, y}). 
2. La primera mitad del resultado de ℘({x, y, z}) está compuesta exactamente de los mismos elementos del resultado de ℘({x, y}).
3. La segunda mitad del resultado de ℘({x, y, z}) está compuesta por los mismos elementos del resultado de ℘({x, y}) pero añadiendo al final de cada subconjunto el elemento z, que fue el que eliminamos al momento de hacer la llamada recursiva.

Podemos plantear lo anterior en Python así:

return r + [s + [c[-1]] for s in r]

Recordemos que r almacena la lista obtenida como resultado de la llamada recursiva. La expresión [... for s in r] es una lista por comprensión que permite procesar cada subconjunto s contenido en el conjunto r. La expresión s + [c[-1]] dentro de la lista por comprensión añade el último elemento del conjunto c (el elemento del índice -1) a cada subconjunto s. Finalmente, la expresión completa r + [...] concatena las dos mitades del resultado descritas en los puntos 2 y 3 de arriba.

El código completo quedaría así:

def potencia(c):
    """Calcula y devuelve el conjunto potencia del 
       conjunto c.
    """
    if len(c) == 0:
        return [[]]
    r = potencia(c[:-1])
    return r + [s + [c[-1]] for s in r]

Probando el código desde el shell de Python:

>>> potencia([])
[[]]
>>> potencia([1])
[[], [1]]
>>> potencia([1, 2])
[[], [1], [2], [1, 2]]

Obteniendo las soluciones de los problemas que originalmente planteamos en la sección anterior:

>>> potencia(['x', 'y', 'z'])
[[], ['x'], ['y'], ['x', 'y'], ['z'], ['x', 'z'], 
['y', 'z'], ['x', 'y', 'z']]
>>> for t in potencia([1, 2, 5, 10]): 
...     print('{0:13} = ${1:2d}'.format(t, sum(t)))
...
[]            = $ 0
[1]           = $ 1
[2]           = $ 2
[1, 2]        = $ 3
[5]           = $ 5
[1, 5]        = $ 6
[2, 5]        = $ 7
[1, 2, 5]     = $ 8
[10]          = $10
[1, 10]       = $11
[2, 10]       = $12
[1, 2, 10]    = $13
[5, 10]       = $15
[1, 5, 10]    = $16
[2, 5, 10]    = $17
[1, 2, 5, 10] = $18

Cardinalidad del conjunto potencia

La cardinalidad de un conjunto C es el número de elementos que lo conforman, y se denota así: |C|. En Python usamos la función len() para determinar el número de elementos que tiene una lista.

Como vimos en nuestra implementación del conjunto potencia, el caso base devuelve el conjunto {{}}, el cual tiene una cardinalidad igual a 1. En cualquier otro caso se devuelve un conjunto cuya cardinalidad es el doble del resultado de su correspondiente llamada recursiva. En otras palabras, la cardinalidad del conjunto potencia de un conjunto C, denotado como |℘(C)|, es: 2ⁿ, donde n = |C|. Podemos verificar que nuestra función potencia() cumple con esta fórmula:

>>> len(potencia([]))
1
>>> len(potencia([])) == 2 ** 0
True
>>> len(potencia(range(4)))
16
>>> len(potencia(range(4))) == 2 ** 4
True
>>> len(potencia(range(20)))
1048576
>>> len(potencia(range(20))) == 2 ** 20
True

Recordemos que x ** n en Python significa elevar x al exponente n.

Imprimiendo bonito el resultado

El algoritmo que implementamos para calcular el conjunto potencia devuelve una lista de listas que, a primera vista, aparenta no tener orden alguno. Afortunadamente es sencillo definir una función que permita imprimir los elementos de nuestra lista resultante ordenada de manera más conveniente:

def imprime_ordenado(c):
    """Imprime en la salida estándar todos los
       subconjuntos del conjunto c (una lista de
       listas) ordenados primero por tamaño y
       luego lexicográficamente. Cada subconjunto
       se imprime en su propia línea. Los
       elementos de los subconjuntos deben ser
       comparables entre sí, de otra forma puede
       ocurrir un TypeError.
    """
    for e in sorted(c, key=lambda s: (len(s), s)):
        print(e)

La función imprime_ordenado() hace uso a su vez de la función sorted(), la cual está predefinida en Python. A sorted() le mandamos como argumentos la lista a ordenar c (un conjunto de conjuntos) y un objeto lambda (función anónima) que se usa para calcular la llave (key) a utilizar al momento de realizar las comparaciones. En este caso nuestra función anónima crea una tupla con el tamaño de cada sublista s siendo comparada y el valor propio de s. Con esto garantizamos efectivamente que primero se compara la cardinalidad de cada subconjunto y en caso de un empate se compara entonces usando su valor. Por ejemplo, la lista [2, 3] se considera menor que [1, 2, 3] por ser más pequeña. Sin embargo, la lista [2, 3] es menor que la lista [2, 4], ya que son del mismo tamaño pero lexicográficamente (según el orden de diccionario) la primera lista debe ir antes que la segunda.

Ahora ya podemos usar la función imprime_ordenado() de la siguiente manera para producir un resultado más agradable a la vista:

>>> imprime_ordenado(potencia([1, 2, 3, 4]))
[]
[1]
[2]
[3]
[4]
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]
[1, 2, 3]
[1, 2, 4]
[1, 3, 4]
[2, 3, 4]
[1, 2, 3, 4]

Conclusión

En esta entrada diseñamos un algoritmo recursivo para calcular el conjunto potencia. Utilizamos una lista por comprensión para simplificar la manera de procesar listas anidadas dentro una lista más grande. En la próxima entrada del blog de EduPython veremos más aplicaciones de las listas por comprensión y descubriremos que el conjunto potencia se puede utilizar también para generar combinaciones y permutaciones.

Listas por comprensión

Las listas por comprensión (list comprehensions en inglés) son una característica interesante pero sobre todo muy útil que fue incorporada al lenguaje Python en su versión 2.0 en el año 2000 (ver: PEP 202). Las listas por comprensión debutaron por vez primera en un lenguaje de programación en 1977, cuando Rod Burstall y John Darlington diseñaron el lenguaje funcional NPL. Hoy en día un número importante de lenguajes de programación soportan listas por comprensión, incluyendo: Haskell, JavaScript 1.7, CoffeeScript, Erlang, F#, Scala, Clojure y Racket.

Conjuntos

Para entender qué son las listas por comprensión tomémonos primero un momento para responder a las siguientes dos preguntas:

¿Qué es un conjunto?

¿Cómo se define un conjunto?

En matemáticas, un conjunto es simplemente una colección de elementos bien definidos. Los elementos pueden ser cualquier cosa: números, letras, nombres, figuras, etc.

Un conjunto de figuras geométricas.

Hay dos manera de definir los elementos que pertenecen a un conjunto: por extensión o por comprensión. Cuando se define un conjunto por extensión cada elemento se enumera de manera explícita. Tomando un ejemplo inspirado en el Señor de los Anillos de J. R. R. Tolkien:

A = { Frodo, Sam, Pippin, Merry, Legolas,
          Gimli, Aragorn, Boromir, Gandalf }

La expresión anterior se lee así: A es el conjunto formado por los elementos: Frodo, Sam, Pippin, Merry, Legolas, Gimli, Aragorn, Boromir y Gandalf.

Por otro lado, cuando un conjunto se define por comprensión no se mencionan los elementos uno por uno sino que se indica una propiedad que todos éstos cumplen, por ejemplo:

B = { x | x ∈ Comunidad del Anillo }

La expresión de arriba se lee así: B es el conjunto de elementos x tales que x  pertenece a la Comunidad del Anillo.

Los miembros de la Comunidad del Anillo.
Imagen de haleyhss.

En los ejemplos anteriores podemos decir que el conjunto A es igual al conjunto B debido a que ambos tienen exactamente los mismos elementos. Esto es cierto aún a pesar de que el conjunto A se definió por extensión y B se definió por comprensión.

Listas en Python

En Python se utilizan listas para representar colecciones de elementos (a decir verdad, Python cuenta también con conjuntos [sets] y otros tipos de datos, pero con el fin de simplificar la discusión usaremos aquí solo listas). Ahora bien, así como matemáticamente se pueden definir los conjuntos por extensión y por comprensión, en Python se pueden definir las listas también por extensión y por comprensión.

Supongamos, como ejemplo, que deseamos tener una lista con las primeras diez potencias de dos. Una potencia de dos es cualquiera de los números obtenidos al elevar el número dos a una potencia entera no negativa. El siguiente código de Python cumple con este cometido usando una lista por extensión (enumerando todos los elementos de manera explícita):

c = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512]

Matemáticamente, un conjunto equivalente a la lista anterior se puede definir por comprensión así:

C = { 2^x | x ∈ ℤ ∧ 0 ≤ x < 10 }

Esta expresión se lee así: C es el conjunto formado por los elementos “dos elevado a la x”, tales que x pertenece al conjunto de los números enteros y además x es mayor o igual a 0 pero menor que 10. Dicha expresión se puede traducir a Python directamente usando la notación de listas por comprensión:

c = [2 ** x for x in range(10)]

La sintaxis de listas por comprensión consiste en colocar entre corchetes una expresión (2 ** x) seguida de una cláusula for. Dicha cláusula es muy similar en intención a un ciclo for convencional. En este caso estamos indicando que la variable x tomará los valores devueltos por la función range(10) (los enteros del 0 al 9). A partir de cada valor que toma x se calcula el resultado de la expresión 2 ** x y con eso se determinan los valores finales de la lista resultante. Esencialmente es como si ejecutáramos el siguiente código:

c = []
for x in range(10):
    c.append(2 ** x)

El método append() se encarga de ir añadiendo un nuevo elemento (2 ** x) al final de la lista c (inicialmente vacía) en cada iteración del ciclo for.

Las listas por comprensión también pueden incluir una expresión condicional que nos permite quedarnos con ciertos elementos y eliminar los restantes. Para ello se debe utilizar una cláusula if después de la cláusula for. Por ejemplo, si queremos una lista con todos los números entre 1 y 100 que sean múltiplos de 7 o que terminen con el dígito 7, podemos escribir la siguiente lista por comprensión:

[n for n in range(1, 101) if n % 7 == 0 or n % 10 == 7]

El resultado es exactamente lo esperado:

[ 7, 14, 17, 21, 27, 28, 35, 37, 42, 47, 
 49, 56, 57, 63, 67, 70, 77, 84, 87, 91,
 97, 98]

Puede haber más de una cláusula for en una lista por comprensión (también puede haber cero o más cláusulas if). Como ejemplo, supongamos que tengo cuatro camisas (de color rojo, amarillo, azul y verde) y dos pantalones (de color negro y blanco). Quiero saber de qué manera puedo combinar mi ropa. Hay ocho formas distintas (4 camisas × 2 pantalones) de hacerlo:

• Camisa roja y pantalón negro.
• Camisa roja y pantalón blanco.
• Camisa amarilla y pantalón negro.
• Camisa amarilla y pantalón blanco.
• Camisa azul y pantalón negro.
• Camisa azul y pantalón blanco.
• Camisa verde y pantalón negro.
• Camisa verde y pantalón blanco.

Podemos dejar que Python calcule lo anterior usando una lista por comprensión con dos cláusulas for:

[(camisa, pantalon) 
    for camisa in ['rojo', 'amarillo', 'azul', 'verde']
    for pantalon in ['negro', 'blanco']]

El resultado es:

[('rojo', 'negro'), ('rojo', 'blanco'),
 ('amarillo', 'negro'), ('amarillo', 'blanco'),
 ('azul', 'negro'), ('azul', 'blanco'),
 ('verde', 'negro'), ('verde', 'blanco')]

Como se puede observar, el ejemplo anterior calcula efectivamente el producto cartesiano de dos conjuntos (los colores de las camisas por los colores de los pantalones).

A continuación presento un par de ejemplos que demuestran el uso de listas por comprensión en contextos más complejos.

Ternas pitagóricas

Una terna pitagórica es un conjunto de tres números naturales (a, b, c) que cumplen con la siguiente ecuación:

 a² + b² = c²

El nombre deriva del teorema de Pitágoras, el cual establece que el cuadrado de la hipotenusa de un triángulo rectángulo es igual a la suma de los cuadrados de sus catetos.

Por ejemplo, si a = 3, b = 4 y c = 5, tenemos:

3² + 4² = 5²
9 + 16 = 25

Queremos calcular todos los valores posibles de a, b y c que sean menores a 20. Traduciendo directamente los requisitos planteados a una lista por comprensión, tenemos:

[(a, b, c)
    for a in range(1, 20)
    for b in range(1, 20)
    for c in range(1, 20)
    if a ** 2 + b ** 2 == c ** 2]

El resultado correspondiente es:

[(3, 4, 5), (4, 3, 5), (5, 12, 13), (6, 8, 10),
 (8, 6, 10), (8, 15, 17), (9, 12, 15), (12, 5, 13),
 (12, 9, 15), (15, 8, 17)]

Si observamos con atención, hay el doble de ternas deseables en esta solución, ya que se incluyen todas las permutaciones de los valores de a y b. Para corregir esta situación podemos cambiar el valor de inicio de los rangos de las variables b y c para que inicien en un valor posterior al contenido en la variable que está a su izquierda inmediata y con ello garantizamos que: a < b < c. El código quedaría así:

[(a, b, c)
    for a in range(1, 20)
    for b in range(a + 1, 20)
    for c in range(b + 1, 20)
    if a ** 2 + b ** 2 == c ** 2]

El resultado final es el buscado:

[(3, 4, 5), (5, 12, 13), (6, 8, 10), (8, 15, 17),
 (9, 12, 15)]

La lista por comprensión anterior es equivalente al resultado que queda en la variable r del siguiente código:

r = []
for a in range(1, 20):
    for b in range(a + 1, 20):
        for c in range(b + 1, 20):
            if a ** 2 + b ** 2 == c ** 2:
                r.append((a, b, c))

Vale la pena notar que se preserva el orden de los fors y el if en ambos códigos.

Creando una matriz identidad

Una matriz identidad de tamaño n es una matriz cuadrada de n renglones por n columnas en donde cada elemento de la diagonal principal es 1 y los elementos restantes son 0. Por ejemplo, la siguiente es una matriz identidad de tamaño n = 4:

En álgebra lineal una matriz identidad sirve como elemento neutro en la multiplicación de matrices. Esto quiere decir que si I es la matriz identidad y A es otra matriz de dimensiones compatibles, entonces A I = A.

La siguiente función de Python crea un matriz identidad de tamaño n usando listas por comprensión:

def matriz_identidad(n):
    """Devuelve una lista de listas que representa una
       matriz identidad de tamaño n."""
    return [[(1 if ren == col else 0) for col in range(n)]
            for ren in range(n)]

En el código anterior hay una lista por comprensión anidada dentro de otra lista por comprensión. La comprensión externa controla que los n renglones sean generados. La comprensión interna crea un renglón particular conformado por n columnas. Cada elemento individual de la matriz es 1 si su número de renglón es igual a su número de columna, de otra forma es 0.

Probando el código:

>>> matriz_identidad(5)
[[1, 0, 0, 0, 0], 
 [0, 1, 0, 0, 0], 
 [0, 0, 1, 0, 0], 
 [0, 0, 0, 1, 0], 
 [0, 0, 0, 0, 1]]
>>> matriz_identidad(1)
[[1]]

Conclusión

Las listas por comprensión en Python se basan en la notación matemática de conjuntos por comprensión. La notación no es muy compleja, y una vez que la entendemos podemos escribir programas más compactos y expresivos.

Pensando en Pi

Esta entrada del blog de EduPython es para conmemorar el Día de Pi del 2015. Como sabemos, π (pi) es una constante matemática que representa la relación que tiene la circunferencia de un círculo con respecto a su diámetro. Esta constante es un número irracional, lo cual significa que posee infinitas cifras decimales sin patrón de repetición alguno. Sus primeras 20 cifras son: 3.1415926535897932384

A partir del formato de fechas que se utiliza en los Estados Unidos (mes/día), el Día de Pi se celebra el 3/14 (14 de marzo) de cada año, dado que 3, 1 y 4 son los tres primeros dígitos de π. Como feliz coincidencia, el 14 de marzo también fue el día que nació Albert Einstein.

Albert Einstein nació el 14 de marzo
de 1879 en Ulm, Alemania.

Cada año el Día de Pi^§ es celebrado por estudiantes, profesores, matemáticos y computólogos alrededor de todo el mundo. Las celebraciones incluyen actividades relacionadas con π, por ejemplo: competencias que consisten en recitar de memoria la mayor cantidad de dígitos de π, hornear y comer pays^¶, o incluso hasta escribir blogs acerca de π.

Un rico pay de moras para celebrar el Día de Pi.

El año 2015 es especial, ya que una sola vez en cada siglo ocurre una combinación de fecha y hora que incluye los primeros 10 dígitos de π: 3/14/15 a las 9:26:53 hrs.

Poniéndonos más técnicos, veremos ahora algunas formas de obtener valores aproximados de π usando Python. Todos los ejemplos de código que se muestran a continuación fueron probados en Python versión 3.4.

π en Python

Cuando un programa en Python requiere del valor de π, lo más recomendable es importarlo directamente del módulo math. Podemos demostrar esto usando el shell de Python:

>>> from math import pi
>>> pi
3.141592653589793

Podemos observar que el valor de pi tiene 16 dígitos de precisión debido a que ésa es la cantidad máxima de dígitos decimales que puede tener un valor de tipo float en una implementación típica de Python. Dicha aproximación es adecuada para la mayoría de los cálculos que comúnmente se requieren en las diversas disciplinas científicas e ingenieriles.

A continuación exploraremos diversas formas programáticas de calcular el valor de π.

π como un número racional

En tiempos pasados era común usar fracciones como 22/7 y 355/113 para aproximar π.  Evaluando dichas fracciones:

>>> 22/7
3.142857142857143
>>> 355/113
3.1415929203539825

Podemos notar que el resultado de la primera división tiene tres dígitos correctos, mientras que la segunda división acierta siete.

Serie de Gregory–Leibniz

En los siglos XVII y XVIII, James Gregory y Gottfried Leibniz descubrieron^† una serie infinita que sirve para calcular π: $$ \begin{align*} \pi &= 4 \left ( \sum_{k=1}^{\infty} \frac{(-1)^{(k + 1)}}{2k-1} \right ) \\ &= 4 \left ( 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \frac{1}{11} \cdots \right ) \end{align*} $$ Resulta relativamente fácil traducir esta fórmula a una función en Python que nos permita aproximar el valor de π tomando los n primeros términos de la serie:

def gregory_leibniz(n):
    """Calcula y devuelve el valor de pi usando
    los primeros n términos de la serie de
    Gregory–Leibniz.
    
    π = 4(1 - 1/3 + 1/5 - 1/7 + ...)
    """
    s = 0
    for k in range(1, n + 1):
        s += (-1)**(k + 1) / (2 * k - 1)
    return 4 * s

Podemos probar desde el shell nuestra función con diferentes valores de n:

>>> gregory_leibniz(1)
4.0
>>> gregory_leibniz(10)
3.0418396189294032
>>> gregory_leibniz(100)
3.1315929035585537
>>> gregory_leibniz(1000)
3.140592653839794
>>> gregory_leibniz(10000)
3.1414926535900345
>>> gregory_leibniz(100000)
3.1415826535897198
>>> gregory_leibniz(1000000)
3.1415916535897743
>>> gregory_leibniz(10000000)
3.1415925535897915

Se puede observar que se obtiene un dígito más de precisión cada vez que multiplicamos por 10 el argumento de la función gregory_leibniz(). Sin embargo, tal como es de esperarse, entre más grande sea el valor de n más tiempo tarda la función en completar.

Método de Montecarlo

Supongamos que tenemos un tablero para jugar a los dardos. Dicho tablero es un cuadrado de un metro de cada lado. El tablero contiene un cuarto de círculo tal como se muestra en la siguiente imagen:

Ahora lanzamos t dardos al tablero. Si todos los dardos caen de manera aleatoria dentro del tablero con una distribución uniforme, algunos dardos caerán en la región gris y otros en la región blanca. Los puntos en la siguiente imagen representan los lugares donde pudieron haber caído los t dardos:

Sabemos que el área del tablero es: 1 metro × 1 metro = 1 metro². El área de la región gris es un cuarto del área de un círculo, es decir: \((\pi \times r^2) \div 4\). Dado que el radio r mide 1 metro (lo que mide un lado del tablero), entonces el área de la región gris es: $$ \frac{\pi \times (1\;\textrm{m}^2)}{4} = \frac{\pi}{4} \textrm{m}^2 $$ Si t es el total de dardos lanzados y g es la cantidad de esos dardos que cayeron dentro de la región gris, entonces podemos considerar que \(g/t\) debe ser aproximadamente igual a la división del área del cuarto de círculo (π/4 metros²) entre el área de todo el tablero (1 metro²): $$ \frac{\frac{\pi}{4} \textrm{m}^2}{1 \; \textrm{m}^2} \approx \frac{g}{t} $$ $$ \frac{\pi}{4} \approx \frac{g}{t} $$ $$ \pi \approx \frac{4 g}{t} $$ Haciendo los despejos necesarios, obtenemos otra forma de aproximar el valor de π. El método de Montecarlo para calcular π quedaría así:

• Inicializar g en 0.
• Repetir t veces lo siguiente:
• Generar de manera aleatoria un punto con coordenadas (x, y) dentro del área del tablero. Dado que que cada lado del tablero mide un metro, los valores de x y y deben ser números reales entre 0 y 1.
• Calcular la distancia d entre el centro del círculo (la esquina inferior izquierda del tablero) y el punto (x, y). Para ello utilizamos el teorema de Pitágoras: \(d = \sqrt{x^2 + y^2}\)
• Si d es menor a un metro (el radio del círculo) entonces el punto (x, y) está dentro del área del círculo. En ese caso incrementamos en uno el valor de g.
• El valor aproximado de π es: \( \frac{4 g}{t} \).

Se le llama Montecarlo a este método en referencia al casino que se encuentra en Mónaco, el cual es considerado por muchos como la capital mundial de los juegos de azar.

Casino de Montecarlo en el Principado
de Mónaco.

La implementación del método de Montecarlo en Python es bastante directa. Para generar los números aleatorios usamos la función random() del módulo random. Dicha función devuelve un número de punto flotante al azar en el intervalo [0, 1), es decir, dicho número es mayor o igual a cero pero menor a uno. Y esto es justo lo que nuestro algoritmo necesita. El código quedaría así:

from random import random
from math import sqrt

def montecarlo(t):
    """Calcula y devuelve el valor aproximado
    de pi usando el método de Montecarlo a
    partir de t puntos.
    """
    g = 0
    for i in range(t):
        x = random()
        y = random()
        d = sqrt(x ** 2 + y ** 2)
        if d < 1:
            g += 1
    return 4 * g / t

Llamando la función con t = 100,000,000 (cien millones) podemos obtener casi cinco dígitos de precisión:

>>> montecarlo(100000000)
3.14166808
>>> montecarlo(100000000)
3.14167544
>>> montecarlo(100000000)
3.1415684

Debido al uso de números aleatorios, cada invocación a la función montecarlo() produce resultados (ligeramente) diferentes.

Calculando muchos dígitos de π

Un problema que tienen todos los esquemas discutidos anteriormente para el cálculo de π es que están limitados a la representación interna del tipo float de Python. La mayoría de las implementaciones de Python usan números de punto flotante de precisión doble de 64 bits tal como lo describe el estándar IEEE 754 (también conocido como IEC 60559). Para fines prácticos, este tipo de dato permite representar valores numéricos de 15 a 17 dígitos decimales significativos. Pero, ¿qué pasa si queremos calcular más dígitos de π? Para eso tenemos los algoritmos de espita (spigot en inglés) que únicamente requieren números enteros. Así como gotea el agua de una espita (grifo, válvula o llave), a los algoritmos de espita se les llama así debido a que producen dígitos individuales de π que no se reusan después de que son calculados. Esto contrasta con las series infinitas o algoritmos iterativos, en los que se retienen y utilizan todos los dígitos intermedios hasta que se produce el resultado final.

Una espita
(grifo, válvula o llave)

El primer algoritmo de espita se le atribuye a A.H.J. Sale, que en 1968 presentó una manera de calcular muchos dígitos de e. En 1995 Stan Wagon y Stanley Rabinowitz publicaron un algoritmo de espita que permite calcular una cantidad arbitraria de dígitos de π.

Boris Gourévitch hace un buen trabajo explicando el algoritmo de espita para calcular π. Dado que el tema es algo complicado, no entraré en más detalles aquí. Solo me limitaré a mostrar un programa escrito por John Zelle en el 2006 que implementa el algoritmo de espita publicado un año antes por Jeremy Gibbons que a su vez mejoró el trabajo original de Wagon y Rabinowitz.

def espita(d):
    """Regresa una lista con los primeros d dígitos
       de pi utilizando el algoritmo de espita
       diseñado por Jeremy Gibbons. Implementación
       de John Zelle con ligeras alteraciones por
       Ariel Ortiz.
    """
    x = []
    q,r,t,k,n,l = 1,0,1,1,3,3
    while len(x) < d:
        if 4*q+r-t < n*t:
            x.append(n)
            q,r,t,k,n,l = (
                10*q,10*(r-n*t),t,k,
                (10*(3*q+r))//t-10*n,l)
        else:
            q,r,t,k,n,l = (
                q*k,(2*q+r)*l,t*l,k+1,
                (q*(7*k+2)+r*l)//(t*l),l+2)
    return x

Probando la función:

>>> espita(1)
[3]
>>> espita(3)
[3, 1, 4]
>>> espita(5)
[3, 1, 4, 1, 5]
>>> espita(10)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3]
>>> espita(20)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 
 3, 2, 3, 8, 4]
>>> espita(100)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 
 3, 2, 3, 8, 4, 6, 2, 6, 4, 3, 3, 8, 3, 2, 7, 
 9, 5, 0, 2, 8, 8, 4, 1, 9, 7, 1, 6, 9, 3, 9, 
 9, 3, 7, 5, 1, 0, 5, 8, 2, 0, 9, 7, 4, 9, 4, 
 4, 5, 9, 2, 3, 0, 7, 8, 1, 6, 4, 0, 6, 2, 8, 
 6, 2, 0, 8, 9, 9, 8, 6, 2, 8, 0, 3, 4, 8, 2, 
 5, 3, 4, 2, 1, 1, 7, 0, 6, 7]
>>> espita(1000)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 
 3, 2, 3, 8, 4, 6, 2, 6, 4, 3, 3, 8, 3, 2, 7, 
 9, 5, 0, 2, 8, 8, 4, 1, 9, 7, 1, 6, 9, 3, 9, 
 9, 3, 7, 5, 1, 0, 5, 8, 2, 0, 9, 7, 4, 9, 4, 
 4, 5, 9, 2, 3, 0, 7, 8, 1, 6, 4, 0, 6, 2, 8, 
 6, 2, 0, 8, 9, 9, 8, 6, 2, 8, 0, 3, 4, 8, 2, 
 5, 3, 4, 2, 1, 1, 7, 0, 6, 7, 9, 8, 2, 1, 4, 
 8, 0, 8, 6, 5, 1, 3, 2, 8, 2, 3, 0, 6, 6, 4, 
 7, 0, 9, 3, 8, 4, 4, 6, 0, 9, 5, 5, 0, 5, 8, 
 2, 2, 3, 1, 7, 2, 5, 3, 5, 9, 4, 0, 8, 1, 2, 
 8, 4, 8, 1, 1, 1, 7, 4, 5, 0, 2, 8, 4, 1, 0, 
 2, 7, 0, 1, 9, 3, 8, 5, 2, 1, 1, 0, 5, 5, 5, 
 9, 6, 4, 4, 6, 2, 2, 9, 4, 8, 9, 5, 4, 9, 3, 
 0, 3, 8, 1, 9, 6, 4, 4, 2, 8, 8, 1, 0, 9, 7, 
 5, 6, 6, 5, 9, 3, 3, 4, 4, 6, 1, 2, 8, 4, 7, 
 5, 6, 4, 8, 2, 3, 3, 7, 8, 6, 7, 8, 3, 1, 6, 
 5, 2, 7, 1, 2, 0, 1, 9, 0, 9, 1, 4, 5, 6, 4, 
 8, 5, 6, 6, 9, 2, 3, 4, 6, 0, 3, 4, 8, 6, 1, 
 0, 4, 5, 4, 3, 2, 6, 6, 4, 8, 2, 1, 3, 3, 9,
 3, 6, 0, 7, 2, 6, 0, 2, 4, 9, 1, 4, 1, 2, 7,
 3, 7, 2, 4, 5, 8, 7, 0, 0, 6, 6, 0, 6, 3, 1,
 5, 5, 8, 8, 1, 7, 4, 8, 8, 1, 5, 2, 0, 9, 2,
 0, 9, 6, 2, 8, 2, 9, 2, 5, 4, 0, 9, 1, 7, 1,
 5, 3, 6, 4, 3, 6, 7, 8, 9, 2, 5, 9, 0, 3, 6,
 0, 0, 1, 1, 3, 3, 0, 5, 3, 0, 5, 4, 8, 8, 2,
 0, 4, 6, 6, 5, 2, 1, 3, 8, 4, 1, 4, 6, 9, 5,
 1, 9, 4, 1, 5, 1, 1, 6, 0, 9, 4, 3, 3, 0, 5,
 7, 2, 7, 0, 3, 6, 5, 7, 5, 9, 5, 9, 1, 9, 5,
 3, 0, 9, 2, 1, 8, 6, 1, 1, 7, 3, 8, 1, 9, 3,
 2, 6, 1, 1, 7, 9, 3, 1, 0, 5, 1, 1, 8, 5, 4,
 8, 0, 7, 4, 4, 6, 2, 3, 7, 9, 9, 6, 2, 7, 4,
 9, 5, 6, 7, 3, 5, 1, 8, 8, 5, 7, 5, 2, 7, 2,
 4, 8, 9, 1, 2, 2, 7, 9, 3, 8, 1, 8, 3, 0, 1,
 1, 9, 4, 9, 1, 2, 9, 8, 3, 3, 6, 7, 3, 3, 6,
 2, 4, 4, 0, 6, 5, 6, 6, 4, 3, 0, 8, 6, 0, 2,
 1, 3, 9, 4, 9, 4, 6, 3, 9, 5, 2, 2, 4, 7, 3,
 7, 1, 9, 0, 7, 0, 2, 1, 7, 9, 8, 6, 0, 9, 4,
 3, 7, 0, 2, 7, 7, 0, 5, 3, 9, 2, 1, 7, 1, 7,
 6, 2, 9, 3, 1, 7, 6, 7, 5, 2, 3, 8, 4, 6, 7,
 4, 8, 1, 8, 4, 6, 7, 6, 6, 9, 4, 0, 5, 1, 3,
 2, 0, 0, 0, 5, 6, 8, 1, 2, 7, 1, 4, 5, 2, 6,
 3, 5, 6, 0, 8, 2, 7, 7, 8, 5, 7, 7, 1, 3, 4,
 2, 7, 5, 7, 7, 8, 9, 6, 0, 9, 1, 7, 3, 6, 3,
 7, 1, 7, 8, 7, 2, 1, 4, 6, 8, 4, 4, 0, 9, 0,
 1, 2, 2, 4, 9, 5, 3, 4, 3, 0, 1, 4, 6, 5, 4,
 9, 5, 8, 5, 3, 7, 1, 0, 5, 0, 7, 9, 2, 2, 7,
 9, 6, 8, 9, 2, 5, 8, 9, 2, 3, 5, 4, 2, 0, 1,
 9, 9, 5, 6, 1, 1, 2, 1, 2, 9, 0, 2, 1, 9, 6,
 0, 8, 6, 4, 0, 3, 4, 4, 1, 8, 1, 5, 9, 8, 1,
 3, 6, 2, 9, 7, 7, 4, 7, 7, 1, 3, 0, 9, 9, 6,
 0, 5, 1, 8, 7, 0, 7, 2, 1, 1, 3, 4, 9, 9, 9,
 9, 9, 9, 8, 3, 7, 2, 9, 7, 8, 0, 4, 9, 9, 5,
 1, 0, 5, 9, 7, 3, 1, 7, 3, 2, 8, 1, 6, 0, 9,
 6, 3, 1, 8, 5, 9, 5, 0, 2, 4, 4, 5, 9, 4, 5,
 5, 3, 4, 6, 9, 0, 8, 3, 0, 2, 6, 4, 2, 5, 2,
 2, 3, 0, 8, 2, 5, 3, 3, 4, 4, 6, 8, 5, 0, 3,
 5, 2, 6, 1, 9, 3, 1, 1, 8, 8, 1, 7, 1, 0, 1,
 0, 0, 0, 3, 1, 3, 7, 8, 3, 8, 7, 5, 2, 8, 8,
 6, 5, 8, 7, 5, 3, 3, 2, 0, 8, 3, 8, 1, 4, 2,
 0, 6, 1, 7, 1, 7, 7, 6, 6, 9, 1, 4, 7, 3, 0,
 3, 5, 9, 8, 2, 5, 3, 4, 9, 0, 4, 2, 8, 7, 5,
 5, 4, 6, 8, 7, 3, 1, 1, 5, 9, 5, 6, 2, 8, 6,
 3, 8, 8, 2, 3, 5, 3, 7, 8, 7, 5, 9, 3, 7, 5,
 1, 9, 5, 7, 7, 8, 1, 8, 5, 7, 7, 8, 0, 5, 3,
 2, 1, 7, 1, 2, 2, 6, 8, 0, 6, 6, 1, 3, 0, 0,
 1, 9, 2, 7, 8, 7, 6, 6, 1, 1, 1, 9, 5, 9, 0,
 9, 2, 1, 6, 4, 2, 0, 1, 9, 8]

Al momento de estar escribiendo esta entrada, el mayor número de dígitos que se han calculado de π son 12.1 billones (1.21 × 10¹³) por Alexander Yee y Shigeru Kondo el 28 de diciembre del 2013.  Aplicando el algoritmo de Chudnovsky tardaron 94 días en total usando un equipo de cómputo con dos procesadores Intel Xeon E5-2690 a 2.9 GHz con 16 núcleos en total, 128 GB de memoria RAM y más de 60 TB de espacio en varios discos duros. Eso sí es como para pensar en π.

¡Feliz día de Pi!

---

Notas

^§ El Día de Pi fue fundado por Larry Shaw y se celebró por primera vez en 1988 en el Exploratorium de San Francisco, California.
Referencia: http://en.wikipedia.org/wiki/Pi_Day.

^¶ En México pay es la castellanización de la palabra inglesa pie. En España se le llama tarta, mientras que en otros países de América Latina se le conoce como pastel. A fin de cuentas la broma consiste en que en inglés pie y π se pronuncian igual.

^† Realmente la serie Gregory–Leibniz fue descubierta tres siglos antes por el matemático indio Madhava de Sangamagrama. Por eso también se le conoce a esta fórmula como la serie de Madhava-Leibniz.