Pilas de peniques

import math
import numpy as np
from prettytable import PrettyTable

# Función para el cálculo de números pentagonales.
def omega(n):
    return (3*n**2 - n) // 2

def p(n):
    #
    # Designaremos por 'P' a la lista cuyo k-ésimo término será 
    # el valor de p(k).
    P = [1]
    # Su inicialización obedece al hecho de que se ha definido 
    # P[0] = 1.
    # -------------------------------------------------------------  
    # Ahora bien, como pretendemos determinar el resto de sus 
    # elementos en virtud de la relación
    #
    # P[n] = P[n-1] + P[n-2] - P[n-5] - P[n-7] + P[n-12] + ...
    #        + P[n-15] - P[n-22] - P[n-26] + ...
    #
    # debemos conocer primero los valores de las diferencias
    # n-1, n-2, n-5, n-7, n-12, n-15, n-22, n-26, ...
    #
    # Así pues, llamaremos 'argumentos_no_negativos' a una lista
    # a su vez conformada por 'n' sublistas, cuyos elementos serán 
    # las sucesiones de las diferencias entre '1, 2, ..., n' y el  
    # k-ésimo número pentagonal generalizado, siempre que estas  
    # sean positivas o, a lo sumo, nulas.
    # -------------------------------------------------------------   
    # Por otro lado, precisamos de un adecuado manejo en la 
    # alternancia de signos. 
    # Como los sucesivos valores de la función de partición yacen 
    # en una lista, bastará multiplicarla por otra cuyos elementos 
    # sean los de la sucesión 1, 1, -1, -1, 1, 1, -1, -1, ...
    # Tal tiene por k-ésimo término a (-1)^(k//2).
    # La lista citada, 'signos' será llamada.

    argumentos_no_negativos = []
    signos = []
    for i in range(1, n + 1):
        sublista_a = []
        sublista_b = []
        for j in range(int((1 + math.sqrt(24*i + 1)) / 6) - 
                       int((1 - math.sqrt(24*i + 1)) / 6)):
            if j % 2 == 0:
                k =  (j + 2) // 2
                sublista_a.append(i - omega(k))
            else:
                k = -(j + 1) // 2
                sublista_a.append(i - omega(k))
            sublista_b.append((-1)**(j // 2))
        argumentos_no_negativos.append(sublista_a)
        signos.append(sublista_b)

    # -------------------------------------------------------------
    m = 1
    while (m <= n):
        # Para poder construir los P[n] en orden creciente, 
        # conviene invertir el orden de cada una de las sublistas  
        # que constituyen las listas declaradas fuera del ciclo.
        # Evitamos así posibles referencias a valores aún no
        # calculados.
        args = argumentos_no_negativos[m-1][::-1]
        sgns = signos[m-1][::-1]

        # ---------------------------------------------------------
        # Finalmente, no resta más que hallar para el valor de 
        # 'm' en cuestión, la suma de todos los productos de la 
        # forma
        # signs[i]*P[args[i], con i = 0, 1, 2, ..., len(args[m]).
        # Dichos productos estarán contenidos en la lista 
        # auxiliar 'elementos_individuales'.

        elementos_individuales = []
        for i in range(len(args)):
            elementos_individuales.append(sgns[i]*P[args[i]])

        P.append(sum(elementos_individuales))
        m += 1

    return P[-1]
  

A000041 = []

for i in range(100):
    A000041.append(p(i))

# Cálculos para especificar el formato de salida.
l = len(str(max(A000041))) + 2
n = 56 // l                    # número de filas.
m = len(A000041) // n       # número de columnas.
if m * n < len(A000041):
    m += 1
M = np.zeros((m, n), dtype = int)
k = 0
for i in range(m):
    for j in range(n):
        if k < len(A000041):
            M[i,j] = A000041[k]
            k += 1
print(M)
  

[[        1         1         2         3         5]
 [        7        11        15        22        30]
 [       42        56        77       101       135]
 [      176       231       297       385       490]
 [      627       792      1002      1255      1575]
 [     1958      2436      3010      3718      4565]
 [     5604      6842      8349     10143     12310]
 [    14883     17977     21637     26015     31185]
 [    37338     44583     53174     63261     75175]
 [    89134    105558    124754    147273    173525]
 [   204226    239943    281589    329931    386155]
 [   451276    526823    614154    715220    831820]
 [   966467   1121505   1300156   1505499   1741630]
 [  2012558   2323520   2679689   3087735   3554345]
 [  4087968   4697205   5392783   6185689   7089500]
 [  8118264   9289091  10619863  12132164  13848650]
 [ 15796476  18004327  20506255  23338469  26543660]
 [ 30167357  34262962  38887673  44108109  49995925]
 [ 56634173  64112359  72533807  82010177  92669720]
 [104651419 118114304 133230930 150198136 169229875]]
  

# Algoritmo de exploración
n = 75000
m = n
limite_sup = 3*n
encontrado = False

while(m <= limite_sup and encontrado == False):
    i = 0
    for elemento in sucesion_de_valores_de_p(n):
        if elemento % 10**6 == 0:
            print(f'i:{i:>5},\np({{i}}) = {{str(elemento)}}')
            encontrado = True
            break
        i += 1
    if m == limite_sup and encontrado == False:
        print(f"No hay ningún valor de n para el que "
              f"1 000 000 | p(n).\n"
              f"Incremente el valor de n o el del límite "
              f"superior de búsqueda.") 
    m += n
  

def mostrar_particiones(l, n, p, v = None, P = None, i = 0):
    '''
        INPUTS:
    
        - "n" es el número de monedas cuyas  particiones deseamos 
          conocer. Se trata, por lo tanto, de un entero positivo.
        - "l", la lista que almacena la sucesión de los enteros 
          positivos menores o iguales que el valor de "n" en la 
          primera llamada de la función.
          Es, consecuentemente, un objeto inmutable.
        - "p" es el valor de p(n); i.e., el número de particiones 
          no restrictivas de "n".
        - "v", el vector que contiene a los escalares a[0], a[1], 
          a[2],..., a[n-1], que satisfacen la ecuación diofántica
          
          a[0] + 2a[1] + 3a[2]... + n·a[n-1] = n
    
          Inicialmente debe ser una lista vacía.
        - "P" es la lista que albergará las posibles 
          representaciones de "n". 
          La conformarán todos los vectores "v" para los que se 
          haya concretado un caso base satisfactorio.
        - "i" es el contador auxiliar empleado para recorrer "l".
            
        OUTPUTS:
            
        - Se devuelve un recuadro que contiene todas las 
          representaciones posibles de "n".
      
    '''
    N = len(l)
    if v is None: v = []
    if P is None: P = []
        
    # CASO BASE SATISFACTORIO: 
    # Este algoritmo recursivo disminuye sucesivamente el valor 
    # primero de 'n', restándole múltiplos de los enteros  
    # contenidos en 'l'. 
    # Si la diferencia resulta nula, habremos dado con una
    # partición de 'n'.
    
    if n == 0:
        # Para expresar la representación recién hallada como 
        # cadena de caracteres (v.g, 5 = 2 + 1 + 1 + 1), debemos 
        # establecer una correspondencia entre las listas 'v' y 
        # 'l' .
        # A título de ejemplo, si fuese n = 5, tras cierto número 
        # de llamadas a la función podríamos obtener la 
        # combinación
        # v = [1, 2, 0, 0, 0]
        # l = [1, 2, 3, 4, 5]
        # Esto supone que, en esta partición, el entero l[0] = 1 
        # puede aparecer 1 vez, pues v[0] = 1, mientras que 
        # l[1] = 2 puede aparecer dos veces. Es decir, que 
        # podemos expresar 5 como 2 + 2 + 1.
        # Nos interesa pues conocer los índices (i. e, las 
        # posiciones) de los elementos no nulos de 'v'.
        # A tal efecto, nos serviremos del método 'nonzero' de 
        # NumPy para identificarlos. 
        # Una vez hecho esto, construiremos una nueva lista en la 
        # que aparecerá repetido cada l[i] no nulos tantas veces 
        # como lo indique v[i]. A tal la llamaremos 'solucion'.
        solucion = []
        for i in range(len(np.nonzero(v)[0])):
            # Llamaremos 'j' al número de veces que aparece el 
            # elemento no nulo  de "l" en la posición 
            # np.nonzero(v)[0][i]
            j = v[np.nonzero(v)[0][i]]
            while(j > 0):
                # En tanto "j" sea positivo, se añadirán copias 
                # del correspondiente valor de 'l'.
                # Insistimos, una copia por iteración.
                solucion.append(l[np.nonzero(v)[0][i]])
                # Tras cada adición, la cuenta disminuye en una 
                # unidad.
                j -= 1
        # Una vez conformada la nueva solución, la dispondremos  
        # en orden decreciente y la incorporamos a la lista de 
        # representaciones, P.
        P.append(solucion[::-1])
        return
        
    # CASOS BASE PRESCINDIBLES:
    # Si las sucesivas sustracciones derivan en un número negativo 
    # o se han efectuado hasta N sustracciones sin llegar al cero, 
    # la combinación hallada es inviable.

    if n < 0 or i == len(l):
        return
    # -------------------------------------------------------------
    
    # Llamaremos 'k_max' al máximo número de veces que l[i] puede 
    # aparecer en una representación de 'n'
    k_max = N // l[i] 
    # Nótese que podríamos haber escrito equivalentemente 
    # l[-1] // l[i], pues tanto N = len(l) como l[-1] se 
    # corresponden con el valor de 'n' pasado en la primera 
    # llamada.

    # Si l = [1,2,3,4,5], entonces los posibles valores de 'k_max' 
    # serán, en virtud del valor del contador 'i' al invocar la 
    # función: 5//1, 5//2, 5//3, 5//4 y 5//5, que, respectivamente, 
    # valen, 5, 2, 1, 1, 1.

    # ALGORITMO RECURSIVO:
    # Las llamadas sucesivas de la función dentro de un ciclo nos 
    # permiten crear listas de posibles soluciones que
    # eventualmente alcanzarán un caso base admisible o uno 
    # descartable.
    # Estas se constituyen incorporando a su estructura al escalar 
    # 'k', que, según la iteración en la que el ciclo se encuentre,
    #  multiplica al elemento l[i].
    nueva_solucion = []
    for k in range(k_max + 1):
        nueva_solucion = v + [k] # Adición de nuevos elementos.
        mostrar_particiones(l, 
                            n - k*l[i], 
                            p, 
                            nueva_solucion, 
                            P, 
                            i + 1)
    
    # -------------------------------------------------------------

    # Construcción de la tabla de salida:
    # Este último proceso inicia en cuanto P tiene tantos  
    # elementos como el valor de p(n).
    if len(P) == p:
        P = sorted(P)
        # Cada sublista_a de P será transformada en una cadena de 
        # caracteres en la que, a fin de conseguir una mayor
        # claridad, se insertarán entre sus elementos signos '+'.
        S = []
        j = 0
        for elemento in P:
            s = str(N) + ' = '
            for i in range(len(elemento)):
                s += str(P[j][i])
                if i < len(elemento) - 1:
                    s += ' + '
            S.append(s)
            j += 1
        tabla = PrettyTable()
        tabla.field_names = ['i', 'i-ésima representación']
        tabla.align['i-ésima representación'] = "l"
        # Finalmente, damos salida a los resultados:
        for i in range(p):
            tabla.add_row([i + 1, S[p - 1 - i]])
        return(tabla)
  

for n in range(1, 11):
    print(mostrar_particiones([i for i in range(1, n + 1)], 
                              n, 
                              p(n), 
                              v = None, 
                              P = None, 
                              i = 0))
  

+---+------------------------+
| i | i-ésima representación |
+---+------------------------+
| 1 | 1 = 1                  |
+---+------------------------+
+---+------------------------+
| i | i-ésima representación |
+---+------------------------+
| 1 | 2 = 2                  |
| 2 | 2 = 1 + 1              |
+---+------------------------+
+---+------------------------+
| i | i-ésima representación |
+---+------------------------+
| 1 | 3 = 3                  |
| 2 | 3 = 2 + 1              |
| 3 | 3 = 1 + 1 + 1          |
+---+------------------------+
+---+------------------------+
| i | i-ésima representación |
+---+------------------------+
| 1 | 4 = 4                  |
| 2 | 4 = 3 + 1              |
| 3 | 4 = 2 + 2              |
| 4 | 4 = 2 + 1 + 1          |
| 5 | 4 = 1 + 1 + 1 + 1      |
+---+------------------------+
+---+------------------------+
| i | i-ésima representación |
+---+------------------------+
| 1 | 5 = 5                  |
| 2 | 5 = 4 + 1              |
| 3 | 5 = 3 + 2              |
| 4 | 5 = 3 + 1 + 1          |
| 5 | 5 = 2 + 2 + 1          |
| 6 | 5 = 2 + 1 + 1 + 1      |
| 7 | 5 = 1 + 1 + 1 + 1 + 1  |
+---+------------------------+
+----+---------------------------+
| i  | i-ésima representación    |
+----+---------------------------+
| 1  | 6 = 6                     |
| 2  | 6 = 5 + 1                 |
| 3  | 6 = 4 + 2                 |
| 4  | 6 = 4 + 1 + 1             |
| 5  | 6 = 3 + 3                 |
| 6  | 6 = 3 + 2 + 1             |
| 7  | 6 = 3 + 1 + 1 + 1         |
| 8  | 6 = 2 + 2 + 2             |
| 9  | 6 = 2 + 2 + 1 + 1         |
| 10 | 6 = 2 + 1 + 1 + 1 + 1     |
| 11 | 6 = 1 + 1 + 1 + 1 + 1 + 1 |
+----+---------------------------+
+----+-------------------------------+
| i  | i-ésima representación        |
+----+-------------------------------+
| 1  | 7 = 7                         |
| 2  | 7 = 6 + 1                     |
| 3  | 7 = 5 + 2                     |
| 4  | 7 = 5 + 1 + 1                 |
| 5  | 7 = 4 + 3                     |
| 6  | 7 = 4 + 2 + 1                 |
| 7  | 7 = 4 + 1 + 1 + 1             |
| 8  | 7 = 3 + 3 + 1                 |
| 9  | 7 = 3 + 2 + 2                 |
| 10 | 7 = 3 + 2 + 1 + 1             |
| 11 | 7 = 3 + 1 + 1 + 1 + 1         |
| 12 | 7 = 2 + 2 + 2 + 1             |
| 13 | 7 = 2 + 2 + 1 + 1 + 1         |
| 14 | 7 = 2 + 1 + 1 + 1 + 1 + 1     |
| 15 | 7 = 1 + 1 + 1 + 1 + 1 + 1 + 1 |
+----+-------------------------------+
+----+-----------------------------------+
| i  | i-ésima representación            |
+----+-----------------------------------+
| 1  | 8 = 8                             |
| 2  | 8 = 7 + 1                         |
| 3  | 8 = 6 + 2                         |
| 4  | 8 = 6 + 1 + 1                     |
| 5  | 8 = 5 + 3                         |
| 6  | 8 = 5 + 2 + 1                     |
| 7  | 8 = 5 + 1 + 1 + 1                 |
| 8  | 8 = 4 + 4                         |
| 9  | 8 = 4 + 3 + 1                     |
| 10 | 8 = 4 + 2 + 2                     |
| 11 | 8 = 4 + 2 + 1 + 1                 |
| 12 | 8 = 4 + 1 + 1 + 1 + 1             |
| 13 | 8 = 3 + 3 + 2                     |
| 14 | 8 = 3 + 3 + 1 + 1                 |
| 15 | 8 = 3 + 2 + 2 + 1                 |
| 16 | 8 = 3 + 2 + 1 + 1 + 1             |
| 17 | 8 = 3 + 1 + 1 + 1 + 1 + 1         |
| 18 | 8 = 2 + 2 + 2 + 2                 |
| 19 | 8 = 2 + 2 + 2 + 1 + 1             |
| 20 | 8 = 2 + 2 + 1 + 1 + 1 + 1         |
| 21 | 8 = 2 + 1 + 1 + 1 + 1 + 1 + 1     |
| 22 | 8 = 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 |
+----+-----------------------------------+
+----+---------------------------------------+
| i  | i-ésima representación                |
+----+---------------------------------------+
| 1  | 9 = 9                                 |
| 2  | 9 = 8 + 1                             |
| 3  | 9 = 7 + 2                             |
| 4  | 9 = 7 + 1 + 1                         |
| 5  | 9 = 6 + 3                             |
| 6  | 9 = 6 + 2 + 1                         |
| 7  | 9 = 6 + 1 + 1 + 1                     |
| 8  | 9 = 5 + 4                             |
| 9  | 9 = 5 + 3 + 1                         |
| 10 | 9 = 5 + 2 + 2                         |
| 11 | 9 = 5 + 2 + 1 + 1                     |
| 12 | 9 = 5 + 1 + 1 + 1 + 1                 |
| 13 | 9 = 4 + 4 + 1                         |
| 14 | 9 = 4 + 3 + 2                         |
| 15 | 9 = 4 + 3 + 1 + 1                     |
| 16 | 9 = 4 + 2 + 2 + 1                     |
| 17 | 9 = 4 + 2 + 1 + 1 + 1                 |
| 18 | 9 = 4 + 1 + 1 + 1 + 1 + 1             |
| 19 | 9 = 3 + 3 + 3                         |
| 20 | 9 = 3 + 3 + 2 + 1                     |
| 21 | 9 = 3 + 3 + 1 + 1 + 1                 |
| 22 | 9 = 3 + 2 + 2 + 2                     |
| 23 | 9 = 3 + 2 + 2 + 1 + 1                 |
| 24 | 9 = 3 + 2 + 1 + 1 + 1 + 1             |
| 25 | 9 = 3 + 1 + 1 + 1 + 1 + 1 + 1         |
| 26 | 9 = 2 + 2 + 2 + 2 + 1                 |
| 27 | 9 = 2 + 2 + 2 + 1 + 1 + 1             |
| 28 | 9 = 2 + 2 + 1 + 1 + 1 + 1 + 1         |
| 29 | 9 = 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1     |
| 30 | 9 = 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 |
+----+---------------------------------------+
+----+--------------------------------------------+
| i  | i-ésima representación                     |
+----+--------------------------------------------+
| 1  | 10 = 10                                    |
| 2  | 10 = 9 + 1                                 |
| 3  | 10 = 8 + 2                                 |
| 4  | 10 = 8 + 1 + 1                             |
| 5  | 10 = 7 + 3                                 |
| 6  | 10 = 7 + 2 + 1                             |
| 7  | 10 = 7 + 1 + 1 + 1                         |
| 8  | 10 = 6 + 4                                 |
| 9  | 10 = 6 + 3 + 1                             |
| 10 | 10 = 6 + 2 + 2                             |
| 11 | 10 = 6 + 2 + 1 + 1                         |
| 12 | 10 = 6 + 1 + 1 + 1 + 1                     |
| 13 | 10 = 5 + 5                                 |
| 14 | 10 = 5 + 4 + 1                             |
| 15 | 10 = 5 + 3 + 2                             |
| 16 | 10 = 5 + 3 + 1 + 1                         |
| 17 | 10 = 5 + 2 + 2 + 1                         |
| 18 | 10 = 5 + 2 + 1 + 1 + 1                     |
| 19 | 10 = 5 + 1 + 1 + 1 + 1 + 1                 |
| 20 | 10 = 4 + 4 + 2                             |
| 21 | 10 = 4 + 4 + 1 + 1                         |
| 22 | 10 = 4 + 3 + 3                             |
| 23 | 10 = 4 + 3 + 2 + 1                         |
| 24 | 10 = 4 + 3 + 1 + 1 + 1                     |
| 25 | 10 = 4 + 2 + 2 + 2                         |
| 26 | 10 = 4 + 2 + 2 + 1 + 1                     |
| 27 | 10 = 4 + 2 + 1 + 1 + 1 + 1                 |
| 28 | 10 = 4 + 1 + 1 + 1 + 1 + 1 + 1             |
| 29 | 10 = 3 + 3 + 3 + 1                         |
| 30 | 10 = 3 + 3 + 2 + 2                         |
| 31 | 10 = 3 + 3 + 2 + 1 + 1                     |
| 32 | 10 = 3 + 3 + 1 + 1 + 1 + 1                 |
| 33 | 10 = 3 + 2 + 2 + 2 + 1                     |
| 34 | 10 = 3 + 2 + 2 + 1 + 1 + 1                 |
| 35 | 10 = 3 + 2 + 1 + 1 + 1 + 1 + 1             |
| 36 | 10 = 3 + 1 + 1 + 1 + 1 + 1 + 1 + 1         |
| 37 | 10 = 2 + 2 + 2 + 2 + 2                     |
| 38 | 10 = 2 + 2 + 2 + 2 + 1 + 1                 |
| 39 | 10 = 2 + 2 + 2 + 1 + 1 + 1 + 1             |
| 40 | 10 = 2 + 2 + 1 + 1 + 1 + 1 + 1 + 1         |
| 41 | 10 = 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1     |
| 42 | 10 = 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 |
+----+--------------------------------------------+