Linearize nested for loops

In this post, I shared a class named multi_index which basically does what you want, i.e.

for(auto m : multi_index(3,3,4))
{
    // now m[i] holds index of i-th loop
    // m[0] goes from 0 to 2
    // m[1] goes from 0 to 2
    // m[2] goes from 0 to 3
    std::cout<<m[0]<<" "<<m[1]<<" "<<m[2]<<std::endl;
}

However, this code is for "normal" loops only, where each dimension runs from 0 to some upper value.

In this post, I'll try to apply this to the antisymmetric case where m[i]<m[j] for i<j. The basic idea of the linked code remains the same, namely to create a class that holds the loop boundaries and provides an iterator which can be used with a range-based for loop. The only difference is that I use a std::vector instead of a std::array as the index array type:

#include <iostream>
#include <numeric>
#include <vector>

struct antisym_index_t
{
    int upper_index;
    int dim;
    antisym_index_t(int upper_index, int dim) : upper_index(upper_index), dim(dim) {}

    struct iterator
    {
        struct sentinel_t {};

        int upper_index;
        int dim;
        std::vector<int> index_array = {};
        bool _end = false;

        iterator(int upper_index, int dim) : upper_index(upper_index), dim(dim), index_array(dim)
        {
            std::iota(std::begin(index_array), std::end(index_array),0);
        }

        auto& operator++()
        {
            for (int i = dim-1;i >= 0;--i)
            {
                if (index_array[i] < upper_index - 1 - (dim-1-i))
                {
                    ++index_array[i];
                    for (int j = i+1;j < dim;++j)
                    {
                        index_array[j] = index_array[j-1]+1;
                    }                    
                    return *this;
                }
            }
            _end = true;
            return *this;
        }
        auto& operator*()
        {
            return index_array;
        }
        bool operator!=(sentinel_t) const
        {
            return !_end;
        }
    };

    auto begin() const
    {
        return iterator{ upper_index, dim };
    }
    auto end() const
    {
        return typename iterator::sentinel_t{};
    }
};

auto antisym_index(int upper_index, int dim)
{
    return antisym_index_t(upper_index, dim);
}

Note, however, that this code is untested so far (written on top of my head). You can use it as

for(auto m : antisym_index(5,3))
{
    // now m[i] holds index of i-th loop
    std::cout<<m[0]<<" "<<m[1]<<" "<<m[2]<<std::endl;
}

EDIT: by, now, I've tested and corrected the code, see here. Memo to myself: don't publish untested code.

EDIT2: by the way, this answers your question inside the question. It's not clear to me, how this should help with multitasking.

I don't have a full answer, but a solution for 2 loops. My sleep-deprived mind cannot generalize this to 3 loops but maybe someone else can.

In 2D the problem becomes figuring out the row and column index of a triangular matrix from a flattened index. This makes it easy to see that the end that "tapers off" is contained in the larger one. In ASCII art something like this:

      n
 ___________
|_          |
| |_        |
|   |_      |
|   | |_    |
|   |   |_  |
|___|_____|_|
  i   ^
      |
     binom(n-i, 2)

So, let's define

• n loop end index (number of matrix rows/columns)
• i outer loop counter range [0, n). As drawn: column index
• j inner loop counter range [0, i). As drawn: row index from bottom up
• a flattened loop counter range [0, binom(n, 2))

Then i can be computed from binom(n, 2) - binom(n-i, 2) = a. One round-trip through Wolfram Alpha gives us:

• i = trunc(-0.5 * sqrt((1 - 2 n)**2 - 8 a) + n - 0.5).

The truncation (=cast to int) "rounds down" to the last full column. So the row index j can be computed from as

• j = a - (binom(n, 2) - binom(n-i, 2))
• j = a - i*(-i + 2 n - 1) / 2

Yet another take on your problem. As said in the comments, what you are looking for is basically finding the successor and the unranking of combinations. For this I use the algorithms from the book 'Combinatorial algorithms' of Kreher and Stinson.

