Diese Herausforderung ist teilweise eine Algorithmus-Herausforderung, beinhaltet etwas Mathematik und ist teilweise einfach eine schnellste Code-Herausforderung.

nBetrachten Sie für eine positive ganze Zahl eine gleichmäßig zufällige Zeichenfolge mit 1s und 0s Länge nund nennen Sie sie A. Betrachten sie nun auch einen zweiten einheitlich gewählt zufälligen Zeichenfolge der Länge , nderen Werte -1, 0,oder 1es nennen B_pre. Nun wollen wir Bsein B_pre+ B_pre. Das ist B_premit sich selbst verknüpft.

Betrachten Sie nun das innere Produkt von Aund B[j,...,j+n-1]und nennen Sie es Z_jund index from 1.

Aufgabe

Die Ausgabe sollte eine Liste von n+1Brüchen sein. Der idritte Term in der Ausgabe sollte die genaue Wahrscheinlichkeit sein, dass alle ersten iTerms Z_jmitj <= i gleich 0.

Ergebnis

Das größte n für die Ihr Code auf meinem Computer in weniger als 10 Minuten die richtige Ausgabe liefert.

Kabelbinder

Wenn zwei Antworten die gleiche Punktzahl haben, gewinnt die zuerst eingereichte.

In dem (sehr sehr) unwahrscheinlichen Fall, dass jemand eine Methode findet, um unbegrenzte Punktzahlen zu erhalten, wird der erste gültige Beweis für eine solche Lösung akzeptiert.

Hinweis

Versuchen Sie nicht, dieses Problem mathematisch zu lösen, es ist zu schwer. Meiner Ansicht nach ist es der beste Weg, zu den grundlegenden Definitionen der Wahrscheinlichkeit von der High School zurückzukehren und kluge Wege zu finden, um den Code dazu zu bringen, eine erschöpfende Aufzählung der Möglichkeiten vorzunehmen.

Sprachen und Bibliotheken

Sie können jede Sprache verwenden, die über einen frei verfügbaren Compiler / Interpreter / etc. Verfügt. für Linux und alle Bibliotheken, die auch für Linux frei verfügbar sind.

Meine Maschine Die Timings werden auf meinem Computer ausgeführt. Dies ist eine Ubuntu-Standardinstallation auf einem AMD FX-8350 Eight-Core-Prozessor. Dies bedeutet auch, dass ich in der Lage sein muss, Ihren Code auszuführen. Verwenden Sie daher nur leicht verfügbare kostenlose Software und fügen Sie vollständige Anweisungen zum Kompilieren und Ausführen Ihres Codes bei.

Einige Testausgaben. Betrachten Sie jeweils nur die erste Ausgabe n. Dann ist es soweit i=1. Für n1 bis 13 sollten sie sein.

 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

Eine allgemeine Formel für finden Sie auch i=1unter http://oeis.org/A081671 .

Leaderboard (aufgeschlüsselt nach Sprache)

• n = 15. Python + Parallelpython + Pypy in 1min49s von Jakube
• n = 17. C ++ in 3min37s von Keith Randall
• n = 16. C ++ in 2min38s von kuroi neko

C ++, n = 18 in 9 min auf 8 Threads

(Lassen Sie mich wissen, ob es auf Ihrem Computer weniger als 10 Minuten dauert.)

Ich nutze verschiedene Formen der Symmetrie im B-Array. Diese sind zyklisch (Verschiebung um eine Position), umgekehrt (Reihenfolge der Elemente umkehren) und vorzeichenbehaftet (Negativ jedes Elements nehmen). Zuerst berechne ich die Liste der Bs, die wir versuchen müssen, und ihr Gewicht. Dann durchläuft jedes B eine schnelle Routine (unter Verwendung von Bitcount-Befehlen) für alle 2 ^ n-Werte von A.

Hier ist das Ergebnis für 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

Kompilieren Sie das folgende Programm mit 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 mit pypy und pp: n = 15 in 3 Minuten

Auch nur eine einfache rohe Kraft. Interessant zu sehen, dass ich mit C ++ fast die gleiche Geschwindigkeit wie kuroi neko bekomme. Mein Code kann erreichenn = 12 in ca. 5 Minuten erreicht werden. Und ich führe es nur auf einem virtuellen Kern aus.

Bearbeiten: Reduziert den Suchraum um den Faktor n

