23 added 17 characters in body

# woolly bully - C++ - n=16 in 8:19 on my i3-2100 (2 cores / 4 CPU)way too slow

A brute withWell since a couple of tricks up its sleeve better programmer took on the C++ implementation, I'm calling quits for this one.

# woolly bully - C++ - n=16 in 8:19 on my i3-2100 (2 cores / 4 CPU)

A brute with a couple of tricks up its sleeve

# woolly bully - C++ - way too slow

Well since a better programmer took on the C++ implementation, I'm calling quits for this one.

22 added 547 characters in body

# woolly bully - C++ - n=16 in 8:2719 on my i3-2100 (2 cores / 4 CPU)

EDIT3: Some more computation space reduction
EDIT4: number of B values cut in half, punching through the n=15 barrier

#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)
{
return res == 0 ? 8 : res;
}

#ifdef __linux__
#include <sched.h>
void lock_on_CPU(int cpu)
{
}
#elif defined (_WIN32)
#include <Windows.h>
#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
{
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++)
{
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;
}

// regroup 2 bits patterns with 1 bit patterns
if (vectorB(a).count() == 2)
{
A_weight[1] += A_weight[a];
A_weight[a] = 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]));
}
}

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

1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.

2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its possible rotations. Rotation groups can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6n) to o(6n/(n-1)), basically going a step further for the same computation time.

3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.

4. For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.WRONG For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.

5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance [ 1, 0,-1, 1] and [-1, 0, 1,-1]).
Except for the null vector (all components equal to 0), B and -B form a pair of distinct vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known global contribution of null B at the end of global computationto each probability only once.

C:\Dev\PHP\_StackOverflow\C++\VectorCrunch>release\VectorCrunch.exe 16
done in 8:2719 by 4 worker threads
0  483937576006487610895942 / 2821109907456 0.1715417284
1   9642768625097652126058 / 2821109907456 0.0341803461
2   2025119057820659337010 / 2821109907456 0.0071800732
3    44954858504631534490 / 2821109907456 0.0015900164
4    10544130181099762394 / 2821109907456 0.0003700039
5     286885626302001914 / 2821109907456 0.0001000011
6     110046266115084858 / 2821109907456 0.00004
7      6855642670235786 / 2821109907456 0.00002
8      5852709859121706 / 2821109907456 0.00002
9      5615140256384426 / 2821109907456 0.00002
10      5557442655686922 / 2821109907456 0.00002
11      5543588255508202 / 2821109907456 0.00002
12      5540220255461994 / 2821109907456 0.00002
13      5539553055451146 / 2821109907456 0.00002
14      5539468255449098 / 2821109907456 0.00002
15      5539485855449002 / 2821109907456 0.00002
16      5539485855449002 / 2821109907456 0.00002


# woolly bully - C++ - n=16 in 8:27 on my i3-2100 (2 cores / 4 CPU)

EDIT3: Some more computation space reduction
EDIT4: number of B values cut in half, punching through the n=15 barrier

#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)
{
return res == 0 ? 8 : res;
}

#ifdef __linux__
#include <sched.h>
void lock_on_CPU(int cpu)
{
}
#elif defined (_WIN32)
#include <Windows.h>
#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
{
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++)
{
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;
}

// regroup 2 bits patterns with 1 bit patterns
if (vectorB(a).count() == 2)
{
A_weight[1] += A_weight[a];
A_weight[a] = 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]));
}
}

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

1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.

2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its possible rotations. Rotation groups can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6n) to o(6n/(n-1)), basically going a step further for the same computation time.

3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.

4. For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.

5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance [ 1, 0,-1, 1] and [-1, 0, 1,-1]).
Except for the null vector (all components equal to 0), B and -B form a pair of distinct vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known contribution of null B at the end of global computation.

