Почему генератор случайных чисел не является потокобезопасным?

Я использовал функцию rand() для генерации псевдослучайных чисел между 0,1 для целей симуляции, но когда я решил сделать мой C++ код запущенным параллельно (через OpenMP), я заметил, что rand() безопасный, а также не очень однородный.

Поэтому я переключился на использование (так называемого) более однородного генератора, представленного во многих ответах по другим вопросам. Что выглядит так

double rnd(const double & min, const double & max) {
    static thread_local mt19937* generator = nullptr;
    if (!generator) generator = new mt19937(clock() + omp_get_thread_num());
    uniform_real_distribution<double> distribution(min, max);
    return fabs(distribution(*generator));
}

Но я видел много научных ошибок в моей оригинальной проблеме, которую я имитировал. Проблемы, которые были как против результатов rand() и против здравого смысла.

Поэтому я написал код для генерации случайных чисел 500k с этой функцией, вычислил их среднее значение и сделал это в 200 раз и построил результаты.

double SUM=0;
for(r=0; r<=10; r+=0.05){   
    #pragma omp parallel for ordered schedule(static)
    for(w=1; w<=500000; w++){   
        double a;
        a=rnd(0,1);
        SUM=SUM+a;
    } 
    SUM=SUM/w_max;
    ft<<r<<'\t'<<SUM<<'\n';
    SUM=0;
}

Мы знаем, что вместо 500k я мог бы делать это бесконечно, это должна быть простая строка со значением 0.5. Но с 500 тыс. Мы будем иметь колебания около 0,5.

При запуске кода с одним потоком результат допустим:

Но вот результат с двумя потоками:

3 темы:

4 темы:

У меня нет моего 8-поточного процессора прямо сейчас, но результаты даже стоили там.

Как вы можете видеть, они оба неравномерны и очень колеблются вокруг их среднего.

Так же ли это псевдослучайный генераторный поток - небезопасный?

Или я что-то ошибаюсь?

Ответ 1

Есть три замечания о вашем тестовом выходе, который я бы сделал:

Он имеет гораздо более сильную дисперсию, чем должен обеспечить хороший случайный источник. Вы сами это заметили по сравнению с результатами одного потока.
Вычисленное среднее уменьшается с количеством потоков и никогда не достигает исходного 0,5 (т.е. Это не только более высокая дисперсия, но также и среднее значение).
В данных есть временное отношение, особенно видимое в случае с 4 потоками.

Все это объясняется состоянием гонки, присутствующим в вашем коде: вы назначаете SUM из нескольких потоков. Приращение двойника не является атомной операцией (даже на x86, где вы, вероятно, получите атомные чтения и записи на регистрах). Два потока могут считывать текущее значение (например, 10), увеличивать его (например, как добавлять 0,5), а затем записывать значение обратно в память. Теперь у вас есть два потока, которые пишут 10.5 вместо правильного 11.

Чем больше потоков пытается написать SUM одновременно (без синхронизации), тем больше их изменений теряется. Этим объясняются все наблюдения:

Насколько жесткие нити расходятся друг с другом в каждом отдельном прогоне, определяет, сколько результатов будет потеряно, и это может варьироваться от прогона до запуска.
Среднее значение ниже с большим количеством рас (например, больше потоков), потому что больше результатов теряется. Вы никогда не сможете превзойти статистический 0,5 средний, потому что вы только теряете запись.
Когда потоки и планировщик "поселяются", дисперсия уменьшается. Это по той же причине, почему вы должны "разогревать" ваши тесты при бенчмаркинге.

Излишне говорить, что это неопределенное поведение. Он просто показывает доброкачественное поведение на процессорах x86, но это не то, что гарантирует вам стандарт C++. Насколько вам известно, отдельные байты двойника могут записываться разными потоками одновременно, что приводит к полному мусору.

Правильным решением будет добавление дублирования нити локально, а затем (с синхронизацией), добавляя их вместе в конце. OMP имеет положения о сокращении для этой конкретной цели.

Для интегральных типов вы можете использовать std::atomic<IntegralType>::fetch_add(). std::atomic<double> существует, но (до C++ 20) упомянутая функция (и другие) доступны только для интегральных типов.

Ответ 2

Проблема не в вашем RNG, а в вашем суммировании. На SUM есть просто состояние гонки. Чтобы исправить это, используйте сокращение, например, измените прагму на:

#pragma omp parallel for ordered schedule(static) reduction(+:SUM)

Обратите внимание, что использование thread_local с OpenMP технически не определено поведением. Это, вероятно, будет работать на практике, но взаимодействие между концепциями потоковой передачи OpenMP и С++ 11 недостаточно четко определено (см. Также этот вопрос). Таким образом, безопасная альтернатива OpenMP для вас будет:

static mt19937* generator = nullptr;
#pragma omp threadprivate(generator)