Here is the corresponding code consisting of the two functions next and unrank as well as a helper for the binomial coefficient which is required in the unranking function:

int binomial ( int n, int k )
{
    int mn = k;
    if ( n - k < mn )
    {
        mn = n - k;
    }

    if ( mn < 0 ) { return 0; }
    if ( mn == 0 ) { return 1; }

    int mx = k;
    if ( mx < n - k )
    {
        mx = n - k;
    }
    int value = mx + 1;

    for (int i = 2; i <= mn; ++i)
    {
        value = ( value * ( mx + i ) ) / i;
    }

    return value;
}         
        
auto unrank(int rank, int n, int k)
{    
    std::vector<int> t(k);
    
    int x = 1;
    for (int i = 0; i < k; ++i)
    {
        while (true)
        {
            int b = binomial ( n - x, k - i - 1);
            if (b > rank) break;
            rank -= b;
            ++x;
        }
    
        t[i] = x;
        ++x;
    }

  return t;
}

auto next(std::vector<int>& index, int n, int k)
{
    for (int i = k-1; i >= 0; --i)
    {
        if (index[i] < n - (k-1) + i)
        {
            ++index[i];
            for (int j = i+1; j < k; ++j)
            {
                index[j] = index[j-1]+1;
            }                    
            return true;
        }
    }
    return false;
}

The idea is then to generate the initial index-configuration from a given start address, and then compute the successor of this index (end-start) times. Here is an example:

int main()
{
    int n = 7;
    int k = 4;
    
    int start = 3;
    int end = 10;

    auto index = unrank(start,n,k);        
    auto print_index = [&]()
    {
        for(auto const& ind : index)
        {
            std::cout<<ind<<"  ";  
        }
        std::cout<<std::endl;
    };

    print_index();
    for(int i=start; i<end; ++i)
    {
        next(index, n, k);
        print_index();            
    }
}

which prints

1  2  3  7  
1  2  4  5  
1  2  4  6  
1  2  4  7  
1  2  5  6  
1  2  5  7  
1  2  6  7  
1  3  4  5 

And here is the Demo. Enjoy!

A third attempt:

I've taken your code, and at last got it to run properly (in python):

def get_k(n):
    total = 0
    for i in range(3, n):
        for j in range(i + 1, n):
            for k in range(j + 1, n):
                total += 1
            
    V = total // 2 # for 2 threads
    V_tmp = 0          
    for i in range(3, n):
        if(V_tmp > V):
            return i
        for j in range(i + 1, n):
            for k in range(j + 1, n):
                V_tmp += 1

def pseudo_thread(start, end, n):
    counter = 0

    for i in range(start, end):
        for j in range(i + 1, n):
            for k in range(j + 1, n):
                counter += 1
    print(counter)


n = 145
k = get_k(n)

pseudo_thread(3, k, n)
pseudo_thread(k, n, n)

This should finally give you a relatively good split. Even with n=145, we get 239260 and 227920 for our counter values. It's obviously not an elegant solution, and also not perfect, but it gives you the right answer without much reference to detailed mathematics.

Depending on how you want to parallelize, you could also use an atomic struct and implement the iteration with compare-and-swap operations. There is a 16 byte CAS on most platforms. Link with -latomic on GCC. Clang inlines the CAS calls if we ensure proper alignment.

#include <atomic>
#include <type_traits>
#include <cstdio>

/**
 * Index for a nested loop
 *
 * Index for loop in style
 * for(i = 0; i < n; ++i)
 *   for(j = 0; j < i; ++j)
 *     for(k = 0; k < j; ++k);
 *
 * The total number of iterations is binom(n, 3)
 *
 * Indices are int for two reasons:
 * 1. Keep overall size at or below 16 byte to allow atomic operations
 * 2. The total number of iterations reaches 2^64 at n ~ 4.8 million
 */
struct Index {
  int i, j, k;

  constexpr Index() noexcept
  : i(2), j(1), k(0)
  {}
  Index& operator++() noexcept
  {
    if(k + 1 < j) {
      ++k;
      return *this;
    }
    k = 0;
    if(j + 1 < i) {
      ++j;
      return *this;
    }
    j = 0;
    ++i;
    return *this;
  }
};

