Ce défi est en partie un défi d'algorithmes, implique des mathématiques et est en partie simplement un défi de code le plus rapide.

Pour un entier positif n, considérez une chaîne uniformément aléatoire de 1s et 0de longueur net appelez-la A. Maintenant, considérons également une deuxième chaîne aléatoire de longueur uniformément choisie ndont les valeurs sont -1, 0,ou 1et appelez-la B_pre. Maintenant , nous allons Bêtre B_pre+ B_pre. Cela est B_preconcaténé à lui-même.

Considérez maintenant le produit intérieur de Aet B[j,...,j+n-1]et appelez-le Z_jet indexez à partir de 1.

Tâche

La sortie doit être une liste de n+1fractions. Le ie terme de la production devrait être la exacte probabilité que tous les premiers itermes Z_javec j <= iégale 0.

But

Le plus grand npour lequel votre code donne la sortie correcte en moins de 10 minutes sur ma machine.

Tie break

Si deux réponses ont le même score, celle soumise en premier l'emporte.

Dans le cas (très très) peu probable où quelqu'un trouve une méthode pour obtenir des scores illimités, la première preuve valable d'une telle solution sera acceptée.

Allusion

N'essayez pas de résoudre ce problème mathématiquement, c'est trop difficile. À mon avis, la meilleure façon est de revenir aux définitions de base des probabilités au secondaire et de trouver des façons intelligentes d'obtenir le code pour faire une énumération exhaustive des possibilités.

Langues et bibliothèques

Vous pouvez utiliser n'importe quelle langue disposant d'un compilateur / interprète / etc disponible gratuitement. pour Linux et toutes les bibliothèques qui sont également disponibles gratuitement pour Linux.

Ma machine Les horaires seront exécutés sur ma machine. Il s'agit d'une installation ubuntu standard sur un processeur AMD FX-8350 à huit cœurs. Cela signifie également que je dois pouvoir exécuter votre code. Par conséquent, n'utilisez que des logiciels gratuits facilement disponibles et veuillez inclure des instructions complètes sur la façon de compiler et d'exécuter votre code.

Quelques sorties de test. Considérez juste la première sortie pour chacun n. C'est à ce moment-là i=1. Pour n1 à 13, ils devraient l'être.

 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

Vous pouvez également trouver une formule générale pour i=1à http://oeis.org/A081671 .

Classement (divisé par langue)

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

C ++, n = 18 en 9 min sur 8 threads

(Faites-moi savoir si cela fait moins de 10 minutes sur votre machine.)

Je profite de plusieurs formes de symétrie dans le tableau B. Ceux-ci sont cycliques (décalage d'une position), inversion (inverser l'ordre des éléments) et signe (prendre le négatif de chaque élément). Je calcule d'abord la liste des Bs que nous devons essayer et leur poids. Ensuite, chaque B est exécuté à travers une routine rapide (en utilisant les instructions de bitcount) pour les 2 ^ n valeurs de A.

Voici le résultat pour 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

Compilez le programme ci-dessous avec 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 utilisant pypy et pp: n = 15 en 3 minutes

Aussi juste une simple force brute. Il est intéressant de voir que j'obtiens presque la même vitesse que kuroi neko avec C ++. Mon code peut atteindre n = 12environ 5 minutes. Et je ne l'exécute que sur un seul cœur virtuel.

modifier: réduire l'espace de recherche d'un facteur de n

J'ai remarqué qu'un vecteur cyclé A*de Aproduit les mêmes nombres que les probabilités (mêmes nombres) que le vecteur d'origine Alorsque je répète B. Par exemple, le vecteur(1, 1, 0, 1, 0, 0) a les mêmes probabilités que chacun des vecteurs (1, 0, 1, 0, 0, 1), (0, 1, 0, 0, 1, 1), (1, 0, 0, 1, 1, 0), (0, 0, 1, 1, 0, 1)et au (0, 1, 1, 0, 1, 0)moment de choisir un échantillon aléatoire B. Par conséquent, je n'ai pas à répéter sur chacun de ces 6 vecteurs, mais seulement environ 1 et à remplacer count[i] += 1par count[i] += cycle_number.