C:\Dev\PHP\_StackOverflow\C++\VectorCrunch>release\VectorCrunch.exe 16
done in 8:27 by 4 worker threads
0  483937576006 / 2821109907456 0.17154
1   96427686250 / 2821109907456 0.03418
2   20251190578 / 2821109907456 0.00718
3    4495485850 / 2821109907456 0.00159
4    1054413018 / 2821109907456 0.00037
5     286885626 / 2821109907456 0.00010
6     110046266 / 2821109907456 0.00004
7      68556426 / 2821109907456 0.00002
8      58527098 / 2821109907456 0.00002
9      56151402 / 2821109907456 0.00002
10      55574426 / 2821109907456 0.00002
11      55435882 / 2821109907456 0.00002
12      55402202 / 2821109907456 0.00002
13      55395530 / 2821109907456 0.00002
14      55394682 / 2821109907456 0.00002
15      55394858 / 2821109907456 0.00002
16      55394858 / 2821109907456 0.00002


# woolly bully - C++ - n=16 in 8:19 on my i3-2100 (2 cores / 4 CPU)

#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)
{
return res == 0 ? 8 : res;
}

#ifdef __linux__
#include <sched.h>
void lock_on_CPU(int cpu)
{
}
#elif defined (_WIN32)
#include <Windows.h>
#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
{
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++)
{
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]));
}
}

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

1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.

2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its possible rotations. Rotation groups can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6n) to o(6n/(n-1)), basically going a step further for the same computation time.

3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.

4. WRONG For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.

5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance [ 1, 0,-1, 1] and [-1, 0, 1,-1]).
Except for the null vector (all components equal to 0), B and -B form a pair of distinct vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known global contribution of null B to each probability only once.

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

21 added 5 characters in body
1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.

2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its permutationspossible rotations. Permutation setsRotation groups can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6n) to o(6n/(n-1)), basically going a step further for the same computation time.

3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.

4. For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.

5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance [ 1, 0,-1, 1] and [-1, 0, 1,-1]).
Except for the null vector (all components equal to 0), B and -B form a pair of distinct vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known contribution of null B at the end of global computation.

1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.

2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its permutations. Permutation sets can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6n) to o(6n/(n-1)), basically going a step further for the same computation time.

3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.

4. For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.

5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance [ 1, 0,-1, 1] and [-1, 0, 1,-1]).
Except for the null vector (all components equal to 0), B and -B form a pair of distinct vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known contribution of null B at the end of global computation.

1. The last Z term is equal to the first (Bpre x A in both cases), so the last two results are always equal, which dispenses of computing the last Z value.
The gain is neglectible, but coding it costs nothing so you might as well use it.

2. As Jakube discovered, all cyclic values of a given A vector produce the same probabilities.
You can compute these with a single instance of A and multiply the result by the number of its possible rotations. Rotation groups can easily be pre-computed in a neglectible amount of time, so this is a huge net speed gain.
Since the number of permutations of an n length vector is n-1, the complexity drops from o(6n) to o(6n/(n-1)), basically going a step further for the same computation time.

3. It appears pairs of symetric patterns also generate the same probabilities. For instance, 100101 and 101001.
I have no mathematical proof of that, but intuitively when presented with all possible B patterns, each symetric A value will be convoluted with the corresponding symetric B value for the same global outcome.
This allows to regroup some more A vectors, for an approximative 30% reduction of A groups number.

4. For some semi-mysterious reason, all patterns with only one or two bits set produce the same result. This does not represent that many distinct groups, but still they can be merged at virtually no cost.

5. The vectors B and -B (B with all components multiplied by -1) produce the same probabilities.
(for instance [ 1, 0,-1, 1] and [-1, 0, 1,-1]).
Except for the null vector (all components equal to 0), B and -B form a pair of distinct vectors.
This allows to cut the number of B values in half by considering only one of each pair and multiplying its contribution by 2, adding the known contribution of null B at the end of global computation.

20 added 3844 characters in body
19 deleted 22 characters in body
18 added 431 characters in body
17 added 502 characters in body
16 deleted 1 character in body
15 deleted 1 character in body
14 added 1412 characters in body
13 added 1412 characters in body
12 added 1412 characters in body
11 deleted 1 character in body
10 edited body
9 added 61 characters in body
8 added 1451 characters in body
7 deleted 10 characters in body
6 added 178 characters in body
5 added 220 characters in body
4 added 1143 characters in body
3 deleted 4 characters in body
2 added 69 characters in body
1