To wyzwanie jest częściowo wyzwaniem algorytmicznym, wymaga pewnej matematyki, a częściowo jest po prostu najszybszym wyzwaniem w kodzie.

Dla pewnej dodatniej liczby całkowitej nrozważ jednolicie losowy ciąg 1si i 0s długości ni wywołaj go A. Teraz rozważ także drugi losowo wybrany losowo ciąg długości, nktórego wartości są -1, 0,lub 1wywołaj go B_pre. Teraz niech Bbędzie B_pre+ B_pre. To jest B_prepołączone z samym sobą.

Rozważmy teraz wewnętrzną produkt Ai B[j,...,j+n-1]i nazwać Z_ji indeks z 1.

Zadanie

Dane wyjściowe powinny być listą n+1ułamków. iTh termin na wyjściu powinien być dokładny prawdopodobieństwo, że wszystko z pierwszych iterminów Z_jz j <= irównym 0.

Wynik

Największy, ndla którego Twój kod zapewnia prawidłowe wyjście w mniej niż 10 minut na moim komputerze.

Tie Breaker

Jeśli dwie odpowiedzi mają ten sam wynik, wygrywa ten, który przesłał pierwszy.

W (bardzo, bardzo) mało prawdopodobnym przypadku, gdy ktoś znajdzie sposób na uzyskanie nieograniczonej liczby punktów, zostanie zaakceptowany pierwszy ważny dowód takiego rozwiązania.

Wskazówka

Nie próbuj rozwiązywać tego problemu matematycznie, jest to zbyt trudne. Moim zdaniem najlepszym sposobem jest powrót do podstawowych definicji prawdopodobieństwa z liceum i znalezienie sprytnych sposobów, aby kod wykonał wyczerpujące wyliczenie możliwości.

Języki i biblioteki

Możesz używać dowolnego języka, który ma swobodnie dostępny kompilator / tłumacz / etc. dla systemu Linux i dowolnych bibliotek, które są również bezpłatnie dostępne dla systemu Linux.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350. Oznacza to również, że muszę być w stanie uruchomić Twój kod. W związku z tym korzystaj tylko z łatwo dostępnego bezpłatnego oprogramowania i dołącz pełne instrukcje, jak skompilować i uruchomić kod.

Niektóre wyjścia testowe. Rozważ tylko pierwsze wyjście dla każdego n. Wtedy to i=1. Dla nod 1 do 13 powinny być.

 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

Możesz również znaleźć ogólną formułę i=1na http://oeis.org/A081671 .

Tabela liderów (w podziale na język)

• n = 15. Python + równoległy python + pypy w 1min49s przez Jakube
• n = 17. C ++ w 3min37s autorstwa Keitha Randalla
• n = 16. C ++ w 2min38s autorstwa kuroi neko

C ++, n = 18 w 9 minut na 8 wątków

(Daj mi znać, jeśli na twoim komputerze będzie to mniej niż 10 minut).

Korzystam z kilku form symetrii w tablicy B. Są one cykliczne (przesunięcie o jedną pozycję), odwrócenie (odwróć kolejność elementów) i znak (weź minus każdego elementu). Najpierw obliczam listę Bs, którą musimy wypróbować i ich wagę. Następnie każdy B jest uruchamiany przez szybką procedurę (przy użyciu instrukcji bitcount) dla wszystkich 2 ^ n wartości A.

Oto wynik dla 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

Skompiluj poniższy program za pomocą 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 używa pypy i pp: n = 15 w 3 minuty

Również zwykła brutalna siła. Ciekawe, że mam prawie taką samą prędkość jak kuroi neko z C ++. Mój kod może dotrzeć n = 12w ciągu około 5 minut. I uruchamiam go tylko na jednym wirtualnym rdzeniu.

edycja: Zmniejsz przestrzeń wyszukiwania o współczynnik n

