Ejemplo práctico. Potencias de Fermi-Dirac - Evaluación perezosa en python - Parte 6

Se llaman potencias de Fermi-Dirac a los números de la forma \(p^{2^k}\), ordenados de menor a mayor, donde p es un número primo y k es un número natural.

Vamos a ver cómo crear la sucesión de potencias Fermi-Dirac. Realizaremos las siguientes comprobaciones:

potencias: list[int]

potencias[:14]    ==  [2,3,4,5,7,9,11,13,16,17,19,23,25,29]
potencias[60]     ==  241
potencias[10**6]  ==  15476303

Estudio previo

Si sacamos la lista de potencias en función del exponente k tendríamos las siguientes sucesiones:

Necesitamos combinar estas sucesiones en una sola. A priori, no sabemos cuántos elementos vamos a necesitar de cada sucesión. Como máximo, para sacar las primeras 14 potencias nos basta con los primeros 14 números primos y crear 14 secuencias, de \(P_0\) a \(P_{13}\), ordenarlos sus elementos en una única lista y escoger los primeros 14 elementos. Con este proceso habremos calculado 196 potencias para sólo 14 elementos que necesitamos al final.

from primes import primes

potencias = sorted(p**2**k for p in primes[:14] for k in range(0, 14))
print(potencias[:14])

[2, 3, 4, 5, 7, 9, 11, 13, 16, 17, 19, 23, 25, 29]

Aún en el caso de que tuviéramos algún medio de reducir el número de elementos a usar de cada secuencia, seguimos sin saber cuántos números primos serán necesarios. Para sacar los 14 primeros elementos de las potencias de Fermi-Dirac sólo se necesitaban los 10 primeros números primos.

Es evidente que una estrategia por fuerza bruta es complicada y termina por hacer muchos cálculos innecesarios, una complejidad del \(O({n^2})\) no resoluble con un ordenador normal. Veamos cómo nos puede ayudar la evaluación perezosa.

Modelos

Por intentar crear un modelo, intentemos ver las sucesiones como un iterador de iteradores:

from itertools import count

from primes import primes

potencias = ((p**2**k for p in primes) for k in count())

Pero el problema con las expresiones generadora es similar al que tienen las expresiones lambda: carecen de su propia clausura y cualquier variable libre queda alterada por el entorno donde se evalúan.

Se puede comprobar el fallo si intentamos extraer dos iteradores:

p0 = next(potencias)
p1 = next(potencias)
next(p1)  # --> 4
next(p0)  # --> 4
next(p0)  # --> 9

El exponente k ha cambiado de valor con el segundo iterador, lo que afecta a las potencias del primero. Tenemos que dotar a los iteradores de su propia clausura:

from collections.abc import Iterator
from itertools import count

from primes import primes

def potencias_gen(k: int) -> Iterator[int]:
    yield from (p**2**k for p in primes)

potencias = (potencias_gen(k) for k in count())

Para obtener una única secuencia a partir de este iterador de iteradores en un único iterador, operación que se conoce como “aplanar la secuencia”.

Definimos la siguiente función para mezclar dos listas ordenadas:

# tipo para secuencias ordenadas
SortedIterator = Iterator[int]

def zipsort(s1: SortedIterator, s2: SortedIterator) -> SortedIterator:
    x = next(s1)
    y = next(s2)
    while True:
        if x <= y:
            yield x
            x = next(s1)
        else:
            yield y
            y = next(s2)

La función zipsort combina dos listas ordenadas SortedIterator para devolver otra lista ordenada SortedIterator. Si quisiéramos combinar tres listas, bastaría con volver repetir con zipsort:

zipsort(zipsort(s1, s2), s3)

En general, podríamos combinar todas las listas de esta manera:

def flat(iterators: Iterator[SortedIterator]) -> SortedIterator:
    it1 = next(iterators)
    it2 = flat(iterators)
    yield from zipsort(it1, it2)

potencias = flat(potencias_gen(k) for k in count())

El problema es que se entra en un bucle infinito de llamadas recursivas a flat que habrá que evitar.

Si observamos las sucesiones \(P_0\), \(P_1\), \(P_2\),…, el primer elemento de una sucesión es siempre inferior a cualquier elemento de sucesiones posteriores. Usando esta propiedad, podemos redefinir nuestra función aplanadora:

def flat(iterators: Iterator[SortedIterator]) -> SortedIterator:
    it1 = next(iterators)
    yield next(it1)
    yield from zipsort(it1, flat(iterators))

potencias = flat(potencias_gen(k) for k in count())

La función flat devuelve siempre un elemento antes de invocarse por recursividad, suficiente para frenar la cadena de llamadas recursivas. Visto de otro modo, se ha convertido la función en perezosa, devolviendo elementos a medida que sean necesarios. De todos modos, seguimos limitados por el nivel de recursividad en python (~3000 niveles en CPython), aunque no vamos a superar este límite en las pruebas¹.

Código final

Descarga: potencias.py

from collections.abc import Iterator
from itertools import count
from typing import TypeVar

from lazyseq import LazySortedSequence
from primes import primes


SortedIterator = Iterator[int]


def join(s1: SortedIterator, s2: SortedIterator) -> SortedIterator:
    x = next(s1)
    y = next(s2)
    while True:
        if x <= y:
            yield x
            x = next(s1)
        else:
            yield y
            y = next(s2)


def flat(it: Iterator[SortedIterator]) -> SortedIterator:
    s1 = next(it)
    yield next(s1)
    yield from join(s1, flat(it))


def mkiter(k):
    yield from (p ** 2 ** k for p in primes)

potencias = LazySortedSequence(flat(mkiter(k) for k in count()))

Para las comprobaciones:

>>> potencias[:14]
[2, 3, 4, 5, 7, 9, 11, 13, 16, 17, 19, 23, 25, 29]
>>> potencias[60]
241
>>> potencias[10 ** 6]
15476303
>>> primes.size
999432

Para obtener el elemento \(10^6\) tarda bastante al necesitar obtener casi un millón de números primos. Una vez obtenidos, el cálculo es bastante rápido.

Serie Evaluación Perezosa en Python

• Parte 1 - Introducción a la evaluación perezosa
• Parte 2 - Secuencias infinitas
• Parte 3 - Memoización
• Parte 4 - Evaluación perezosa avanzada
• Parte 5 - Formalización de la Secuencia Perezosa
• Parte 6 - Ejemplo práctico: Potencias de Fermi-Dirac
• Apéndice: sobre el tipado de datos utilizado

ANOTACIONES: