He calculado el primer n=22 coeficientes de dominó basados en las ideas de Vepir. El algoritmo tiene una complejidad de tiempo O(4n) y la complejidad espacial O(2n) (cálculo para n=22 me llevó alrededor de 5 horas), porque no almacenamos la matriz M y calcular sus valores siempre que sea necesario (la comprobación se realiza en la sentencia if, tiene una forma sencilla gracias a la implementación de operaciones bitwise). El código está en C++ y estoy usando boost para poder usar tipos de datos más grandes en lugar de sólo 64 bit long int (de hecho, 256 bits es el tamaño más pequeño que se puede utilizar para calcular el primer 22 coeficientes e incluso se agotará pronto - D22≈1.04159∗1053 por lo que la secuencia efectivamente crece rápidamente).
#include <boost/multiprecision/cpp_int.hpp>
#include <iostream>
#include <vector>
int main() {
const int n = 22;
std::vector<boost::multiprecision::int256_t> v(1 << n), w = v;
v[0] = 1;
v[1] = 1;
std::cout << "0: " << 1 << "\n1: " << 2 << "\n";
for (int k = 2; k <= n; k++) {
for (int j = 0; j < 1 << k; j++) {
w[j] = 0;
for (int i = 0; i < 1 << (k - 1); i++)
if (!(j & ~i & ~(i << 1)))
w[j] += v[i];
}
boost::multiprecision::int256_t s = 0;
for (int j = 0; j < 1 << k; j++) {
v[j] = w[j];
s += v[j];
}
std::cout << k << ": " << s << "\n";
}
return 0;
}
Los coeficientes resultantes son:
D0=1
D1=2
D2=5
D3=18
D4=97
D5=802
D6=10565
D7=228850
D8=8289217
D9=506526530
D10=52501381765
D11=9260170733266
D12=2784551512218145
D13=1429063630276963426
D14=1252517782851235507141
D15=1875484239084442842046130
D16=4798818821638537354534159233
D17=20984654018757393270224583817858
D18=156836820191186149182649796402611461
D19=2003511561206789377122778581243926768018
D20=43746600983206157305472349537031148375171041
D21=1632720578386908900755016228846582339781524154786
D22=104159355027420388820987121345303301272029996166458949
D23=11358102129403917231432165494822142898736362479589146061682
D24=2117071962387596476689143406365411645227548761608194942194112321
D25=674510308152303132671249373016113929650762465949837199637457187249602
D26=367337379241994889203376240956745359336562483947308610993230043730270521221
El logaritmo de Dn parece crecer de forma cuadrática (lo que se puede apreciar incluso mirando la forma de la lista anterior). Tenemos aproximadamente: Dn≈20.386712n2−0.667456n+3.63662=exp(0.268048n2−0.462645n+2.52072) Pero no sé qué precisión tiene esta aproximación ni cuáles podrían ser los coeficientes reales. Y tampoco dice mucho sobre la posible fórmula recursiva: parece que sigue estando fuera de alcance.
Este problema es un caso especial del más general - dado un grafo dirigido G Cuántas formas hay de asignar números 0 et 1 a sus vértices de forma que las flechas siempre apunten en la dirección no decreciente (es decir, no de 1 a 0 )? Podemos suponer, sin pérdida de generalidad, que el grafo es conexo (si no, la respuesta es sólo un producto de las respuestas de los componentes) y que es acíclico (el ciclo implica que todos sus tienen que tener el mismo número, por lo que podemos colapsarlos en un vértice). Hay un algoritmo fácil para responder a este problema en el caso de un árbol dirigido (sólo la recursión de las hojas hacia arriba), pero con vértices que interactúan de una manera más compleja, como en este problema, parece difícil, si no imposible, para llegar a un algoritmo rápido, sobre todo teniendo en cuenta, lo rápido que estos coeficientes crecen (aunque, por ejemplo, la fórmula para los árboles binarios perfectos de profundidad n crece aún más rápido y tiene una simple recursión Bn+1=B2n+1 con B0=1 ).
EDITAR
La primera 5 vectores vi generando estos Di son (se omiten los ceros finales):
v0=(0,1)
v1=(1,1)
v2=(2,1,1,1)
v3=(5,2,3,2,2,1,2,1)
v4=(18,6,11,6,11,4,9,4,6,2,4,2,6,2,4,2)
v5=(97,28,56,28,65,20,48,20,56,16,33,16,48,14,31,14,28,8,16,8,20,6,14,6,28,8,16,8,20,6,14,6)
Podemos encontrar la interpretación de cada coeficiente en este vector basándonos en las propiedades fractales de M - por ejemplo el primer elemento del vector v es sólo la suma de todos los elementos del vector anterior (y por lo tanto es sólo Dn−1 . El segundo elemento es la suma de cada elemento en la posición impar (si contamos las posiciones desde la izquierda empezando por 0 ) y la tercera es la suma en cada posición no divisible por 4 . Sin embargo, a medida que avanzamos los patrones se vuelven más complejos, pero analizándolos se puede obtener tal vez un algoritmo más rápido (por ejemplo en lugar de sumar 34 de los elementos del vector anterior para obtener el tercer elemento se puede simplemente restar del primer elemento previamente calculado la suma de los 14 elementos con índices divisibles por 4 aunque esta mejora puede producir como máximo el doble de velocidad).
EDITAR 2
Pude calcular más términos de la secuencia (hasta 26 en lugar de 22 ) utilizando el siguiente algoritmo:
#include <boost/multiprecision/cpp_int.hpp>
#include <iostream>
#include <vector>
int main() {
int k = 0, n = 26;
std::vector<boost::multiprecision::int256_t> v(1 << n - 2, 1), a(n + 1, 1);
for (int i = 2; i < a.size(); i++)
a[i] = 3 * a[i - 1] + 2 * a[i - 2];
for (int i = 0; i < v.size(); i++)
for (int j = i | (i << 1); j; j >>= 1) {
k += j & 1;
if (!(j & 1 && j - 1)){
v[i] *= a[k + 1];
k = 0;
}
}
for (int i = n - 3; i > 0; i--) {
std::vector<boost::multiprecision::int256_t> w(1 << i, 0);
for (int j = 0; j < 1 << i; j++)
for (int k = 0; k < 1 << i + 1; k++)
if (!(~(j | j << 1) & k))
w[j] += v[k];
v = w;
std::cout << n - i + 1 << " " << v[1] + 1 << "\n";
}
}
Lo que puede ser fácilmente paralelizado acelerando el proceso aún más. La forma de calcular los números también es diferente. Dejemos que D(s,n) denotan la configuración válida del dominó con la cadena binaria s en la parte superior y teniendo n capas - el problema inicial es evaluar D(′0′,n)+D(′1′,n)=D(′1′,n)+1 . Para obtener D(s,n) observamos que para n=1 todos son iguales a uno y para mayor n podemos considerar todas las cadenas s′ directamente debajo n , calcula recursivamente todos los valores D(s′,n−1) y sumarlos. La complejidad asintótica es la misma, pero podemos acelerar el programa partiendo de D(s,2) que es igual a 2 a la potencia del número de unos en s o (como en el algoritmo anterior) de D(s,3) observando que el número D(sk,3) donde sk es la cadena formada por k es igual a la OEIS A007483 que tiene una recursión simple y que para todos los demás s es sólo el producto de estos números para cada bloque consecutivo de unos en s . ak=D(sk,4) , bk=D(sk,5) et ck=(sk,6) también tienen sus recurrencias: ak+3=4ak+2+12ak+1+8ak bk+5=5bk+4+48bk+3+108bk+2+64bk+1−32bk ck+8=6ck+7+164ck+6+896ck+5+1584ck+4−992ck+3−3328ck+2+3584ck+1−1024ck Con a0=b0=c0=1 , a1=18=D3 , a2=96=D4−1 , b1=97=D4 , b2=801=D5−1 , b3=8961 , b4=93697 , c1=802=D5 , c2=10564=D6−1 , c3=207720 , c4=3692192 , c5=66924256 , c6=1209149120 , c7=21854550528 . Para D(sk,7) No he encontrado tal serie - o no existe o (más probablemente) no hay suficientes términos para detectarla. Tampoco he encontrado una forma de calcular D(s,n) para cualquier n≥4 para todas las cadenas s por lo que el algoritmo sólo utiliza D(s,3) - si uno es capaz de calcularlos directamente puede obtener términos más grandes de la secuencia. Sólo quiero mencionar que D26≈2247.7 es el término más grande por debajo de 2255 por lo que habrá que utilizar tipos enteros más grandes cuando se intente ir más allá.