Zauważyłem, że cyklicznie wektor A*z Aprodukuje te same numery jak prawdopodobieństw (te same numery) jako oryginalnego wektora Akiedy iteracyjnego B. Np Wektor (1, 1, 0, 1, 0, 0)ma takie same prawdopodobieństwa jak każdy z wektorów (1, 0, 1, 0, 0, 1), (0, 1, 0, 0, 1, 1), (1, 0, 0, 1, 1, 0), (0, 0, 1, 1, 0, 1)i (0, 1, 1, 0, 1, 0)przy wyborze losowym B. Dlatego nie muszę iterować każdego z tych 6 wektorów, ale tylko około 1 i zastąpić count[i] += 1je count[i] += cycle_number.

Zmniejsza to złożoność od Theta(n) = 6^ndo Theta(n) = 6^n / n. Dlatego n = 13jest to około 13 razy szybsze niż moja poprzednia wersja. Oblicza się n = 13w ciągu około 2 minut i 20 sekund. Bo n = 14wciąż jest trochę za wolno. Zajmuje to około 13 minut.

edycja 2: Programowanie wielordzeniowe

Niezbyt zadowolony z następnej poprawy. Postanowiłem także spróbować uruchomić mój program na wielu rdzeniach. Na moich rdzeniach 2 + 2 mogę teraz obliczyć n = 14w około 7 minut. Tylko współczynnik poprawy 2.

Kod jest dostępny w tym repozytorium github: Link . Wielordzeniowe programowanie sprawia, że ​​jest trochę brzydkie.

edycja 3: Zmniejszenie przestrzeni wyszukiwania Awektorów i Bwektorów

Zauważyłem tę samą lustrzaną symetrię wektorów, Aco kuroi neko. Nadal nie jestem pewien, dlaczego to działa (i czy działa dla każdego n).

Zmniejszenie przestrzeni wyszukiwania Bwektorów jest nieco sprytniejsze. Generację wektorów ( itertools.product) zastąpiłem własną funkcją. Zasadniczo zaczynam od pustej listy i umieszczam ją na stosie. Dopóki stos nie będzie pusty, usuwam listę, jeśli nie ma takiej samej długości jak n, generuję 3 inne listy (dodając -1, 0, 1) i wypychając je na stos. Jeśli lista ma taką samą długość n, mogę ocenić sumy.

Teraz, gdy sam generuję wektory, mogę je filtrować w zależności od tego, czy mogę osiągnąć sumę = 0, czy nie. Np. Jeśli mój wektor Ajest (1, 1, 1, 0, 0)i mój wektor Bwygląda (1, 1, ?, ?, ?), wiem, że nie mogę wypełnić ?wartości, więc to A*B = 0. Więc nie muszę powtarzać tych wszystkich 6 wektorów Bformularza (1, 1, ?, ?, ?).

Możemy to poprawić, jeśli zignorujemy wartości dla 1. Jak zauważono w pytaniu, wartościami i = 1są sekwencja A081671 . Istnieje wiele sposobów ich obliczania. Wybieram proste nawrotom: a(n) = (4*(2*n-1)*a(n-1) - 12*(n-1)*a(n-2)) / n. Ponieważ możemy obliczyć i = 1w zasadzie w krótkim czasie, możemy filtrować więcej wektorów B. Np A = (0, 1, 0, 1, 1)a B = (1, -1, ?, ?, ?). Możemy zignorować wektory, gdzie pierwszy ? = 1, ponieważ A * cycled(B) > 0, dla wszystkich tych wektorów. Mam nadzieję, że możesz śledzić. To chyba nie najlepszy przykład.

Dzięki temu mogę obliczyć n = 15w 6 minut.

edycja 4:

Szybko wprowadzono świetny pomysł kuroi neko, który mówi, że Bi -Bdaje te same wyniki. Przyspieszenie x2. Wdrożenie to tylko szybki hack. n = 15za 3 minuty.

Kod:

Aby uzyskać pełny kod, odwiedź Github . Poniższy kod stanowi jedynie przedstawienie głównych funkcji. Pominąłem importowanie, programowanie wielordzeniowe, drukowanie wyników, ...

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

Stosowanie:

Musisz zainstalować pypy (dla Python 2 !!!). Równoległy moduł python nie jest portowany dla Pythona 3. Następnie musisz zainstalować równoległy moduł python pp-1.6.4.zip . Wyodrębnij go cddo folderu i zadzwoń pypy setup.py install.

Następnie możesz wywołać mój program za pomocą

pypy you-do-the-math.py 15

