Este desafío es en parte un desafío de algoritmos, involucra algunas matemáticas y es en parte simplemente un desafío de código más rápido.

Para algún entero positivo n, considere una cadena aleatoria de forma uniforme 1s y 0s de longitud ny llamarlo A. Ahora también considere una segunda cadena aleatoria de longitud elegida uniformemente ncuyos valores son -1, 0,o 1y llámela B_pre. Ahora deja Bser B_pre+ B_pre. Eso se B_preconcatena a sí mismo.

Consideremos ahora el producto interno de Ay B[j,...,j+n-1]y llamarlo Z_jy el índice de 1.

Tarea

El resultado debe ser una lista de n+1fracciones. El iésimo término en la salida debe ser la exacta probabilidad de que todo de los primeros itérminos Z_jcon j <= iigual 0.

Puntuación

El más grande npara el que su código proporciona la salida correcta en menos de 10 minutos en mi máquina.

Desempate

Si dos respuestas tienen el mismo puntaje, gana la que presentó primero.

En el caso (muy, muy poco probable) de que alguien encuentre un método para obtener puntajes ilimitados, se aceptará la primera prueba válida de tal solución.

Insinuación

No intentes resolver este problema matemáticamente, es demasiado difícil. En mi opinión, la mejor manera es volver a las definiciones básicas de probabilidad de la escuela secundaria y encontrar formas inteligentes de obtener el código para hacer una enumeración exhaustiva de las posibilidades.

Idiomas y bibliotecas

Puede usar cualquier idioma que tenga un compilador / intérprete / etc. para Linux y cualquier biblioteca que también esté disponible gratuitamente para Linux.

Mi máquina Los tiempos se ejecutarán en mi máquina. Esta es una instalación estándar de ubuntu en un procesador AMD FX-8350 de ocho núcleos. Esto también significa que necesito poder ejecutar su código. Como consecuencia, solo use software gratuito fácilmente disponible e incluya instrucciones completas sobre cómo compilar y ejecutar su código.

Algunas salidas de prueba. Considere solo el primer resultado para cada uno n. Eso es cuando i=1. Para nde 1 a 13 que deberían ser.

 1: 4/6
 2: 18/36
 3: 88/216
 4: 454/1296
 5: 2424/7776
 6: 13236/46656
 7: 73392/279936
 8: 411462/1679616
 9: 2325976/10077696
10: 13233628/60466176
11: 75682512/362797056
12: 434662684/2176782336
13: 2505229744/13060694016

También puede encontrar una fórmula general i=1en http://oeis.org/A081671 .

Tabla de clasificación (dividida por idioma)

• n = 15. Python + python paralelo + pypy en 1min49s por Jakube
• n = 17. C ++ en 3min37s por Keith Randall
• n = 16. C ++ en 2min38s por kuroi neko

C ++, n = 18 en 9 min en 8 hilos

(Avíseme si dura menos de 10 minutos en su máquina).

Aprovecho varias formas de simetría en la matriz B. Esos son cíclicos (desplazamiento en una posición), inversión (invertir el orden de los elementos) y signo (tomar el negativo de cada elemento). Primero calculo la lista de Bs que necesitamos probar y su peso. Luego, cada B se ejecuta a través de una rutina rápida (usando instrucciones de conteo de bits) para todos los 2 ^ n valores de A.

Aquí está el resultado para n == 18:

> time ./a.out 18
 1: 16547996212044 / 101559956668416
 2:  3120508430672 / 101559956668416
 3:   620923097438 / 101559956668416
 4:   129930911672 / 101559956668416
 5:    28197139994 / 101559956668416
 6:     6609438092 / 101559956668416
 7:     1873841888 / 101559956668416
 8:      813806426 / 101559956668416
 9:      569051084 / 101559956668416
10:      510821156 / 101559956668416
11:      496652384 / 101559956668416
12:      493092812 / 101559956668416
13:      492186008 / 101559956668416
14:      491947940 / 101559956668416
15:      491889008 / 101559956668416
16:      449710584 / 101559956668416
17:      418254922 / 101559956668416
18:      409373626 / 101559956668416

real    8m55.854s
user    67m58.336s
sys 0m5.607s

Compile el siguiente programa con g++ --std=c++11 -O3 -mpopcnt dot.cc

