23 added 17 characters in body
source | link

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
source | link

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)
{
    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;
        }

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

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

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

// -----------------------------------------------------------------
// 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;
}
  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)
{
    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;
}
  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
source | link
  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
source | link
19 deleted 22 characters in body
source | link
18 added 431 characters in body
source | link
17 added 502 characters in body
source | link
16 deleted 1 character in body
source | link
15 deleted 1 character in body
source | link
14 added 1412 characters in body
source | link
13 added 1412 characters in body
source | link
12 added 1412 characters in body
source | link
11 deleted 1 character in body
source | link
10 edited body
source | link
9 added 61 characters in body
source | link
8 added 1451 characters in body
source | link
7 deleted 10 characters in body
source | link
6 added 178 characters in body
source | link
5 added 220 characters in body
source | link
4 added 1143 characters in body
source | link
3 deleted 4 characters in body
source | link
2 added 69 characters in body
source | link
1
source | link