Cela réduit la complexité de Theta(n) = 6^nà Theta(n) = 6^n / n. Par conséquent, n = 13c'est environ 13 fois plus rapide que ma version précédente. Il calcule n = 13en environ 2 minutes 20 secondes. Car n = 14c'est encore un peu trop lent. Cela prend environ 13 minutes.

edit 2: Programmation multicœur

Pas vraiment satisfait de la prochaine amélioration. J'ai décidé d'essayer également d'exécuter mon programme sur plusieurs cœurs. Sur mes 2 + 2 cœurs, je peux maintenant calculer n = 14en environ 7 minutes. Seul un facteur d'amélioration de 2.

Le code est disponible dans ce dépôt github: Lien . La programmation multicœur rend un peu moche.

edit 3: Réduction de l'espace de recherche pour les Avecteurs et les Bvecteurs

J'ai remarqué la même symétrie miroir pour les vecteurs A que kuroi neko. Je ne sais toujours pas pourquoi cela fonctionne (et si cela fonctionne pour chacun n).

La réduction de l'espace de recherche de Bvecteurs est un peu plus intelligente. J'ai remplacé la génération des vecteurs ( itertools.product), par une fonction propre. Fondamentalement, je commence avec une liste vide et la mets sur une pile. Jusqu'à ce que la pile soit vide, je supprime une liste, si elle n'a pas la même longueur quen , je génère 3 autres listes (en ajoutant -1, 0, 1) et en les poussant sur la pile. Je une liste a la même longueur que n, je peux évaluer les sommes.

Maintenant que je génère moi-même les vecteurs, je peux les filtrer selon que je peux atteindre la somme = 0 ou non. Par exemple, si mon vecteur Aest (1, 1, 1, 0, 0), et mon vecteur Bsemble (1, 1, ?, ?, ?), je sais, que je ne peux pas remplir les ?valeurs, de sorte que A*B = 0. Je n'ai donc pas à parcourir tous les 6 vecteurs Bde la forme (1, 1, ?, ?, ?).

Nous pouvons améliorer cela si nous ignorons les valeurs de 1. Comme indiqué dans la question, car les valeurs de i = 1sont la séquence A081671 . Il existe de nombreuses façons de les calculer. Je choisis la récurrence simple: a(n) = (4*(2*n-1)*a(n-1) - 12*(n-1)*a(n-2)) / n. Comme nous pouvons calculer i = 1en un rien de temps, nous pouvons filtrer plus de vecteurs pour B. Par exemple A = (0, 1, 0, 1, 1)et B = (1, -1, ?, ?, ?). Nous pouvons ignorer les vecteurs, où le premier ? = 1, parce que A * cycled(B) > 0, pour tous ces vecteurs. J'espère que vous pourrez suivre. Ce n'est probablement pas le meilleur exemple.

Avec cela, je peux calculer n = 15en 6 minutes.

modifier 4:

La grande idée de kuroi neko mise en œuvre rapidement, qui dit cela Bet -Bproduit les mêmes résultats. Accélération x2. La mise en œuvre n'est qu'un rapide hack, cependant. n = 15en 3 minutes.

Code:

Pour le code complet, visitez Github . Le code suivant n'est qu'une représentation des principales fonctionnalités. J'ai laissé de côté les importations, la programmation multicœur, l'impression des résultats, ...

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

Usage:

Vous devez installer pypy (pour Python 2 !!!). Le module python parallèle n'est pas porté pour Python 3. Ensuite, vous devez installer le module python parallèle pp-1.6.4.zip . Extrayez-le cddans le dossier et appelez pypy setup.py install.

Ensuite, vous pouvez appeler mon programme avec

pypy you-do-the-math.py 15