#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <thread>
#include <mutex>
#include <chrono>

using namespace std;

typedef long long word;

word n;

void inner(word bpos, word bneg, word w, word *cnt) {
    word maxi = n-1;
    for(word a = (1<<n)-1; a >= 0; a--) {
        word m = a;
        for(word i = maxi; i >= 0; i--, m <<= 1) {
            if(__builtin_popcount(m&bpos) != __builtin_popcount(m&bneg))
                break;
            cnt[i]+=w;
        }
    }
}

word pow(word n, word e) {
    word r = 1;
    for(word i = 0; i < e; i++) r *= n;
    return r;
}

typedef struct {
    word b;
    word weight;
} Bentry;

mutex block;
Bentry *bqueue;
word bhead;
word btail;
word done = -1;

word maxb;

// compute -1*b
word bneg(word b) {
    word w = 1;
    for(word i = 0; i < n; i++, w *= 3) {
        word d = b / w % 3;
        if(d == 1)
            b += w;
        if(d == 2)
            b -= w;
    }
    return b;
}

// rotate b one position
word brot(word b) {
    b *= 3;
    b += b / maxb;
    b %= maxb;
    return b;
}

// reverse b
word brev(word b) {
    word r = 0;
    for(word i = 0; i < n; i++) {
        r *= 3;
        r += b % 3;
        b /= 3;
    }
    return r;
}

// individual thread's work routine
void work(word *cnt) {
    while(true) {
        // get a queue entry to work on
        block.lock();
        if(btail == done) {
            block.unlock();
            return;
        }
        if(bhead == btail) {
            block.unlock();
            this_thread::sleep_for(chrono::microseconds(10));
            continue;
        }
        word i = btail++;
        block.unlock();

        // thread now owns bqueue[i], work on it
        word b = bqueue[i].b;
        word w = 1;
        word bpos = 0;
        word bneg = 0;
        for(word j = 0; j < n; j++, b /= 3) {
            word d = b % 3;
            if(d == 1)
                bpos |= 1 << j;
            if(d == 2)
                bneg |= 1 << j;
        }
        bpos |= bpos << n;
        bneg |= bneg << n;
        inner(bpos, bneg, bqueue[i].weight, cnt);
    }
}

int main(int argc, char *argv[]) {
    n = atoi(argv[1]);

    // allocate work queue
    maxb = pow(3, n);
    bqueue = (Bentry*)(malloc(maxb*sizeof(Bentry)));

    // start worker threads
    word procs = thread::hardware_concurrency();
    vector<thread> threads;
    vector<word*> counts;
    for(word p = 0; p < procs; p++) {
        word *cnt = (word*)calloc(64+n*sizeof(word), 1);
        threads.push_back(thread(work, cnt));
        counts.push_back(cnt);
    }

    // figure out which Bs we actually want to test, and with which weights
    bool *bmark = (bool*)calloc(maxb, 1);
    for(word i = 0; i < maxb; i++) {
        if(bmark[i]) continue;
        word b = i;
        word w = 0;
        for(word j = 0; j < 2; j++) {
            for(word k = 0; k < 2; k++) {
                for(word l = 0; l < n; l++) {
                    if(!bmark[b]) {
                        bmark[b] = true;
                        w++;
                    }
                    b = brot(b);
                }
                b = bneg(b);
            }
            b = brev(b);
        }
        bqueue[bhead].b = i;
        bqueue[bhead].weight = w;
        block.lock();
        bhead++;
        block.unlock();
    }
    block.lock();
    done = bhead;
    block.unlock();

    // add up results from threads
    word *cnt = (word*)calloc(n,sizeof(word));
    for(word p = 0; p < procs; p++) {
        threads[p].join();
        for(int i = 0; i < n; i++) cnt[i] += counts[p][i];
    }
    for(word i = 0; i < n; i++)
        printf("%2lld: %14lld / %14lld\n", i+1, cnt[n-1-i], maxb<<n);
    return 0;
}

Python 2 usando pypy y pp: n = 15 en 3 minutos

También solo una simple fuerza bruta. Es interesante ver que casi obtengo la misma velocidad que Kuroi Neko con C ++. Mi código puede alcanzarn = 12 en unos 5 minutos. Y solo lo ejecuto en un núcleo virtual.

