Задача N-тел: эффективное распараллеливание двойного цикла for

Очень распространенной проблемой в задаче N тел является использование двойного цикла для расчета взаимодействий между частицами. Рассматривая задачу N тел с n частицами, цикл можно записать как

for (i = 0, i < n; i++)
    for (j = i+1, j < n; j++)
        // calculate interaction

Мой вопрос о том, как можно распараллелить этот цикл, используя разные потоки. Цель состоит в том, чтобы каждый поток «в идеале» должен был рассчитывать одинаковое количество взаимодействий.

Моя идея заключалась в том, чтобы разделить внешний цикл, i-цикл, на разные интервалы, скажем, a_k=a(k), где k = 1,2,...,p, где p — количество потоков, которые мы хотим разделить проблему на.

Таким образом, цикл можно записать как

for (k = 1, k < p; k++)
    for (i = a(k), i < a(k+1); i++)
        for (j = i+1, j < n; j++)
            // calculate interaction

Где самый внешний цикл, k-цикл, должен быть распараллелен.

Поскольку количество взаимодействий в самом внутреннем цикле, j-цикле, равно n-(i+1), число взаимодействий, вычисленное каждым потоком, равно

\sum_{i=a(k)} ^{a(k+1)} n - (i+1)

Это означает, что нужно найти дискретную функцию a_k такую, что функционал

f[a_k] = \sum_{i=a(k )}^{a(k+1)} n - (i+1)

с граничными условиями a(1)=0 и a(p)=n является постоянным функционалом, таким образом заставляя число взаимодействий на каждая нить должна быть одинаковой.

Я пытался использовать разные «эвристики» (например, полином a_k, экспоненциальный, логарифмический), и до сих пор ни один из них не дал удовлетворительного ответа.Прямое решение этой проблемы для меня не очевидно.

Для малых p эту задачу можно отнести к «задачам минимизации мешка», где в основном каждое a_k является переменной для минимизации функции

f(a_1,a_2,a_3,...) = sum(|f [a_k] - n/p|^2)

Но как вы могли догадаться, это неэффективно (или даже не сходится) для более высоких значений p.

Кто-нибудь знает, как решить эту проблему?