Mir ist aufgefallen, dass ein zyklisierter Vektor A*von Adieselben Zahlen wie Wahrscheinlichkeiten (dieselben Zahlen) wie der ursprüngliche Vektor erzeugt, Awenn ich darüber iteriere B. Eg Der Vektor (1, 1, 0, 1, 0, 0)hat die gleichen Wahrscheinlichkeiten als jeder der Vektoren (1, 0, 1, 0, 0, 1), (0, 1, 0, 0, 1, 1), (1, 0, 0, 1, 1, 0), (0, 0, 1, 1, 0, 1)und , (0, 1, 1, 0, 1, 0)wenn eine zufällige Auswahl B. Deshalb muss ich nicht jede dieser sechs Vektoren durchlaufen, aber nur etwa 1 und ersetzen count[i] += 1mit count[i] += cycle_number.

Dies reduziert die Komplexität von Theta(n) = 6^nzu Theta(n) = 6^n / n. Deshalb ist n = 13es ungefähr 13 mal so schnell wie meine Vorgängerversion. Es rechnet n = 13in ca. 2 Minuten 20 Sekunden. Dafür ist n = 14es noch etwas zu langsam. Es dauert ungefähr 13 Minuten.

edit 2: Multi-Core-Programmierung

Mit der nächsten Verbesserung nicht wirklich zufrieden. Ich habe mich entschlossen, mein Programm auch auf mehreren Kernen auszuführen. Auf meinen 2 + 2 Kernen kann ich jetzt n = 14in ca. 7 Minuten rechnen . Nur ein Faktor 2 Verbesserung.

Der Code ist in diesem Github-Repo verfügbar: Link . Die Mehrkernprogrammierung macht das etwas hässlich.

Bearbeiten 3: Reduzieren des Suchraums für AVektoren und BVektoren

Ich bemerkte die gleiche Spiegelsymmetrie für die Vektoren Awie kuroi neko. Immer noch nicht sicher, warum das funktioniert (und ob es für jeden funktioniert)n ).

Die Reduzierung des Suchraums für BVektoren ist etwas cleverer. Ich habe die Erzeugung der Vektoren ( itertools.product) durch eine eigene Funktion ersetzt. Grundsätzlich beginne ich mit einer leeren Liste und lege sie auf einen Stapel. Bis der Stapel leer ist, entferne ich eine Liste. Wenn sie nicht dieselbe Länge hat wie n, erstelle ich 3 weitere Listen (indem ich -1, 0, 1 anhänge) und schiebe sie auf den Stapel. Wenn eine Liste die gleiche Länge hat wie n, kann ich die Summen auswerten.

Jetzt, da ich die Vektoren selbst generiere, kann ich sie filtern, je nachdem, ob ich die Summe = 0 erreichen kann oder nicht. Wenn z. B. mein Vektor Aist (1, 1, 1, 0, 0)und mein Vektor Baussieht (1, 1, ?, ?, ?), weiß ich, dass ich die nicht ?mit Werten füllen kann , so dass A*B = 0. Ich muss also nicht alle 6 Vektoren Bdes Formulars durchlaufen (1, 1, ?, ?, ?).

Wir können dies verbessern, wenn wir die Werte für 1 ignorieren. Wie in der Frage erwähnt, sind die Werte für i = 1die Sequenz A081671 . Es gibt viele Möglichkeiten, diese zu berechnen. Ich wähle die einfache Wiederholung: a(n) = (4*(2*n-1)*a(n-1) - 12*(n-1)*a(n-2)) / n. Da wir i = 1im Grunde genommen in kürzester Zeit rechnen können, können wir weitere Vektoren nach filtern B. ZB A = (0, 1, 0, 1, 1)und B = (1, -1, ?, ?, ?). Wir können Vektoren ignorieren, bei denen der erste ? = 1, weil der A * cycled(B) > 0, für alle diese Vektoren gilt. Ich hoffe du kannst folgen. Es ist wahrscheinlich nicht das beste Beispiel.

Damit kann ich n = 15in 6 Minuten rechnen .

bearbeiten 4:

Schnell umgesetzt kuroi Nekos große Idee, die besagt, dass Bund -Berzeugt die gleichen Ergebnisse. Beschleunigung x2. Die Implementierung ist jedoch nur ein kurzer Hack. n = 15in 3 minuten.

Code:

Den vollständigen Code finden Sie unter Github . Der folgende Code ist nur eine Darstellung der Hauptfunktionen. Ich habe Importe weggelassen, Multicore-Programmierung, Drucken der Ergebnisse, ...

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

Verwendung:

Du musst pypy installieren (für Python 2 !!!). Das parallele Python-Modul ist nicht für Python 3 portiert. Anschließend müssen Sie das parallele Python-Modul pp-1.6.4.zip installieren . Extrahieren Sie es cdin den Ordner und rufen Sie pypy setup.py install.

