C++0x threads, случайные числа и метод Монте-Карло

В основе моделирования методом Монте-Карло (МК), как известно, лежит генерация большого числа случайных чисел (СЧ) с заданными вероятностными характеристиками. По сути, метод состоит в проведении очень большого числа независимых экспериментов, а т.к. они независимые, значит их можно выполнять параллельно, и распараллеливание алгоритмов, основанных на методе МК, должно быть очень эффективно.
Я решил проверить это на классическом алгоритме для расчета числа π. Алгоритм состоит в следующем. Берется квадрат со стороной единичной длины и в нем рисуется сектор круга единичного радиуса с центром в одном из углов.



Затем в этот квадрат начинают многократно бросать дротик. Отношение количества попаданий в сектор круга к общему количеству бросков будет приблизительно соответствовать отношению площадей четверти круга и квадрата. Чем больше бросков, тем выше точность. Если площадь четверти круга рассчитывается по формуле S = ¼πR², а площадь квадрата — S' = a², где a — длина стороны квадрата, а R — радиус круга, то при a = R = 1, π = 4(S/S').
Броски дротика можно моделировать путем генерации равномерно распределенных случайных чисел (x, y) в интервале [0; 1], тогда их принадлежность сектору круга будет определятся по формуле x² + y² ≥ 1.

Последовательная версия.

#include <iostream>
#include <cstdlib> // тут генератор СЧ
#include <time.h> // системное время для инициализации генератора

using namespace std;

const long N = 1000000000; // количество бросков в мешень

int main()
{
srand(time(0)); // инициализация генератора
long hit = 0; // количество попаданий
double x, y; // координаты попаданий, оказывается, так на 2 с быстрее,
// чем если объявлять непосредственно в цикле
for(long i = 0; i < N; ++i) {
    x = double(rand()) / RAND_MAX; // получаем новую пару координат
    y = double(rand()) / RAND_MAX; // в промежутке 0..1
    if (x * x + y * y < 1) ++hit; // +1, если попали
}
double pi = 4 * (double(hit) / N); // т.к. бросали только в 1/4 круга
cout << "pi = " << pi << endl;
return 0;
}

Эта программа позволяет получить значение π = 3.14157 для 1.000.000.000 бросков и работает у меня на 4-ядерном Athlon 2,7 МГц 1 минуту и 4 секунды. При этом, загрузка процессора составляет около 25%, т.к. используется только 1 ядро. Естественно, возникает желание ее ускорить за счет распараллеливания: ничего не мешает нам бросать по 4 дротика за раз. Ядер-то 4.
Но это оказывается не так просто: генератор СЧ, который, по сути, является генератором псевдо-СЧ, представляет собой функцию, рассчитывающую каждое новое число на основе предыдущего. Предыдущее число принято называть seed, и оно, фактически, является скрытым состоянием функции-генератора. Это означает, что такая функция не потокобезопасна, и в многопоточной среде ею пользоваться нельзя.
Способом решения этой проблемы является использование функций, которые получают seed в качестве параметра. К таким функциям относится rand_r(), однако тут возникает следующая проблема. Эта функция возвращает СЧ типа int, которых на платформе x86 может быть всего 2³². Это значит, что в лучшем случае через 4.294.967.295 чисел последовательность чисел начнет повторяться. Для 1.000.000.000 итераций в последовательной реализации это будет означать, что все СЧ будут разными (на каждый бросок 2 итерации, т.е. 2.000.000.000 чисел). Для 4-х потоков нужно 4 генератора, а т.к. алгоритм генерации одинаков, то при близко расположенных seed'ах в параллельных потоках мы будем получать большой процент одинаковых чисел, что приведет к смещению мат. ожидания генератора.
В приведенной ниже реализации (конечно же, с использованием потоков C++0x, кстати, компилировать ее нужно с ключами -std=c++0x и -pthread) я постарался разнести как можно дальше seed'ы каждого потока (time(0) * i), чтобы уменьшить наложение последовательностей одинаковых значений, но значение π колебалось от 3.1234 до 3.1478. Зато время вычисления составило 15 с, что в 4,4 раза быстрее последовательного варианта. Кроме rand_r() с 32-битным seed'ом, есть еще ряд функций с 48-битным (закомментированный вариант). Использование чисел с большей точностью приводит к снижению производительности — до 22 с, всего в 3 раза быстрее последовательной версии, но зато точность остается на уровне π = 3.14149 (ее, вероятно, можно еще улучшить за счет увеличения количества бросков).

Параллельная версия.

#include <iostream>
#include <cstdlib>
#include <thread> // тут работа с потоками
#include <time.h>

using namespace std;

const long N = 1000000000;
const int P = 4; // количество процессоров

int main()
{
thread* t[P];
mutex m; // надо для синхронизации, когда будем собирать общую сумму
long hit = 0;

for(int i = 0; i < P; ++i) { // запускаем P потоков
    t[i] = new thread([&hit, &m, N, i]() { // это лямбда-функция
        long hit_ = 0; // локальный счетчик попаданий
        double x, y;
        unsigned int r = time(0) * i; // i чтобы гарантировать разные seed
                                      // во всех потоках
// --- это вариант с nrand48() ---------------
//            unsigned short seed[3];
//            unsigned int r = time(0) * i;
//            copy(seed, seed + 1, &r);
// -------------------------------------------
        for(long i = 0; i < N / P; ++i) { // N / P в P потоках будет ~ N
            x = double(r = rand_r(&r)) / RAND_MAX;
            y = double(r = rand_r(&r)) / RAND_MAX;
// --- это вариант с nrand48() ---------------
//                x = erand48(seed);
//                y = erand48(seed);
// -------------------------------------------
            if (x * x + y * y < 1) ++hit_;
        }
        lock_guard<mutex> lg(m); // блокировка до конца блока, чтобы не было
        hit += hit_;             // конкуренции за hit
    });
}

for(int i = 0; i <>join();
}
double pi = 4 * (double(hit) / N);
cout << "pi = " << pi << endl;
return 0;
}

Еще о C++0x:

• C++0x: lambda functions

Еще о параллелизме:

• Синхронизация в потоках POSIX: семафоры, мютексы, спин-локи и атомарные операции
• CUDA + Ubuntu 9.10