Automatycznie określi liczbę procesorów. Po zakończeniu programu mogą pojawić się komunikaty o błędach, po prostu je zignoruj. n = 16powinno być możliwe na twoim komputerze.

Wynik:

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%

Uwagi i pomysły:

• Mam procesor i7-4600m z 2 rdzeniami i 4 wątkami. Nie ma znaczenia, jeśli użyję 2 lub 4 wątków. Zużycie procesora wynosi 50% przy 2 wątkach i 100% przy 4 wątkach, ale nadal zajmuje to tyle samo czasu. Nie wiem dlaczego. Sprawdziłem, że każdy wątek ma tylko połowę ilości danych, gdy są 4 wątki, sprawdziłem wyniki ...
• Używam wielu list. Python nie jest dość wydajny w przechowywaniu, muszę skopiować wiele list, ... Więc pomyślałem o użyciu liczby całkowitej. Mógłbym użyć bitów 00 (dla 0) i 11 (dla 1) w wektorze A, a bitów 10 (dla -1), 00 (dla 0) i 01 (dla 1) w wektorze B. Dla produktu z A i B, musiałbym tylko obliczyć A & Bi policzyć 01 i 10 bloków. Cyklowanie można wykonać, przesuwając wektor i używając masek ... Właściwie to wszystko zaimplementowałem, można to znaleźć w niektórych moich starszych wersjach na Githubie. Okazało się jednak, że jest wolniejszy niż w przypadku list. Chyba pypy naprawdę optymalizuje operacje na listach.

woolly bully - C ++ - zdecydowanie za wolny

Cóż, ponieważ lepszy programista przyjął implementację C ++, wzywam do tego.

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

Budowanie pliku wykonywalnego

Jest to samodzielne źródło C ++ 11, które kompiluje się bez ostrzeżeń i działa płynnie na:

• Win7 i MSVC2013
• Win7 i MinGW - g ++ 4.7
• Ubuntu i g ++ 4.8 (w maszynie wirtualnej VirtualBox z przydzielonymi 2 procesorami)

Jeśli kompilujesz za pomocą g ++, użyj: g ++ -O3 -pthread -std = c ++ 11
zapomnienie o -pthreadutworzy miły i przyjazny zrzut pamięci .

Optymalizacje

1. Ostatni człon Z jest równy pierwszemu (Bpre x A w obu przypadkach), więc dwa ostatnie wyniki są zawsze równe, co zrezygnuje z obliczenia ostatniej wartości Z.
   Zysk jest pomijalny, ale kodowanie nic nie kosztuje, więc równie dobrze możesz go użyć.

2. Jak odkrył Jakube, wszystkie wartości cykliczne danego wektora A dają takie same prawdopodobieństwa.
   Można je obliczyć za pomocą pojedynczego wystąpienia A i pomnożyć wynik przez liczbę możliwych obrotów. Grupy rotacji można łatwo wstępnie obliczyć w nieistotnym czasie, więc jest to ogromny wzrost prędkości netto.
   Ponieważ liczba permutacji wektora długości n wynosi n-1, złożoność spada z o (6 ⁿ ) do o (6 ⁿ / (n-1)), zasadniczo idąc krok dalej w tym samym czasie obliczeń.

3. Wydaje się, że pary wzorów symetrycznych również generują takie same prawdopodobieństwa. Na przykład 100101 i 101001.
   Nie mam matematycznego dowodu na to, ale intuicyjnie, gdy zostaną przedstawione wszystkie możliwe wzorce B, każda symetryczna wartość A zostanie spleciona z odpowiednią symetryczną wartością B dla tego samego wyniku globalnego.
   Pozwala to na zgrupowanie kilku wektorów A, w celu około 30% zmniejszenia liczby grup A.

4. ŹLE Z jakiegoś na wpół tajemniczego powodu wszystkie wzory zawierające tylko jeden lub dwa bity dają ten sam rezultat. Nie reprezentuje to tak wielu odrębnych grup, ale nadal można je łączyć praktycznie bez żadnych kosztów.