editar: reducir el espacio de búsqueda por un factor de n

Noté que un vector cíclico A*de Aproduce los mismos números que las probabilidades (los mismos números) que el vector original Acuando repito B. Ej El vector (1, 1, 0, 1, 0, 0)tiene las mismas probabilidades como cada uno de los vectores (1, 0, 1, 0, 0, 1), (0, 1, 0, 0, 1, 1), (1, 0, 0, 1, 1, 0), (0, 0, 1, 1, 0, 1)y (0, 1, 1, 0, 1, 0)al elegir una al azar B. Por lo tanto, no tengo que iterar sobre cada uno de estos 6 vectores, sino solo alrededor de 1 y reemplazar count[i] += 1con count[i] += cycle_number.

Esto reduce la complejidad de Theta(n) = 6^na Theta(n) = 6^n / n. Por n = 13lo tanto , es aproximadamente 13 veces más rápido que mi versión anterior. Se calcula n = 13en aproximadamente 2 minutos y 20 segundos. porn = 14 todavía es un poco lento. Tarda unos 13 minutos.

editar 2: programación multi-core

No estoy muy contento con la próxima mejora. Decidí también intentar ejecutar mi programa en múltiples núcleos. En mis núcleos 2 + 2 ahora puedo calcular n = 14en aproximadamente 7 minutos. Solo un factor de mejora 2.

El código está disponible en este repositorio de github: Link . La programación multinúcleo es un poco fea.

editar 3: Reducción del espacio de búsqueda para Avectores y Bvectores

Noté la misma simetría de espejo para los vectores Aque Kuroi Neko. Todavía no estoy seguro, por qué esto funciona (y si funciona para cada uno n).

La reducción del espacio de búsqueda de Bvectores es un poco más inteligente. Reemplacé la generación de los vectores ( itertools.product), con una función propia. Básicamente, comienzo con una lista vacía y la pongo en una pila. Hasta que la pila esté vacía, elimino una lista, si no tiene la misma longitud n, genero otras 3 listas (agregando -1, 0, 1) y empujándolas a la pila. Si una lista tiene la misma longitud que n, puedo evaluar las sumas.

Ahora que genero los vectores yo mismo, puedo filtrarlos dependiendo de si puedo alcanzar la suma = 0 o no. Por ejemplo, si mi vector Aes (1, 1, 1, 0, 0), y mi vector se Bve (1, 1, ?, ?, ?), lo sé, que no puedo llenar los ?valores, de modo que A*B = 0. Así que no tengo que iterar sobre todos esos 6 vectores Bde la forma(1, 1, ?, ?, ?) .

Podemos mejorar esto si ignoramos los valores para 1. Como se señaló en la pregunta, los valores para i = 1son la secuencia A081671 . Hay muchas formas de calcularlos. Elijo la repetición simple: a(n) = (4*(2*n-1)*a(n-1) - 12*(n-1)*a(n-2)) / n. Como i = 1básicamente podemos calcular en poco tiempo, podemos filtrar más vectores para B. Por ejemplo A = (0, 1, 0, 1, 1)y B = (1, -1, ?, ?, ?). Podemos ignorar los vectores, donde los primeros ? = 1, porque los A * cycled(B) > 0, para todos estos vectores. Espero que puedas seguir. Probablemente no sea el mejor ejemplo.

Con esto puedo calcular n = 15en 6 minutos.

editar 4:

Implementé rápidamente la gran idea de kuroi neko, que dice eso By -Bproduce los mismos resultados. Aceleración x2. Sin embargo, la implementación es solo un truco rápido. n = 15en 3 minutos

Código:

Para ver el código completo, visite Github . El siguiente código es solo una representación de las características principales. Omití importaciones, programación multinúcleo, impresión de resultados, ...

count = [0] * n
count[0] = oeis_A081671(n)

#generating all important vector A
visited = set(); todo = dict()
for A in product((0, 1), repeat=n):
    if A not in visited:
        # generate all vectors, which have the same probability
        # mirrored and cycled vectors
        same_probability_set = set()
        for i in range(n):
            tmp = [A[(i+j) % n] for j in range(n)]
            same_probability_set.add(tuple(tmp))
            same_probability_set.add(tuple(tmp[::-1]))
        visited.update(same_probability_set)
        todo[A] = len(same_probability_set)