Il déterminera automatiquement le nombre de processeurs. Il peut y avoir des messages d'erreur après avoir terminé le programme, ignorez-les. n = 16devrait être possible sur votre machine.

Production:

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%

Notes et idées:

• J'ai un processeur i7-4600m avec 2 cœurs et 4 threads. Peu importe si j'utilise 2 ou 4 fils. L'utilisation du processeur est de 50% avec 2 threads et 100% avec 4 threads, mais cela prend toujours le même temps. Je ne sais pas pourquoi. J'ai vérifié que chaque thread n'a que la moitié de la quantité de données, lorsqu'il y a 4 threads, j'ai vérifié les résultats, ...
• J'utilise beaucoup de listes. Python n'est pas assez efficace pour stocker, je dois copier beaucoup de listes, ... J'ai donc pensé à utiliser un entier à la place. Je pourrais utiliser les bits 00 (pour 0) et 11 (pour 1) dans le vecteur A, et les bits 10 (pour -1), 00 (pour 0) et 01 (pour 1) dans le vecteur B. Pour le produit de A et B, je n'aurais qu'à calculer A & Bet compter les blocs 01 et 10. Le cyclisme peut être fait en déplaçant le vecteur et en utilisant des masques, ... J'ai en fait implémenté tout cela, vous pouvez le trouver dans certaines de mes anciennes validations sur Github. Mais il s'est avéré être plus lent qu'avec les listes. Je suppose que pypy optimise vraiment les opérations de liste.

intimidateur laineux - C ++ - beaucoup trop lent

Eh bien, depuis qu'un meilleur programmeur a pris en charge l'implémentation C ++, j'appelle quits pour celle-ci.

#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;
}

Construire l'exécutable

Il s'agit d'une source C ++ 11 autonome qui se compile sans avertissements et s'exécute sans problème sur:

• Win7 et MSVC2013
• Win7 et MinGW - g ++ 4.7
• Ubuntu & g ++ 4.8 (dans une VM VirtualBox avec 2 CPU alloués)

Si vous compilez avec g ++, utilisez: g ++ -O3 -pthread -std = c ++ 11
oublier -pthreadproduira un vidage de mémoire agréable et convivial.

Optimisations

1. Le dernier terme Z est égal au premier (Bpre x A dans les deux cas), donc les deux derniers résultats sont toujours égaux, ce qui dispense de calculer la dernière valeur Z.
   Le gain est négligeable, mais le coder ne coûte rien donc vous pouvez aussi bien l'utiliser.

2. Comme Jakube l'a découvert, toutes les valeurs cycliques d'un vecteur A donné produisent les mêmes probabilités.
   Vous pouvez les calculer avec une seule instance de A et multiplier le résultat par le nombre de ses rotations possibles. Les groupes de rotation peuvent facilement être précalculés en un temps négligeable, il s'agit donc d'un gain de vitesse net énorme.
   Étant donné que le nombre de permutations d'un vecteur de longueur n est n-1, la complexité passe de o (6 ⁿ ) à o (6 ⁿ / (n-1)), allant fondamentalement plus loin pour le même temps de calcul.

3. Il semble que des paires de motifs symétriques génèrent également les mêmes probabilités. Par exemple, 100101 et 101001.
   Je n'ai aucune preuve mathématique de cela, mais intuitivement lorsqu'il est présenté avec tous les modèles B possibles, chaque valeur A symétrique sera alambiquée avec la valeur B symétrique correspondante pour le même résultat global.
   Cela permet de regrouper quelques vecteurs A supplémentaires, pour une réduction approximative de 30% du nombre de groupes A.

4. FAUX Pour une raison semi-mystérieuse, tous les motifs avec seulement un ou deux bits définis produisent le même résultat. Cela ne représente pas autant de groupes distincts, mais ils peuvent néanmoins être fusionnés pratiquement sans frais.