/**
 * Padds Index to power of 2 alignment up to 16 byte
 *
 * This improves atomic operation performance because it avoids
 * split-locks. Not sure if GCC's std::atomic makes actual use of this
 * but clang does.
 */
struct AlignedIndex
{
private:
  static constexpr std::size_t alignment =
    sizeof(Index) < 2 ? 1 :
    sizeof(Index) < 3 ? 2 :
    sizeof(Index) < 5 ? 4 :
    sizeof(Index) < 9 ? 8 :
    16;
public:
  union {
    std::aligned_storage<sizeof(Index), alignment>::type pod;
    Index index;
  };
  constexpr AlignedIndex() noexcept
  : index()
  {}
};

Index increment(std::atomic<AlignedIndex>& index) noexcept
{
  AlignedIndex last = index.load(std::memory_order_relaxed);
  AlignedIndex next;
  do {
    next = last;
    ++next.index;
  } while(! index.compare_exchange_weak(last, next, std::memory_order_relaxed));
  return last.index;
}

int main()
{
  std::atomic<AlignedIndex> index(AlignedIndex{});
  int n = 5;
  for(Index cur; (cur = increment(index)).i < n; ) {
    std::printf("%d %d %d\n", cur.i, cur.j, cur.k);
  }
}

Instead of iterating from 1..binom(n, 3), just iterate from 1..n^3 (conceptually the cartesian product of the set of numbers 1..n with itself 2x, as opposed to the combinations of 3 elements without repetition). In doing so, we can easily calculate i/j/k from M:

k = (M / N^0) % N = M % N
j = (M / N^1) % N
i = (M / N^2) % N = M / N^2

Sure, this will result in duplicates, but we're not going to skip the duplicates one by one. Once we reach a number where k>=j, we need to increment b by (N-k)*N^0 = N-k to cause it to "wrap around" to 0 again. Same goes for j>=i - increment b by (N-j)*N^1, to wrap around.

In doing so, we arrive back at just the original set of numbers. There's some overhead of division and modulus calculations, which can repeat up to once per variable (minus the first variable), so yes, there's some overhead, but its constant, for a constant number of variables.

Here is another solution based on Dillon Davis comments.

auto divide = [](float pos, int len) -> float {
    auto n = static_cast<float>(len);

    if (pos == 1) {
        return n;
    }

    if (pos == 0) {
        return 0;
    }

    // solve   x * (x - 1) * (x - 2) = n * (n - 1) * (n - 2) * pos   for x
    // https://en.wikipedia.org/wiki/Bisection_method

    float d = n * (n - 1) * (n - 2) * (1 - pos);

    auto f = [d](float x) {
        return std::pow(x, 3) - 3 * std::pow(x, 2) + 2 * x - d;
    };

    float a = 0;
    float b = n;
    float epsilon = 0.1f;

    float x = 0;

    while (std::abs(a - b) > epsilon) {
        x = (a + b) / 2;

        if (std::abs(f(x)) <= epsilon) {
            break;
        } else if (f(x) * f(a) < 0) {
            b = x;
        } else {
            a = x;
        }
    }

    return std::ceil(n - x);
};

It's pretty fast and with minimal overhead, but it's not as accurate as davidhigh's solution, which allows splitting 'work' into equal parts.

Example:

auto testRun = [](int begin, int end, int n) {
    int counter = 0;

    for (int i = begin; i < end; ++i) {
        for (int j = i + 1; j < n; ++j) {
            for (int k = j + 1; k < n; ++k) {
                ++counter;
            }
        }
    }

    std::cout << counter << "\n";
};

int n = 1200;
int ranges = 4;

for (int i = 0; i < ranges; ++i) {
    auto begin = static_cast<int>(divide((float) i / (float) ranges, n));
    auto end = static_cast<int>(divide((float) (i + 1) / (float) ranges, n));
    testRun(begin, end, n);
}

Output:

72035920
71897080
71619380
71728020