Dann kannst du mein Programm mit aufrufen

pypy you-do-the-math.py 15

Es wird automatisch die Anzahl der CPUs ermittelt. Nach dem Beenden des Programms kann es zu Fehlermeldungen kommen. Ignorieren Sie diese einfach. n = 16sollte auf Ihrem Rechner möglich sein.

Ausgabe:

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%

Notizen und Ideen:

• Ich habe einen i7-4600m Prozessor mit 2 Kernen und 4 Threads. Es spielt keine Rolle, ob ich 2 oder 4 Threads verwende. Die CPU-Auslastung beträgt 50% bei 2 Threads und 100% bei 4 Threads, dauert aber immer noch genauso lange. Ich weiß nicht warum. Ich habe überprüft, dass jeder Thread nur die Hälfte der Datenmenge hat, wenn es 4 Threads gibt, die Ergebnisse überprüft, ...
• Ich benutze viele Listen. Python ist nicht sehr effizient im Speichern, ich muss viele Listen kopieren, ... Also habe ich mir überlegt, stattdessen eine Ganzzahl zu verwenden. Ich könnte die Bits 00 (für 0) und 11 (für 1) im Vektor A und die Bits 10 (für -1), 00 (für 0) und 01 (für 1) im Vektor B verwenden von A und B müsste ich nur die A & BBlöcke 01 und 10 berechnen und zählen. Radfahren kann durch Verschieben des Vektors und Verwenden von Masken erfolgen. Ich habe das alles tatsächlich implementiert. Sie finden es in einigen meiner älteren Commits auf Github. Es stellte sich aber heraus, langsamer zu sein als mit Listen. Ich denke, pypy optimiert die Listenoperationen wirklich.

Wollmobber - C ++ - viel zu langsam

Nun, da ein besserer Programmierer die Implementierung von C ++ übernommen hat, rufe ich für dieses Programm das Beenden auf.

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

Erstellen der ausführbaren Datei

Es ist eine eigenständige C ++ 11-Quelle, die ohne Warnungen kompiliert und problemlos ausgeführt werden kann:

• Win7 & MSVC2013
• Win7 & MinGW - g ++ 4.7
• Ubuntu & g ++ 4.8 (in einer VirtualBox-VM mit 2 zugeordneten CPUs)

Wenn Sie mit g ++ kompilieren, geben Sie Folgendes ein : g ++ -O3 -pthread -std = c ++ 11 Wenn Sie etwas
vergessen, -pthreadwird ein netter und benutzerfreundlicher Core-Dump erstellt.

Optimierungen

1. Der letzte Z-Term ist gleich dem ersten (Bpre x A in beiden Fällen), so dass die letzten beiden Ergebnisse immer gleich sind, wodurch der letzte Z-Wert nicht berechnet werden muss.
   Der Gewinn ist vernachlässigbar, aber die Codierung kostet nichts, so dass Sie es genauso gut verwenden können.

2. Wie Jakube entdeckte, ergeben alle zyklischen Werte eines gegebenen A-Vektors die gleichen Wahrscheinlichkeiten.
   Sie können diese mit einer einzelnen Instanz von A berechnen und das Ergebnis mit der Anzahl der möglichen Umdrehungen multiplizieren. Rotationsgruppen können leicht in vernachlässigbarer Zeit vorberechnet werden, was einen enormen Netto-Geschwindigkeitsgewinn darstellt.
   Da die Anzahl der Permutationen eines n-Längenvektors n-1 beträgt, sinkt die Komplexität von o (6 ⁿ ) auf o (6 ⁿ / (n-1)), was im Grunde genommen für dieselbe Rechenzeit einen Schritt weiter geht.

3. Es scheint, dass Paare von symmetrischen Mustern auch die gleichen Wahrscheinlichkeiten erzeugen. Zum Beispiel 100101 und 101001.
   Ich habe keinen mathematischen Beweis dafür, aber intuitiv, wenn alle möglichen B-Muster dargestellt werden, wird jeder symetrische A-Wert mit dem entsprechenden symetrischen B-Wert für dasselbe globale Ergebnis verknüpft.
   Dies ermöglicht es, einige weitere A-Vektoren neu zu gruppieren, um die Anzahl der A-Gruppen um ca. 30% zu verringern.

4. FALSCH Aus einem halb mysteriösen Grund führen alle Muster mit nur einem oder zwei gesetzten Bits zum gleichen Ergebnis. Dies stellt nicht so viele unterschiedliche Gruppen dar, aber dennoch können sie praktisch kostenlos zusammengelegt werden.