5. Les vecteurs B et -B (B avec toutes les composantes multipliées par -1) produisent les mêmes probabilités.
   (par exemple [1, 0, -1, 1] et [-1, 0, 1, -1]).
   À l'exception du vecteur nul (toutes les composantes égales à 0), B et -B forment une paire de vecteurs distincts .
   Cela permet de réduire de moitié le nombre de valeurs B en ne considérant qu'une seule de chaque paire et en multipliant sa contribution par 2, en ajoutant une seule fois la contribution globale connue de B nul à chaque probabilité.

Comment ça fonctionne

Le nombre de valeurs B est énorme (3 ⁿ ), donc leur pré-calcul nécessiterait des quantités de mémoire indécentes, ce qui ralentirait le calcul et finirait par épuiser la RAM disponible.
Malheureusement, je n'ai pas pu trouver un moyen simple d'énumérer le demi-ensemble de valeurs B optimisées, j'ai donc eu recours au codage d'un générateur dédié.

Le puissant générateur B était très amusant à coder, bien que les langages qui prennent en charge les mécanismes de rendement auraient permis de le programmer de manière beaucoup plus élégante.
Ce qu'il fait en un mot, c'est considérer le "squelette" d'un vecteur Bpre comme un vecteur binaire où les 1 représentent les valeurs réelles de -1 ou +1.
Parmi toutes ces valeurs potentielles + 1 / -1, la première est fixée à +1 (sélectionnant ainsi l'un des vecteurs B / -B possibles), et toutes les combinaisons possibles + 1 / -1 restantes sont énumérées.
Enfin, un système d'étalonnage simple garantit que chaque thread de travail traitera une plage de valeurs d'environ la même taille.

Les valeurs A sont fortement filtrées pour le regroupement en morceaux équiprobables.
Cela se fait dans une phase de pré-calcul où la force brute examine toutes les valeurs possibles.
Cette partie a un temps d'exécution O (2 ⁿ ) négligeable et n'a pas besoin d'être optimisée (le code étant déjà suffisamment lisible tel quel!).

Pour évaluer le produit interne (qui n'a besoin d'être testé que contre zéro), les composantes -1 et 1 de B sont regroupées en vecteurs binaires.
Le produit interne est nul si (et seulement si) il y a un nombre égal de + 1 et -1 parmi les valeurs B correspondant à des valeurs A non nulles.
Cela peut être calculé avec des opérations de masquage et de comptage de bits simples, aidées par std::bitsetcela produira un code de comptage de bits raisonnablement efficace sans avoir à recourir à des instructions intrinsèques laides.

Le travail est également réparti entre les cœurs, avec une affinité CPU forcée (chaque petit geste aide, ou du moins disent-ils).

Exemple de résultat

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

Les performances

Le multithreading devrait fonctionner parfaitement à ce sujet, bien que seuls les "vrais" cœurs contribuent pleinement à la vitesse de calcul. Mon processeur n'a que 2 cœurs pour 4 processeurs, et le gain par rapport à une version à thread unique est "seulement" d'environ 3,5.

Compilateurs

Un premier problème avec le multithreading m'a amené à croire que les compilateurs GNU fonctionnaient moins bien que Microsoft.

Après des tests plus approfondis, il semble que g ++ gagne encore une fois la journée, produisant un code environ 30% plus rapide (le même ratio que j'ai remarqué sur deux autres projets lourds en calcul).

Notamment, la std::bitsetbibliothèque est implémentée avec des instructions de comptage de bits dédiées par g ++ 4.8, tandis que MSVC 2013 utilise uniquement des boucles de décalages de bits conventionnels.

Comme on aurait pu s'y attendre, la compilation en 32 ou 64 bits ne fait aucune différence.

Améliorations supplémentaires

J'ai remarqué quelques groupes A produisant les mêmes probabilités après toutes les opérations de réduction, mais je n'ai pas pu identifier un schéma qui permettrait de les regrouper.

Voici les paires que j'ai remarquées pour 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