Computer: you do the math

Question

This challenge is partly an algorithms challenge, involves some math and is partly simply a fastest code challenge.

For some positive integer n, consider a uniformly random string of 1s and 0s of length n and call it A. Now also consider a second uniformly chosen random string of length n whose values are -1, 0, or 1 and call it B_pre. Now let B be B_pre+B_pre. That is B_pre concatenated to itself.

Now consider the inner product of A and B[j,...,j+n-1] and call it Z_j and index from 1.

Task

The output should be a list of n+1 fractions. The ith term in the output should be the exact probability that all of the first i terms Z_j with j <= i equal 0.

Score

The largest n for which your code gives the correct output in under 10 minutes on my machine.

Tie Breaker

If two answers have the same score, the one submitted first wins.

In the (very very) unlikely event that someone finds a method to get unlimited scores, the first valid proof of such a solution will be accepted.

Hint

Do not try to solve this problem mathematically, it's too hard. The best way in my view is to go back to basic definitions of probability from high school and find smart ways to get the code to do an exhaustive enumeration of the possibilities.

Languages and libraries

You can use any language which has a freely available compiler/interpreter/etc. for Linux and any libraries which are also freely available for Linux.

My machine The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code. As a consequence, only use easily available free software and please include full instructions how to compile and run your code.

---

Some test outputs. Consider just the first output for each n. That is when i=1. For n from 1 to 13 they should be.

 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

You can also find a general formula for i=1 at http://oeis.org/A081671 .

Leaderboard (split by language)

• n = 15. Python + parallel python + pypy in 1min49s by Jakube
• n = 17. C++ in 3min37s by Keith Randall
• n = 16. C++ in 2min38s by kuroi neko

Answer 1

C++, n=18 in 9 min on 8 threads

(Let me know if it makes it under 10 minutes on your machine.)

I take advantage of several forms of symmetry in the B array. Those are cyclic (shift by one position), reversal (reverse the order of the elements), and sign (take the negative of each element). First I compute the list of Bs we need to try and their weight. Then each B is run through a fast routine (using bitcount instructions) for all 2^n values of A.

Here's the result for n==18:

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

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

Compile the program below with 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;
}

Answer 2

Python 2 using pypy and pp: n = 15 in 3 minutes

Also just a simple brute force. Interesting to see, that I nearly get the same same speed as kuroi neko with C++. My code can reach n = 12 in about 5 minutes. And I only run it on one virtual core.

edit: Reduce search space by a factor of n

I noticed, that a cycled vector A* of A produces the same numbers as probabilities (same numbers) as the original vector A when I iterate over B. E.g. The vector (1, 1, 0, 1, 0, 0) has the same probabilities as each of the vectors (1, 0, 1, 0, 0, 1), (0, 1, 0, 0, 1, 1), (1, 0, 0, 1, 1, 0), (0, 0, 1, 1, 0, 1) and (0, 1, 1, 0, 1, 0) when choosing a random B. Therefore I don't have to iterate over each of these 6 vectors, but only about 1 and replace count[i] += 1 with count[i] += cycle_number.

This reduces the complexity from Theta(n) = 6^n to Theta(n) = 6^n / n. Therefore for n = 13 it's about 13 times as fast as my previous version. It calculates n = 13 in about 2 minutes 20 seconds. For n = 14 it is still a little bit too slow. It takes about 13 minutes.

edit 2: Multi-core-programming

Not really happy with the next improvement. I decided to also try to execute my program on multiple cores. On my 2 + 2 cores I now can calculate n = 14 in about 7 minutes. Only a factor of 2 improvement.

The code is available in this github repo: Link. The multi-core programming makes is a little bit ugly.

edit 3: Reducing search space for A vectors and B vectors

I noticed the same mirror-symmetry for the vectors A as kuroi neko did. Still not sure, why this works (and if it works for each n).

The reducing of the search space for B vectors is a bit cleverer. I replaced the generation of the vectors (itertools.product), with an own function. Basically, I start with an empty list and put it on a stack. Until the stack is empty, I remove a list, if it has not the same lenght as n, I generate 3 other lists (by appending -1, 0, 1) and pushing them onto the stack. I a list has the same length as n, I can evaluate the sums.

Now that I generate the vectors myself, I can filter them depending on if I can reach the sum = 0 or not. E.g. if my vector A is (1, 1, 1, 0, 0), and my vector B looks (1, 1, ?, ?, ?), I know, that I cannot fill the ? with values, so that A*B = 0. So I don't have to iterate over all those 6 vectors B of the form (1, 1, ?, ?, ?).

We can improve on this, if we ignore the values for 1. As noted in the question, for the values for i = 1 are the sequence A081671. There are many ways to calculate those. I choose the simple recurrence: a(n) = (4*(2*n-1)*a(n-1) - 12*(n-1)*a(n-2)) / n. Since we can calculate i = 1 in basically no time, we can filter more vectors for B. E.g. A = (0, 1, 0, 1, 1) and B = (1, -1, ?, ?, ?). We can ignore vectors, where the first ? = 1, because the A * cycled(B) > 0, for all these vectors. I hope you can follow. It's probably not the best example.