# for each vector A, create all possible vectors B
stack = []
for A, cycled_count in dict_A.iteritems():
    ones = [sum(A[i:]) for i in range(n)] + [0]
    # + [0], so that later ones[n] doesn't throw a exception
    stack.append(([0] * n, 0, 0, 0, False))

    while stack:
        B, index, sum1, sum2, used_negative = stack.pop()
        if index < n:
            # fill vector B[index] in all possible ways,
            # so that it's still possible to reach 0.
            if used_negative:
                for v in (-1, 0, 1):
                    sum1_new = sum1 + v * A[index]
                    sum2_new = sum2 + v * A[index - 1 if index else n - 1]
                    if abs(sum1_new) <= ones[index+1]:
                        if abs(sum2_new) <= ones[index] - A[n-1]:
                            C = B[:]
                            C[index] = v
                            stack.append((C, index + 1, sum1_new, sum2_new, True))
            else:
                for v in (0, 1):
                    sum1_new = sum1 + v * A[index]
                    sum2_new = sum2 + v * A[index - 1 if index else n - 1]
                    if abs(sum1_new) <= ones[index+1]:
                        if abs(sum2_new) <= ones[index] - A[n-1]:
                            C = B[:]
                            C[index] = v
                            stack.append((C, index + 1, sum1_new, sum2_new, v == 1))
        else:
            # B is complete, calculate the sums
            count[1] += cycled_count  # we know that the sum = 0 for i = 1
            for i in range(2, n):
                sum_prod = 0
                for j in range(n-i):
                    sum_prod += A[j] * B[i+j]
                for j in range(i):
                    sum_prod += A[n-i+j] * B[j]
                if sum_prod:
                    break
                else:
                    if used_negative:
                        count[i] += 2*cycled_count
                    else:
                        count[i] += cycled_count

Uso:

Tienes que instalar pypy (para Python 2 !!!). El módulo paralelo de Python no está portado para Python 3. Luego debe instalar el módulo paralelo de Python pp-1.6.4.zip . Extraerlo cden la carpeta y llamar pypy setup.py install.

Entonces puedes llamar a mi programa con

pypy you-do-the-math.py 15

Determinará automáticamente el número de CPU. Puede haber algunos mensajes de error después de finalizar el programa, simplemente ignórelos. n = 16debería ser posible en su máquina.

Salida:

Calculation for n = 15 took 2:50 minutes

 1  83940771168 / 470184984576  17.85%
 2  17379109692 / 470184984576   3.70%
 3   3805906050 / 470184984576   0.81%
 4    887959110 / 470184984576   0.19%
 5    223260870 / 470184984576   0.05%
 6     67664580 / 470184984576   0.01%
 7     30019950 / 470184984576   0.01%
 8     20720730 / 470184984576   0.00%
 9     18352740 / 470184984576   0.00%
10     17730480 / 470184984576   0.00%
11     17566920 / 470184984576   0.00%
12     17521470 / 470184984576   0.00%
13     17510280 / 470184984576   0.00%
14     17507100 / 470184984576   0.00%
15     17506680 / 470184984576   0.00%

Notas e ideas:

• Tengo un procesador i7-4600m con 2 núcleos y 4 hilos. No importa si uso 2 o 4 hilos. El uso de la CPU es del 50% con 2 subprocesos y del 100% con 4 subprocesos, pero aún requiere la misma cantidad de tiempo. No se porque. Verifiqué que cada subproceso solo tiene la mitad de la cantidad de datos, cuando hay 4 subprocesos, verifiqué los resultados, ...
• Yo uso muchas listas. Python no es muy eficiente en el almacenamiento, tengo que copiar muchas listas, así que pensé en usar un número entero. Podría usar los bits 00 (para 0) y 11 (para 1) en el vector A, y los bits 10 (para -1), 00 (para 0) y 01 (para 1) en el vector B. Para el producto de A y B, solo tendría que calcular A & By contar los bloques 01 y 10. El ciclismo se puede hacer cambiando el vector y usando máscaras, ... De hecho, implementé todo esto, puedes encontrarlo en algunos de mis compromisos más antiguos en Github. Pero resultó ser más lento que con las listas. Supongo que pypy realmente optimiza las operaciones de la lista.