5. Die Vektoren B und -B (B mit allen Komponenten multipliziert mit -1) ergeben die gleichen Wahrscheinlichkeiten.
   (zum Beispiel [1, 0, -1, 1] und [-1, 0, 1, -1]).
   Mit Ausnahme des Nullvektors (alle Komponenten sind gleich 0) bilden B und -B ein Paar verschiedener Vektoren.
   Dies ermöglicht es, die Anzahl der B-Werte zu halbieren, indem nur einer von jedem Paar betrachtet und sein Beitrag mit 2 multipliziert wird, wobei der bekannte globale Beitrag von Null B zu jeder Wahrscheinlichkeit nur einmal addiert wird.

Wie es funktioniert

Die Anzahl der B-Werte ist sehr groß (3 ⁿ ), sodass für die Vorberechnung unzureichender Speicher erforderlich ist, wodurch die Berechnung verlangsamt und der verfügbare RAM-Speicher möglicherweise erschöpft wird.
Leider konnte ich keinen einfachen Weg finden, um die Hälfte der optimierten B-Werte aufzulisten. Deshalb habe ich auf die Codierung eines dedizierten Generators zurückgegriffen.

Der mächtige B-Generator hat eine Menge Spaß beim Programmieren gemacht, obwohl Sprachen, die Ertragsmechanismen unterstützen, es erlaubt hätten, ihn viel eleganter zu programmieren.
Kurz gesagt: Betrachten Sie das "Skelett" eines Bpre-Vektors als einen binären Vektor, bei dem Einsen tatsächliche -1- oder +1-Werte darstellen.
Unter all diesen +1 / -1-Potentialwerten ist der erste fest auf +1 (wodurch einer der möglichen B / -B-Vektoren ausgewählt wird), und alle verbleibenden möglichen +1 / -1-Kombinationen werden aufgelistet.
Schließlich stellt ein einfaches Kalibrierungssystem sicher, dass jeder Worker-Thread einen Wertebereich von ungefähr derselben Größe verarbeitet.

A-Werte werden stark gefiltert, um sie in gleich wahrscheinlichen Blöcken neu zu gruppieren.
Dies geschieht in einer Vorberechnungsphase, in der Brute-Force alle möglichen Werte untersucht.
Dieser Teil hat eine vernachlässigbare O (2 ⁿ ) -Ausführungszeit und muss nicht optimiert werden (der Code ist bereits unlesbar genug!).

Um das innere Produkt (das nur gegen Null getestet werden muss) zu bewerten, werden die Komponenten -1 und 1 von B in binäre Vektoren umgruppiert.
Das innere Produkt ist null, wenn (und nur wenn) es eine gleiche Anzahl von +1 und -1 unter den B-Werten gibt, die Nicht-Null-A-Werten entsprechen.
Dies kann mit einfachen Maskierungs- und Bitzähloperationen berechnet std::bitsetwerden, wodurch ein einigermaßen effizienter Bitzählcode erzeugt wird, ohne auf hässliche intrinsische Anweisungen zurückgreifen zu müssen.

Die Arbeit ist zu gleichen Teilen auf die Kerne aufgeteilt, mit erzwungener CPU-Affinität (jedes kleine bisschen hilft, so heißt es).

Beispielergebnis

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

Vorstellungen

Multithreading sollte perfekt funktionieren, obwohl nur "echte" Kerne zur Rechengeschwindigkeit beitragen. Meine CPU hat nur 2 Kerne für 4 CPUs und der Gewinn gegenüber einer Singlethread-Version beträgt "nur" etwa 3,5.

Compiler

Ein anfängliches Problem mit Multithreading ließ mich glauben, dass die GNU-Compiler schlechter abschnitten als Microsoft.

Nach eingehenderen Tests scheint es, als ob g ++ den Tag erneut gewinnt und ungefähr 30% schnelleren Code produziert (das gleiche Verhältnis, das ich bei zwei anderen rechenintensiven Projekten festgestellt habe).

Insbesondere ist die std::bitsetBibliothek mit dedizierten Bitzählanweisungen von g ++ 4.8 implementiert, während MSVC 2013 nur Schleifen herkömmlicher Bitverschiebungen verwendet.

Wie zu erwarten war, macht das Kompilieren in 32 oder 64 Bit keinen Unterschied.

Weitere Verfeinerungen

Ich bemerkte ein paar A-Gruppen, die nach allen Reduktionsoperationen die gleichen Wahrscheinlichkeiten aufwiesen, konnte jedoch kein Muster identifizieren, das es mir ermöglichte, sie neu zu gruppieren.

Hier sind die Paare, die ich für n = 11 bemerkt habe:

  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