With this I can calculate n = 15 in 6 minutes.

edit 4:

Quickly implemented kuroi neko's great idea, which says, that B and -B produces the same results. Speedup x2. Implementation is only a quick hack, though. n = 15 in 3 minutes.

Code:

For the complete code visit Github. The following code is only a representation of the main features. I left out imports, multicore programming, printing the results, ...

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

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

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

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

Usage:

You have to install pypy (for Python 2!!!). The parallel python module isn't ported for Python 3. Then you have to install the parallel python module pp-1.6.4.zip. Extract it, cd into the folder and call pypy setup.py install.

Then you can call my program with

pypy you-do-the-math.py 15

It will automatically determine the number of cpu's. There may be some error messages after finishing the program, just ignore them. n = 16 should be possible on your machine.

Output:

Calculation for n = 15 took 2:50 minutes

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

Notes and ideas:

• I have a i7-4600m processor with 2 cores and 4 threads. It doesn't matter If I use 2 or 4 threads. The cpu-usage is 50% with 2 threads and 100% with 4 threads, but it still takes the same amount of time. I don't know why. I checked, that each thread only has to the the half amout of data, when there are 4 threads, checked the results, ...
• I use a lot of lists. Python isn't quite efficient in storing, I have to copy lots of lists, ... So I thought of using an integer instead. I could use the bits 00 (for 0) and 11 (for 1) in the vector A, and the bits 10 (for -1), 00 (for 0) and 01 (for 1) in the vector B. For the product of A and B, I would only have to calculate A & B and count the 01 and 10 blocks. Cycling can be done with shifting the vector and using masks, ... I actually implemented all this, you can find it in some of my older commits on Github. But it turned out, to be slower than with lists. I guess, pypy really optimizes list operations.

Answer 3

woolly bully - C++ - way too slow

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

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

Building the executable

It's a standalone C++11 source that compiles without warnings and runs smoothly on:

• Win7 & MSVC2013
• Win7 & MinGW - g++4.7
• Ubuntu & g++4.8 (in a VirtualBox VM with 2 CPUs allocated)

If you compile with g++, use : g++ -O3 -pthread -std=c++11
forgetting about -pthread will produce a nice and friendly core dump.

Optimizations

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(6ⁿ) to o(6ⁿ/(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.

How it works

The number of B values is huge (3ⁿ), so pre-computing them would require indecent amounts of memory, which would slow down computation and eventually exhaust available RAM.
Unfortunately, I could not find a simple way of enumerating the half set of optimized B values, so I resorted to coding a dedicated generator.

The mighty B generator was a lot of fun to code, though languages that support yield mechanisms would have allowed to program it in a much more elegant way.
What it does in a nutshell is consider the "skeleton" of a Bpre vector as a binary vector where 1s represent actual -1 or +1 values.
Among all these +1/-1 potential values, the first one is fixed to +1 (thus selecting one of the possible B/-B vector), and all remaining possible +1/-1 combinations are enumerated.
Lastly, a simple calibration system ensures each worker thread will process a range of values of approximately the same size.

A values are heavily filtered for regrouping in equiprobable chunks.
This is done in a pre-computing phase that brute-force examines all possible values.
This part has a neglectible O(2ⁿ) execution time and does not need to be optimized (the code being already unreadable enough as it is!).

To evaluate the inner product (which needs only be tested against zero), the -1 and 1 components of B are regrouped into binary vectors.
The inner product is null if (and only if) there is an equal number of +1s and -1s among the B values corresponding to non-zero A values.
This can be computed with simple masking and bit count operations, helped by std::bitset that will produce reasonably efficient bit count code without having to resort to ugly intrinsic instructions.

The work is equally divided among cores, with forced CPU affinity (every little bit helps, or so they say).

Example result

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

Performances

Multithreading should work perfectly on this, though only "true" cores will fully contribute to computation speed. My CPU has only 2 cores for 4 CPUs, and the gain over a single-threaded version is "only" about 3.5.

Compilers

An initial problem with multithreading led me to believe the GNU compilers were performing worse than Microsoft.

After more thorough testing, it appears g++ wins the day yet again, producing approximatively 30% faster code (the same ratio I noticed on two other computation-heavy projects).

Notably, the std::bitset library is implemented with dedicated bit count instructions by g++ 4.8, while MSVC 2013 uses only loops of conventional bit shifts.

As one could have expected, compiling in 32 or 64 bits makes no difference.

Further refinements

I noticed a few A groups producing the same probabilities after all the reduction operations, but I could not identify a pattern that would allow to regroup them.

Here are the pairs I noticed for 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