matón lanudo - C ++ - demasiado lento

Bueno, ya que un mejor programador asumió la implementación de C ++, estoy llamando a dejar de fumar para este.

#include <cstdlib>
#include <cmath>
#include <vector>
#include <bitset>
#include <future>
#include <iostream>
#include <iomanip>

using namespace std;

/*
6^^n events will be generated, so the absolute max
that can be counted by a b bits integer is
E(b*log(2)/log(6)), i.e. n=24 for a 64 bits counter

To enumerate 3 possible values of a size n vector we need
E(n*log(3)/log(2))+1 bits, i.e. 39 bits
*/
typedef unsigned long long Counter; // counts up to 6^^24

typedef unsigned long long Benumerator; // 39 bits
typedef unsigned long      Aenumerator; // 24 bits

#define log2_over_log6 0.3869

#define A_LENGTH ((size_t)(8*sizeof(Counter)*log2_over_log6))
#define B_LENGTH (2*A_LENGTH)

typedef bitset<B_LENGTH> vectorB;

typedef vector<Counter> OccurenceCounters;

// -----------------------------------------------------------------
// multithreading junk for CPUs detection and allocation
// -----------------------------------------------------------------
int number_of_CPUs(void)
{
    int res = thread::hardware_concurrency();
    return res == 0 ? 8 : res;
}

#ifdef __linux__
#include <sched.h>
void lock_on_CPU(int cpu)
{
    cpu_set_t mask;
    CPU_ZERO(&mask);
    CPU_SET(cpu, &mask);
    sched_setaffinity(0, sizeof(mask), &mask);
}
#elif defined (_WIN32)
#include <Windows.h>
#define lock_on_CPU(cpu) SetThreadAffinityMask(GetCurrentThread(), 1 << cpu)
#else
// #warning is not really standard, so this might still cause compiler errors on some platforms. Sorry about that.
#warning "Thread processor affinity settings not supported. Performances might be improved by providing a suitable alternative for your platform"
#define lock_on_CPU(cpu)
#endif

// -----------------------------------------------------------------
// B values generator
// -----------------------------------------------------------------
struct Bvalue {
    vectorB p1;
    vectorB m1;
};

struct Bgenerator {
    int n;                 // A length
    Aenumerator stop;      // computation limit
    Aenumerator zeroes;    // current zeroes pattern
    Aenumerator plusminus; // current +1/-1 pattern
    Aenumerator pm_limit;  // upper bound of +1/-1 pattern

    Bgenerator(int n, Aenumerator start=0, Aenumerator stop=0) : n(n), stop(stop)
    {
        // initialize generator so that first call to next() will generate first value
        zeroes    = start - 1;
        plusminus = -1;
        pm_limit  = 0;
    }

    // compute current B value
    Bvalue value(void)
    {
        Bvalue res;
        Aenumerator pm = plusminus;
        Aenumerator position = 1;
        int i_pm = 0;
        for (int i = 0; i != n; i++)
        {
            if (zeroes & position)
            {
                if (i_pm == 0)  res.p1 |= position; // first non-zero value fixed to +1
                else         
                {
                    if (pm & 1) res.m1 |= position; // next non-zero values
                    else        res.p1 |= position;
                    pm >>= 1;
                }
                i_pm++;
            }
            position <<= 1;
        }
        res.p1 |= (res.p1 << n); // concatenate 2 Bpre instances
        res.m1 |= (res.m1 << n);
        return res;
    }

    // next value
    bool next(void)
    {
        if (++plusminus == pm_limit)
        {
            if (++zeroes == stop) return false;
            plusminus = 0;
            pm_limit = (1 << vectorB(zeroes).count()) >> 1;
        }
        return true;
    }

    // calibration: produces ranges that will yield the approximate same number of B values
    vector<Aenumerator> calibrate(int segments)
    {
        // setup generator for the whole B range
        zeroes = 0;
        stop = 1 << n;
        plusminus = -1;
        pm_limit = 0;

        // divide range into (nearly) equal chunks
        Aenumerator chunk_size = ((Aenumerator)pow (3,n)-1) / 2 / segments;

        // generate bounds for zeroes values
        vector<Aenumerator> res(segments + 1);
        int bound = 0;
        res[bound] = 1;
        Aenumerator count = 0;
        while (next()) if (++count % chunk_size == 0) res[++bound] = zeroes;
        res[bound] = stop;
        return res;
    }
};