5. Wektory B i -B (B ze wszystkimi składnikami pomnożonymi przez -1) dają takie same prawdopodobieństwa.
   (na przykład [1, 0, -1, 1] i [-1, 0, 1, -1]).
   Z wyjątkiem wektora zerowego (wszystkie składniki równe 0), B i -B tworzą parę różnych wektorów.
   Pozwala to zmniejszyć liczbę wartości B o połowę, biorąc pod uwagę tylko jedną z każdej pary i mnożąc jej udział przez 2, dodając znany globalny udział zerowy B do każdego prawdopodobieństwa tylko raz.

Jak to działa

Liczba wartości B jest ogromna (3 ⁿ ), więc ich wstępne obliczenie wymagałoby nieprzyzwoitej ilości pamięci, co spowolniłoby obliczenia i ostatecznie wyczerpało dostępną pamięć RAM.
Niestety nie mogłem znaleźć prostego sposobu wyliczenia połowy zestawu zoptymalizowanych wartości B, więc postanowiłem zakodować dedykowany generator.

Potężny generator B był świetną zabawą w kodowaniu, chociaż języki obsługujące mechanizmy wydajności pozwoliłyby na programowanie go w znacznie bardziej elegancki sposób.
To, co robi w skrócie, to „szkielet” wektora Bpre jako wektora binarnego, gdzie 1s reprezentują rzeczywiste wartości -1 lub +1.
Spośród wszystkich tych wartości potencjału + 1 / -1 pierwsza jest ustalona na +1 (w ten sposób wybierając jeden z możliwych wektorów B / -B), a wszystkie pozostałe możliwe kombinacje + 1 / -1 są wyliczone.
Wreszcie, prosty system kalibracji zapewnia, że ​​każdy wątek roboczy przetworzy zakres wartości w przybliżeniu tego samego rozmiaru.

Wartości są silnie filtrowane w celu przegrupowania w odpowiednie części.
Odbywa się to w fazie obliczeń wstępnych, w których brutalna siła sprawdza wszystkie możliwe wartości.
Ta część ma pomijalny czas wykonania O (2 ⁿ ) i nie musi być optymalizowana (kod jest już wystarczająco nieczytelny, jak jest!).

Aby ocenić iloczyn wewnętrzny (który należy przetestować tylko w odniesieniu do zera), komponenty -1 i 1 B są zgrupowane w wektorach binarnych.
Iloczyn wewnętrzny jest zerowy, jeśli (i tylko wtedy) istnieje równa liczba + 1s i -1s wśród wartości B odpowiadających niezerowym wartościom A.
Można to obliczyć za pomocą prostych operacji maskowania i liczenia bitów, std::bitsetco pomoże wygenerować dość wydajny kod liczby bitów bez konieczności uciekania się do brzydkich instrukcji wewnętrznych.

Praca jest równo podzielona między rdzenie, z wymuszonym powinowactwem do procesora (co nieco pomaga, a przynajmniej tak mówią).

Przykładowy wynik

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

Występy

Wielowątkowość powinna na tym doskonale działać, chociaż tylko „prawdziwe” rdzenie w pełni przyczynią się do szybkości obliczeń. Mój procesor ma tylko 2 rdzenie dla 4 procesorów, a zysk w stosunku do wersji jednowątkowej wynosi „tylko” około 3,5.

Kompilatory

Początkowy problem z wielowątkowością doprowadził mnie do wniosku, że kompilatory GNU działały gorzej niż Microsoft.

Po dokładniejszych testach okazuje się, że g ++ po raz kolejny wygrywa, generując około 30% szybszy kod (ten sam współczynnik, który zauważyłem w dwóch innych projektach wymagających dużej mocy obliczeniowej).

W szczególności std::bitsetbiblioteka jest implementowana z dedykowanymi instrukcjami liczenia bitów przez g ++ 4.8, podczas gdy MSVC 2013 wykorzystuje tylko pętle konwencjonalnych przesunięć bitów.

Jak można się było spodziewać, kompilacja w 32 lub 64 bitach nie ma znaczenia.

Dalsze udoskonalenia

Zauważyłem kilka grup A wytwarzających takie same prawdopodobieństwa po wszystkich operacjach redukcji, ale nie mogłem zidentyfikować wzoru, który pozwoliłby je zgrupować.

Oto pary, które zauważyłem dla 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