// -----------------------------------------------------------------
// equiprobable A values merging
// -----------------------------------------------------------------
static char A_weight[1 << A_LENGTH];
struct Agroup {
    vectorB value;
    int     count;
    Agroup(Aenumerator a = 0, int length = 0) : value(a), count(length) {}
};
static vector<Agroup> A_groups;

Aenumerator reverse(Aenumerator n) // this works on N-1 bits for a N bits word
{
    Aenumerator res = 0;
    if (n != 0) // must have at least one bit set for the rest to work
    {
        // construct left-padded reverse value
        for (int i = 0; i != 8 * sizeof(n)-1; i++)
        {
            res |= (n & 1);
            res <<= 1;
            n >>= 1;
        }

        // shift right to elimitate trailing zeroes
        while (!(res & 1)) res >>= 1;
    }
    return res;
}

void generate_A_groups(int n)
{
    static bitset<1 << A_LENGTH> lookup(0);
    Aenumerator limit_A = (Aenumerator)pow(2, n);
    Aenumerator overflow = 1 << n;
    for (char & w : A_weight) w = 0;

    // gather rotation cycles
    for (Aenumerator a = 0; a != limit_A; a++)
    {
        Aenumerator rotated = a;
        int cycle_length = 0;
        for (int i = 0; i != n; i++)
        {
            // check for new cycles
            if (!lookup[rotated])
            {
                cycle_length++;
                lookup[rotated] = 1;
            }

            // rotate current value
            rotated <<= 1;
            if (rotated & overflow) rotated |= 1;
            rotated &= (overflow - 1);
        }

        // store new cycle
        if (cycle_length > 0) A_weight[a] = cycle_length;
    }

    // merge symetric groups
    for (Aenumerator a = 0; a != limit_A; a++)
    {
        // skip already grouped values
        if (A_weight[a] == 0) continue;

        // regroup a symetric pair
        Aenumerator r = reverse(a);
        if (r != a)
        {
            A_weight[a] += A_weight[r];
            A_weight[r] = 0;
        }  
    }

    // generate groups
    for (Aenumerator a = 0; a != limit_A; a++)
    {
        if (A_weight[a] != 0) A_groups.push_back(Agroup(a, A_weight[a]));
    }
}

// -----------------------------------------------------------------
// worker thread
// -----------------------------------------------------------------
OccurenceCounters solve(int n, int index, Aenumerator Bstart, Aenumerator Bstop)
{
    OccurenceCounters consecutive_zero_Z(n, 0);  // counts occurences of the first i terms of Z being 0

    // lock on assigned CPU
    lock_on_CPU(index);

    // enumerate B vectors
    Bgenerator Bgen(n, Bstart, Bstop);
    while (Bgen.next())
    {
        // get next B value
        Bvalue B = Bgen.value();

        // enumerate A vector groups
        for (const auto & group : A_groups)
        {
            // count consecutive occurences of inner product equal to zero
            vectorB sliding_A(group.value);
            for (int i = 0; i != n; i++)
            {
                if ((sliding_A & B.p1).count() != (sliding_A & B.m1).count()) break;
                consecutive_zero_Z[i] += group.count;
                sliding_A <<= 1;
            }
        }
    }
    return consecutive_zero_Z;
}

// -----------------------------------------------------------------
// main
// -----------------------------------------------------------------
#define die(msg) { cout << msg << endl; exit (-1); }

int main(int argc, char * argv[])
{
    int n = argc == 2 ? atoi(argv[1]) : 16; // arbitray value for debugging
    if (n < 1 || n > 24) die("vectors of lenght between 1 and 24 is all I can (try to) compute, guv");

    auto begin = time(NULL);

    // one worker thread per CPU
    int num_workers = number_of_CPUs();

    // regroup equiprobable A values
    generate_A_groups(n);

    // compute B generation ranges for proper load balancing
    vector<Aenumerator> ranges = Bgenerator(n).calibrate(num_workers);

    // set workers to work
    vector<future<OccurenceCounters>> workers(num_workers);
    for (int i = 0; i != num_workers; i++)
    {
        workers[i] = async(
            launch::async, // without this parameter, C++ will decide whether execution shall be sequential or asynchronous (isn't C++ fun?).
            solve, n, i, ranges[i], ranges[i+1]); 
    }

    // collect results
    OccurenceCounters result(n + 1, 0);
    for (auto& worker : workers)
    {
        OccurenceCounters partial = worker.get();
        for (size_t i = 0; i != partial.size(); i++) result[i] += partial[i]*2; // each result counts for a symetric B pair
    }
    for (Counter & res : result) res += (Counter)1 << n; // add null B vector contribution
    result[n] = result[n - 1];                           // the last two probabilities are equal by construction

    auto duration = time(NULL) - begin;

    // output
    cout << "done in " << duration / 60 << ":" << setw(2) << setfill('0') << duration % 60 << setfill(' ')
        << " by " << num_workers << " worker thread" << ((num_workers > 1) ? "s" : "") << endl;
    Counter events = (Counter)pow(6, n);
    int width = (int)log10(events) + 2;
    cout.precision(5);
    for (int i = 0; i <= n; i++) cout << setw(2) << i << setw(width) << result[i] << " / " << events << " " << fixed << (float)result[i] / events << endl;

    return 0;
}

Construyendo el ejecutable

Es una fuente independiente de C ++ 11 que se compila sin advertencias y se ejecuta sin problemas en:

• Win7 y MSVC2013
• Win7 y MinGW - g ++ 4.7
• Ubuntu & g ++ 4.8 (en una VM VirtualBox con 2 CPU asignadas)

Si compila con g ++, use: g ++ -O3 -pthread -std = c ++ 11
olvidarse -pthreadproducirá un volcado de núcleo agradable y amigable.

Optimizaciones

1. El último término Z es igual al primero (Bpre x A en ambos casos), por lo que los dos últimos resultados son siempre iguales, lo que dispensa el cálculo del último valor Z.
   La ganancia es despreciable, pero codificarla no cuesta nada, por lo que bien podría usarla.

2. Como descubrió Jakube, todos los valores cíclicos de un vector A dado producen las mismas probabilidades.
   Puede calcularlos con una sola instancia de A y multiplicar el resultado por el número de sus posibles rotaciones. Los grupos de rotación se pueden precalcular fácilmente en un tiempo descuidado, por lo que esta es una gran ganancia de velocidad neta.
   Dado que el número de permutaciones de un vector de longitud n es n-1, la complejidad cae de o (6 ⁿ ) a o (6 ⁿ / (n-1)), básicamente yendo un paso más allá para el mismo tiempo de cálculo.

3. Parece que los pares de patrones simétricos también generan las mismas probabilidades. Por ejemplo, 100101 y 101001.
   No tengo ninguna prueba matemática de eso, pero intuitivamente cuando se me presentan todos los patrones B posibles, cada valor A simétrico se enredará con el valor B simétrico correspondiente para el mismo resultado global.
   Esto permite reagrupar algunos vectores A más, para una reducción aproximada del 30% del número de grupos A.

4. Incorrecto Por alguna razón semi-misteriosa, todos los patrones con solo uno o dos bits establecidos producen el mismo resultado. Esto no representa tantos grupos distintos, pero aun así pueden fusionarse prácticamente sin costo.

5. Los vectores B y -B (B con todos los componentes multiplicados por -1) producen las mismas probabilidades.
   (por ejemplo [1, 0, -1, 1] y [-1, 0, 1, -1]).
   Excepto por el vector nulo (todos los componentes son iguales a 0), B y -B forman un par de vectores distintos .
   Esto permite reducir a la mitad el número de valores de B considerando solo uno de cada par y multiplicando su contribución por 2, agregando la contribución global conocida de B nulo a cada probabilidad solo una vez.

Cómo funciona

El número de valores B es enorme (3 ⁿ ), por lo que precalcularlos requeriría cantidades indecentes de memoria, lo que ralentizaría el cálculo y eventualmente agotaría la RAM disponible.
Desafortunadamente, no pude encontrar una manera simple de enumerar el medio conjunto de valores B optimizados, así que recurrí a la codificación de un generador dedicado.

El poderoso generador B fue muy divertido de codificar, aunque los lenguajes que admiten mecanismos de rendimiento habrían permitido programarlo de una manera mucho más elegante.
Lo que hace en pocas palabras es considerar el "esqueleto" de un vector Bpre como un vector binario donde los 1 representan valores reales de -1 o +1.
Entre todos estos valores potenciales + 1 / -1, el primero se fija en +1 (seleccionando así uno de los posibles vectores B / -B), y se enumeran todas las combinaciones posibles + 1 / -1 restantes.
Por último, un sistema de calibración simple garantiza que cada hilo de trabajo procesará un rango de valores de aproximadamente el mismo tamaño.

Los valores A están muy filtrados para reagruparse en trozos equiprobables.
Esto se hace en una fase previa a la computación que la fuerza bruta examina todos los valores posibles.
Esta parte tiene un tiempo de ejecución de O (2 ⁿ ) descuidado y no necesita ser optimizado (¡el código ya es lo suficientemente ilegible como es!).

Para evaluar el producto interno (que solo necesita ser probado contra cero), los componentes -1 y 1 de B se reagrupan en vectores binarios.
El producto interno es nulo si (y solo si) hay un número igual de + 1s y -1s entre los valores B correspondientes a valores A distintos de cero.
Esto se puede calcular con operaciones simples de enmascaramiento y conteo de bits, ayudado por std::bitseteso producirá un código de conteo de bits razonablemente eficiente sin tener que recurrir a instrucciones intrínsecas feas.

El trabajo se divide por igual entre los núcleos, con afinidad forzada de la CPU (cada poquito ayuda, o eso dicen).

Resultado de ejemplo

C:\Dev\PHP\_StackOverflow\C++\VectorCrunch>release\VectorCrunch.exe 16
done in 8:19 by 4 worker threads
 0  487610895942 / 2821109907456 0.17284
 1   97652126058 / 2821109907456 0.03461
 2   20659337010 / 2821109907456 0.00732
 3    4631534490 / 2821109907456 0.00164
 4    1099762394 / 2821109907456 0.00039
 5     302001914 / 2821109907456 0.00011
 6     115084858 / 2821109907456 0.00004
 7      70235786 / 2821109907456 0.00002
 8      59121706 / 2821109907456 0.00002
 9      56384426 / 2821109907456 0.00002
10      55686922 / 2821109907456 0.00002
11      55508202 / 2821109907456 0.00002
12      55461994 / 2821109907456 0.00002
13      55451146 / 2821109907456 0.00002
14      55449098 / 2821109907456 0.00002
15      55449002 / 2821109907456 0.00002
16      55449002 / 2821109907456 0.00002

Actuaciones

Multithreading debería funcionar perfectamente en esto, aunque solo los núcleos "verdaderos" contribuirán completamente a la velocidad de cálculo. Mi CPU tiene solo 2 núcleos para 4 CPU, y la ganancia sobre una versión de subproceso único es "solo" aproximadamente 3.5.

Compiladores

Un problema inicial con subprocesos múltiples me llevó a creer que los compiladores de GNU estaban funcionando peor que Microsoft.

Después de pruebas más exhaustivas, parece que g ++ gana el día una vez más, produciendo un código aproximadamente un 30% más rápido (la misma proporción que noté en otros dos proyectos pesados ​​en computación).

En particular, la std::bitsetbiblioteca se implementa con instrucciones de recuento de bits dedicadas por g ++ 4.8, mientras que MSVC 2013 usa solo bucles de cambios de bits convencionales.

Como era de esperar, compilar en 32 o 64 bits no hace ninguna diferencia.

Más refinamientos

Noté algunos grupos A que producen las mismas probabilidades después de todas las operaciones de reducción, pero no pude identificar un patrón que permitiera reagruparlos.

Aquí están los pares que noté para n = 11:

  10001011 and 10001101
 100101011 and 100110101
 100101111 and 100111101
 100110111 and 100111011
 101001011 and 101001101
 101011011 and 101101011
 101100111 and 110100111
1010110111 and 1010111011
1011011111 and 1011111011
1011101111 and 1011110111