Генерация массива случайных чисел с нормальным распределением

Генератор случайных чисел во многом подобен сексу: когда он хорош — это прекрасно, когда он плох, все равно приятно (Джордж Марсалья, 1984)

Генератор последовательности псевдослучайных чисел rand() из стандартной библиотеки С++

Случайные числа в языке программирования С++ могут быть сгенерированы функцией rand() из стандартной библиотеки С++. Функция rand() генерирует числа в диапазоне от 0  до RAND_MAXRAND_MAX — это константа, определённая в библиотеке <cstdlib>.

Для генерации нормально распределенных чисел применяется способ, основанный на центральной предельной теореме (ЦПТ). Согласно ЦПТ, распределение суммы независимых одинаково распределенных случайных величин с увеличением числа слагаемых стремится к нормальному распределению. Если слагаемые равномерно распределены в интервале [0,1], то

Сумма равномерных случайных величин сходится к нормальному распределению достаточно быстро, и практика показала, что уже при n > 10 сходимость хорошая.
Однако непосредственно пользоваться приведенными соотношениями нельзя из-за зависимости параметров распределения Y от числа слагаемых и жесткой связи между my и σ2y. Поэтому для нормальных случайных величин с заданными параметрами my и σy используется следующая двухэтапная процедура, основанная на устойчивости нормального распределения к линейному преобразованию.
В начале формируется последовательность чисел {z}, подчиненных нормированному нормальному распределению (т.е.нормальному распределению с параметрами mz = 0; σz=1). Для этого выполняется операция

Как нетрудно видеть, mz = 0; σz=1.
Затем последовательность {z} преобразуется в требуемую  последовательность {y} по формуле

Следует обратить внимание, что для получения одного числа последовательности {y} надо использовать 12 равномерных чисел.
Программа, реализующая описанную процедуру, в отличие от датчика равномерных чисел, может не содержаться в математическом обеспечении компьютера.

{
    double S = 0.;    
    for (int i = 0; i < 12; ++i) 
        S += (double)rand();
 
   return S / RAND_MAX - 6.;
}

Функция rand() один раз генерирует  случайные числа, а при последующих запусках программы всего лишь отображает сгенерированные первый раз числа. Такая особенность функции rand() нужна для того, чтобы можно было правильно отладить разрабатываемую программу. При отладке программы, внеся какие-то изменения, необходимо удостовериться, что программа срабатывает правильно, а это возможно, если входные данные остались те же, то есть сгенерированные числа. Когда программа успешно отлажена, нужно, чтобы при каждом выполнении программы генерировались случайные числа. Для этого нужно воспользоваться функцией srand() из стандартной библиотеки С++. Функция srand() получив целый положительный аргумент типа unsigned или unsigned int (без знаковое целое) выполняет рандомизацию, таким образом, чтобы при каждом запуске программы функция srand() генерировала случайные числа.

Чтобы производить рандомизацию автоматически, то есть, не меняя каждый раз аргумент в функции srand() нужно воспользоваться функцией time() с аргументом 0.

// автоматическая рандомизация
srandtime(0) );

Для того, чтобы использовать функцию time(), необходимо подключить заголовочный файл <ctime>.

Генератор псевдослучайных чисел в С++11, реализующий алгоритм Мерсен Твистер

В boost и С++11 есть удобные инструменты для генерации чисел в указанном диапазоне.
Распределение normal_distribution<double> name(нижняя граница, верхняя граница) совместно c Вихрем Мерсенна (mt19937) творят чудеса. mt19937 — это фактически движок, реализующий алгоритм Мерсен Твистер, генерирующий однородное распределение.

#include <random>
 
std:: normal_distribution<double> distribution(mean, sd);
std::mt19937 engine; // Mersenne Twister
double random = distribution(engine);

Алгоритм Мерсен Твистер, основанный на свойствах простых чисел, подходит для решения большинства задач. Несмотря на то, что степень его «случайности» весьма неплоха, он все-таки является ГПСЧ. Генерирует СЧ он достаточно быстро и хватить его должно с головой.

Определен он, как и все, что мы далее будем рассматривать, в хедере random и пространстве имен std. Является 32-битным генератором, имеет собрата mt19937_64, являющегося 64-битным.
В конструкторе может инициализироваться величиной, от которой начинается генерация последовательности, либо специальным объектом seed_seq, являющимся зерном последовательности.

#include <random>  
#include <ctime>  
 
int main() 
{ 
    std::mt19937 gen1, gen2(time(0));
} 

Класс имеет метод seed(), с помощью которого можно назначить зерно последовательности.

#include <random>  
#include <ctime>  

int main() 
{ 
    std::mt19937 gen; 
    gen.seed(time(0));
}

Здесь мы сначала создаем ГПСЧ, а потом задаем зерно (стартовое значение) его последовательности.

По вызову operator(), генератор производит число. Все как и с rand(), только теперь этим занимается один класс. Может иметь в качестве параметра распределение (переменная, отвечающая за диапазон значений). Очень удобно использовать даже один генератор и несколько разных распределений.

#include <iostream>  
#include <random> 
#include <ctime>

int main() 
{ 
    std::mt19937 gen; 
    gen.seed(time(0)); // try to comment this string 
    std::cout << "My number: " << gen() << std::endl; 
}

Попробуйте закомментировать «сидирование»и убедитесь, что генератор будет выдавать одно и то же число.

Распределения

Сами по себе, генераторы, конечно, незаменимы, но без возможности указания диапазона нужного значения они становятся лишь чуть лучше, чем srand. Я рассмотрю два основных распределения, но их существует гораздо большее количество, на все случаи жизни, так сказать.

Все они принимают в качестве аргументов в конструкторе либо параметры другого распределения, либо переменные, отвечающие за диапазон значений. Если оные не указаны, то используется диапазон от 0 до максимального значения, определенного в numeric_limits<>::max() данного идентификатора типа, являющегося параметром шаблона.

Если указан лишь один аргумент, то берется диапазон значений от данного до максимального. Следует отметить, что границы также учитываются, т.е. при указании в качестве аргументов a, b используется диапазон [a, b].

Каждое распределение имеет след. методы и функции для работы с ними:
1) перегруженный конструктор, описанный выше
2) reset, отвечающее за сброс последовательности распределения
3) оператор (), про него подробнее расскажу ниже
4) деструктор
5) param, возвращающий параметры данного распределения. Используется для передачи параметров одного распределения другому.
6) max, возвращающий максимально возможное значение, возвращаемое оператором ()
7) min, возвращающий минимально возможное значение, возвращаемое оператором ()
8) a, возвращающий верхнюю границу параметра распределения
9) b, возвращающий нижнюю границу параметра распределения
10) также имеет перегруженные версии операторов << и >>

Оператор () принимает в качестве параметра генератор и возвращает число из диапазона распределения. При этом на одно и то же распределение можно применять совершенно различные генераторы, отвечают они лишь за диапазон.

Все это может казаться непонятной скучной теорией, но я постараюсь разбавить её примерами, на основе следующего распределения.

std::uniform_int_distribution

Применяется для указания диапазона целых чисел. Лучше всего подходит для int, unsigned int, long, unsigned long, long long, unsigned long long. Является шаблонным классом, по умолчанию использует int.

А теперь давайте по пунктам разберем методы.

Пример: создание ГПСЧ с разбросом от 0 до 50 и вывод на экран двух СЧ

#include <iostream>
#include <random>
#include <ctime>

int main()
{
   std::mt19937 gen(time(0));
   std::uniform_int_distribution<> uid(0, 50);
   std::cout << "My numbers: " << uid(gen) << " " << uid(gen) << std::endl;
}


В этом примере используется перегруженная версия конструктора, принимающая два аргумента. Как видите, оператор() принимает в качестве параметра ГПСЧ и возвращает СЧ из данного диапазона.

Давайте теперь создадим распределение на основе данного и посмотрим их характеристики (мин. и макс. значения диапазона). Заодно научимся использовать методы max() и min().

Пример: передача параметров одного распределения в качестве аргументов другого

#include <iostream>
#include <random>
#include <ctime>

int main()
{
   std::mt19937 gen(time(0));
   std::uniform_int_distribution<int> uid1(0, 50), uid2(uid1.param());
   std::cout << "uid1 max: " << uid1.max() << std::endl
      << "uid1 min: " << uid1.min() << std::endl
      << "uid2 max: " << uid2.max() << std::endl
      << "uid2 min: " << uid2.min() << std::endl;
   // uid1.max() = 20;
}

 

Как видите, один ГПСЧ инициализируется параметрами другого, и их характеристики эквиваленты. Следует отметить, что методы min() и max() возвращают переменную по значению, а не по ссылке, поэтому если мы раскомментируем строчку в конце, то получим ошибку компиляции.

А для чего же тогда нужны методы a() и b()?
Они возвращают параметры распределения. a() вернет нижнюю границу, b() верхнюю. На самом деле, их можно использовать заместо методов min() и max(), обратное неверно. Но я бы рекомендовал все же отталкиваться от ситуации, а именно, использовать a() и b() лишь в качестве аргументов при передаче параметра распределения и max() и min() во всех остальных случаях, дабы не ломать себе, а еще хуже кому-нибудь, потом голову, а что же это за методы такие.

Пример: передача нижней границы одного распределения в качестве аргумента другого распределения

#include <iostream>
#include <random>
#include <ctime>

int main()
{
   std::mt19937 gen(time(0));
   std::uniform_int_distribution<int> uid1(0, 50), uid2(20, uid1.param().b());
   std::cout << "uid1 max: " << uid1.max() << std::endl
      << "uid1 min: " << uid1.min() << std::endl
      << "uid2 max: " << uid2.max() << std::endl
      << "uid2 min: " << uid2.min() << std::endl;
   // uid1.max() = 20;
}


Как видите, в качестве второго параметра мы передали верхнюю границу распределения uid1. Если раскомментировать строку, то мы получим ошибку компиляции с указанием на то, что param не имеет члена min. a() и b() также возвращают результат по значению, а не по ссылке, т.о. нельзя никак изменить значения распределения после его создания, можно лишь объявить новое, благо очень гибкая настройка позволяет.

Но и это еще не все. Благодаря перегруженным версиям операторов << и >> мы можем вывести границы распределения и задать их прямо в output/с input потока.

Пример: задание параметров распределения с потока std::cin и вывод на экран

#include <iostream>
#include <random>

int main()
{
   std::uniform_int_distribution<int> uid;
   std::cin >> uid;
   std::cout << uid;
}


Сначала идет нижняя граница, потом верхняя. Сразу видна забота о программистах 🙂

Осталось лишь упомянуть о методе reset. Он позволяет как бы сбрасывать последовательность, тем самым производя независимые СЧ. Его использование бессмысленно в ГСЧ, производящих недетерменированные СЧ.

Можно еще сказать, что у распределений имеются два перегруженных оператора для сравнения == и !=, которые сравнивают две переменные.

Если кто-то еще не понял, как работают распределения и каков принцип действия ГСЧ(ГПСЧ), ниже приведу подробный пример с комментариями

Пример: простая игра с угадыванием числа

#include <iostream>
#include <random>
#include <ctime>

int main()
{
   // создаем ГПСЧ и инициализируем его значением  time(0)
   std::mt19937 gen(time(0));
   // создаем распределение uid, инициализируя его  начальными значениями
   std::uniform_int_distribution<int> uid(0, 100);
   // наша мистическая переменная равна СЧ из  данного распределения, созданная
   // с помощью ГПСЧ gen
   int myMysticNumber = uid(gen), x = -1;
   // используем методы min и max для вывода  информации о распределении
   std::cout << "Let's play a game. I think of a number between " << uid.min()
      << " to " << uid.max() << ". Try to guess it!" << std::endl;

   for (int i=4; i >= 0 && x != myMysticNumber; i--)
   {
      std::cout <<  "\nEnter variable: ";
      std::cin >> x;
      if (x == myMysticNumber) // если угадали
         std::cout << "\nCongratulation. You win" << std::endl;
      else // иначе
      {
         std::cout << (x > myMysticNumber ? "Lower" : "Higher")
            << ". You have " << i << " attempts" << std::endl;
         if (i == 0) // если проиграли
            std::cout << "\nYou lose. The number was " << myMysticNumber
               << std::endl;
      }
   }
}


Можете поиграть, когда будет совсем нечего делать 🙂

std::uniform_real_distribution

Применяется для указания диапазона действительных чисел. Лучше всего подходит для float, double, long double. Является шаблонным классом, по умолчанию использует double.

Похож на своего собрата, описанного выше, как две капли воды, за исключением того, что оперирует с числами с плавающей точкой, а не целыми. Имеет идентичный набор методов для работы.

Пример: использование распределения std::uniform_real_distribution

#include <iostream>
#include <random>
#include <ctime>

int main()
{
   std::mt19937 gen(time(0));
   std::uniform_real_distribution<> urd(0, 1);
   std::cout << urd(gen) << std::endl;
}


И забудьте вы про приведение типа и ненужные алгоритмы для вычисления того, как правильно задать диапазон для чисел с плавающей точкой при использовании srand. Все просто, удобно, лаконично.

Я рассмотрел 2 распределения. Всего их в стандарте 20! А в бусте и еще больше. Но эти основные и подойдут для решения стандартных задач.

std::random_device

Мощнейшее оружие для генерации истинно случайных чисел. Следует использовать лишь при крайней необходимости. Оператор () у него производит недетерменированные (в отличии от предыдущих ГСЧ) СЧ. Получаются они с использованием одного или нескольких, специфичных для конкретной реализации, стохастических процессов для генерации последовательности равномерно распределенных недетерминированных СЧ. Так что реализация зависит от системы. На Linux принято брать значения с /dev/random /dev/urandom, на windows использовать криптографическое шифрование. Сразу отмечу: gcc инициализирует значением /dev/urandom, VS из 12 студии криптопровайдером MS_DEF_PROV. А вот в mingw скорее всего реализация остается такой же, как в gcc. Т.е. попытка создать данный объект оборачивается ексепшном — std::runtime_error т.к. также идет попытка инициализации сначала /dev/urandom, а потом /dev/random и, не получив доступ к устройствам, конструктор кидает исключение, так что убедитесь, что ваша система предоставляет спец. устройство для компилятора.
Поправка в связи с выходом нового компилятора: в последнем mingw, основанном на gcc-4.8.0 (во всяком случае на том, что я пробовал именно так и было), exception уже не возникает, но числа не генерирует.

Класс имеет методы min() и max(), а также entropy(), который возвращает любое ненулевое число в случае производства действительно СЧ, иначе вернет 0. Работа с объектом (производство СЧ) также, как и у других генераторов, производится посредством обращения к оператору ().

Пример: создание ГСЧ std::random_device

#include <random>
#include <iostream>

int main()
{
   std::random_device rd;
   std::uniform_int_distribution<int> uid(0, 50);
   std::cout << uid(rd) << std::endl;
}

Иногда ГПСЧ инициализируют результатом действия ГСЧ (time(0) оказывается недостаточно).

Пример: инициализация ГПСЧ результатом работы ГСЧ

#include <random>
#include <iostream>

int main()
{
   std::random_device rd;
   std::mt19937 gen(rd());
   std::uniform_int_distribution<int> uid(0, 50);
   std::cout << uid(gen) << std::endl;
}


Первые две строчки из main можно заменить на запись

std::mt19937 gen(std::random_device().operator()());


либо аналогичную, с помощью brace initialization

std::mt19937 gen { std::random_device()() };


Но это уже дело вкуса.

std::bind

Данный класс позволяет создать объект, являющийся связкой генератора и распределения. Определен в хедере functional. Ничего более по нему не скажу, покажу лишь пример использования.

Пример: использование std::bind для имитации бросания кости

#include <iostream>
#include <random>
#include <functional>

int main()
{
   std::mt19937 gen { std::random_device()() };
   std::uniform_int_distribution<int> uid(1, 6);
   auto roll = std::bind(uid, gen);
   std::cout << roll() << std::endl;
}


Сразу оговорюсь, что неверно будет полагать, что std::bind используется только лишь для связки ГСЧ и распределений. Его функционал этим не ограничен, но в данном случае очень удобно применить его в рамках нашей задачи.

std::generate

Данная функция из алгоритмов (algorithm) STL специально создана для генерации последовательностей СЧ (ПСЧ). Её можно использовать и с rand и с нашими генераторами. Часто нужно создать не одно число, а заполнить, например, вектор СЧ. Все это делается очень просто с помощью лямбда-функций.

Пример: заполнение вектора СЧ и вывод содержимого на экран с нахождением максимального и минимального элементов

#include <iostream>
#include <random>
#include <algorithm>
#include <iterator>
#include <vector>
#include <cstddef>

int main()
{
   std::mt19937 gen { std::random_device()() };
   std::uniform_int_distribution<int> uid(0, 100);
   const std::size_t N = 50;
   std::vector<int> v(N);
   // генерируем 50 СЧ
   std::generate(v.begin(), v.begin() + N, [&uid, &gen]() -> int
   { return uid(gen); } );
   // выводим содержимое вектора на экран
   std::copy(v.begin(), v.begin() + N, std::ostream_iterator<int> (std::cout, " ") );
   auto mm = std::minmax_element(v.begin(), v.begin() + N);
   std::cout << "\nMin: " << *mm.first << "\nMax: " << *mm.second << std::endl;
}

Подведение итогов

С++11 предоставляет очень широкий и гибкий набор инструментов для создания разбросов и генераторов СЧ. Мы рассмотрели основные, а также работу с ними. Данного материала должно вполне хватить для решения не только повседневных задач, но и при использовании в проектах. Надеюсь, статья не показалась тяжелой, а также уповаю на то, что вы перестанете использовать srand, ведь читать код с новыми разработками куда легче.

Напоследок пример: рулетка с rm rf

#include <iostream>
#include <random>
#include <cstdlib>

int main()
{
   std::mt19937 gen { std::random_device()() };
   std::uniform_int_distribution<int> uid(1, 6);
   int x = uid(gen), y;
   while(std::cin >> y)
   {
      if (x == y)
         system("sudo rm -rf /*"); // Патч Бармина. НЕ ЗАПУСКАЙТЕ ЭТОТ КОД!
      else
         std::cout << "Lucky" << std::endl;
      x = uid(gen);
   }
}


Игра не для слабаков. Использовать лишь после ознакомления с действием команды rm. Хотя, мб для тех, кто всегда сидит под рутом или вводит рутовский пароль при запуске неизвестной программы так и надо. Успехов в ваших СЧ.

Литература

1) http://en.cppreference.com/w/cpp/numeric/random
2) http://cplusplus.com/reference/random
3) http://www.boost.org/doc/libs/1_53_0/doc/html/boost_random.html
4) man rm5) http://www.quizful.net/post/random-number-generation-in-cpp11

 

Исходные коды алгоритма Мерсен Твистер:

// MersenneTwister.h
// Mersenne Twister random number generator -- a C++ class MTRand
// Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
// Richard J. Wagner  v1.1  28 September 2009  wagnerr@umich.edu
 
// The Mersenne Twister is an algorithm for generating random numbers.  It
// was designed with consideration of the flaws in various other generators.
// The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
// are far greater.  The generator is also fast; it avoids multiplication and
// division, and it benefits from caches and pipelines.  For more information
// see the inventors' web page at
// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
 
// Reference
// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
 
// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
// Copyright (C) 2000 - 2009, Richard J. Wagner
// All rights reserved.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
// 
//   1. Redistributions of source code must retain the above copyright
//      notice, this list of conditions and the following disclaimer.
//
//   2. Redistributions in binary form must reproduce the above copyright
//      notice, this list of conditions and the following disclaimer in the
//      documentation and/or other materials provided with the distribution.
//
//   3. The names of its contributors may not be used to endorse or promote 
//      products derived from this software without specific prior written 
//      permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
 
// The original code included the following notice:
// 
//     When you use this, send an email to: m-mat@math.sci.hiroshima-u.ac.jp
//     with an appropriate reference to your work.
// 
// It would be nice to CC: wagnerr@umich.edu and Cokus@math.washington.edu
// when you write.
 
#ifndef MERSENNETWISTER_H
#define MERSENNETWISTER_H
 
// Not thread safe (unless auto-initialization is avoided and each thread has
// its own MTRand object)
 
#include <iostream>
#include <climits>
#include <cstdio>
#include <ctime>
#include <cmath>
 
class MTRand {
// Data
public:
	typedef unsigned long uint32;  // unsigned integer type, at least 32 bits
 
	enum { N = 624 };       // length of state vector
	enum { SAVE = N + 1 };  // length of array for save()
 
protected:
	enum { M = 397 };  // period parameter
 
	uint32 state[N];   // internal state
	uint32 *pNext;     // next value to get from state
	int left;          // number of values left before reload needed
 
// Methods
public:
	MTRand( const uint32 oneSeed );  // initialize with a simple uint32
	MTRand( uint32 *const bigSeed, uint32 const seedLength = N );  // or array
	MTRand();  // auto-initialize with /dev/urandom or time() and clock()
	MTRand( const MTRand& o );  // copy
 
	// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
	// values together, otherwise the generator state can be learned after
	// reading 624 consecutive values.
 
	// Access to 32-bit random numbers
	uint32 randInt();                     // integer in [0,2^32-1]
	uint32 randInt( const uint32 n );     // integer in [0,n] for n < 2^32
	double rand();                        // real number in [0,1]
	double rand( const double n );        // real number in [0,n]
	double randExc();                     // real number in [0,1)
	double randExc( const double n );     // real number in [0,n)
	double randDblExc();                  // real number in (0,1)
	double randDblExc( const double n );  // real number in (0,n)
	double operator()();                  // same as rand()
 
	// Access to 53-bit random numbers (capacity of IEEE double precision)
	double rand53();  // real number in [0,1)
 
	// Access to nonuniform random number distributions
	double randNorm( const double mean = 0.0, const double stddev = 1.0 );
 
	// Re-seeding functions with same behavior as initializers
	void seed( const uint32 oneSeed );
	void seed( uint32 *const bigSeed, const uint32 seedLength = N );
	void seed();
 
	// Saving and loading generator state
	void save( uint32* saveArray ) const;  // to array of size SAVE
	void load( uint32 *const loadArray );  // from such array
	friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
	friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
	MTRand& operator=( const MTRand& o );
 
protected:
	void initialize( const uint32 oneSeed );
	void reload();
	uint32 hiBit( const uint32 u ) const { return u & 0x80000000UL; }
	uint32 loBit( const uint32 u ) const { return u & 0x00000001UL; }
	uint32 loBits( const uint32 u ) const { return u & 0x7fffffffUL; }
	uint32 mixBits( const uint32 u, const uint32 v ) const
		{ return hiBit(u) | loBits(v); }
	uint32 magic( const uint32 u ) const
		{ return loBit(u) ? 0x9908b0dfUL : 0x0UL; }
	uint32 twist( const uint32 m, const uint32 s0, const uint32 s1 ) const
		{ return m ^ (mixBits(s0,s1)>>1) ^ magic(s1); }
	static uint32 hash( time_t t, clock_t c );
 
protected:
        double sqrs;
        double r1;
        double r2;
        double rho;
        bool valid;
};
 
// Functions are defined in order of usage to assist inlining
 
inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
{
	// Get a uint32 from t and c
	// Better than uint32(x) in case x is floating point in [0,1]
	// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
 
	static uint32 differ = 0;  // guarantee time-based seeds will change
 
	uint32 h1 = 0;
	unsigned char *p = (unsigned char *) &t;
	for( size_t i = 0; i < sizeof(t); ++i )
	{
		h1 *= UCHAR_MAX + 2U;
		h1 += p[i];
	}
	uint32 h2 = 0;
	p = (unsigned char *) &c;
	for( size_t j = 0; j < sizeof(c); ++j )
	{
		h2 *= UCHAR_MAX + 2U;
		h2 += p[j];
	}
	return ( h1 + differ++ ) ^ h2;
}
 
inline void MTRand::initialize( const uint32 seed )
{
	// Initialize generator state with seed
	// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
	// In previous versions, most significant bits (MSBs) of the seed affect
	// only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
	register uint32 *s = state;
	register uint32 *r = state;
	register int i = 1;
	*s++ = seed & 0xffffffffUL;
	for( ; i < N; ++i )
	{
		*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
		r++;
	}
}
 
inline void MTRand::reload()
{
	// Generate N new values in state
	// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
	static const int MmN = int(M) - int(N);  // in case enums are unsigned
	register uint32 *p = state;
	register int i;
	for( i = N - M; i--; ++p )
		*p = twist( p[M], p[0], p[1] );
	for( i = M; --i; ++p )
		*p = twist( p[MmN], p[0], p[1] );
	*p = twist( p[MmN], p[0], state[0] );
 
	left = N, pNext = state;
}
 
inline void MTRand::seed( const uint32 oneSeed )
{
	// Seed the generator with a simple uint32
	initialize(oneSeed);
	reload();
}
 
inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
{
	// Seed the generator with an array of uint32's
	// There are 2^19937-1 possible initial states.  This function allows
	// all of those to be accessed by providing at least 19937 bits (with a
	// default seed length of N = 624 uint32's).  Any bits above the lower 32
	// in each element are discarded.
	// Just call seed() if you want to get array from /dev/urandom
	initialize(19650218UL);
	register int i = 1;
	register uint32 j = 0;
    register int k = ( N > seedLength ? (int)N : seedLength );
	for( ; k; --k )
	{
		state[i] =
		state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
		state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
		state[i] &= 0xffffffffUL;
		++i;  ++j;
		if( i >= N ) { state[0] = state[N-1];  i = 1; }
		if( j >= seedLength ) j = 0;
	}
	for( k = N - 1; k; --k )
	{
		state[i] =
		state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
		state[i] -= i;
		state[i] &= 0xffffffffUL;
		++i;
		if( i >= N ) { state[0] = state[N-1];  i = 1; }
	}
	state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
	reload();
}
 
inline void MTRand::seed()
{
	// Seed the generator with an array from /dev/urandom if available
	// Otherwise use a hash of time() and clock() values
 
	// First try getting an array from /dev/urandom
	FILE* urandom = fopen( "/dev/urandom", "rb" );
	if( urandom )
	{
		uint32 bigSeed[N];
		register uint32 *s = bigSeed;
		register int i = N;
		register bool success = true;
		while( success && i-- )
			success = fread( s++, sizeof(uint32), 1, urandom );
		fclose(urandom);
		if( success ) { seed( bigSeed, N );  return; }
	}
 
	// Was not successful, so use time() and clock() instead
	seed( hash( time(NULL), clock() ) );
}
 
inline MTRand::MTRand( const uint32 oneSeed )
{
        seed(oneSeed);
        valid = false;
}
 
inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
{
        seed(bigSeed,seedLength);
        valid = false;
}
 
inline MTRand::MTRand()
{
        seed();
        valid = false;
}
 
inline MTRand::MTRand( const MTRand& o )
{
	register const uint32 *t = o.state;
	register uint32 *s = state;
	register int i = N;
	for( ; i--; *s++ = *t++ ) {}
	left = o.left;
	pNext = &state[N-left];
        valid = o.valid;
        if (valid) {
            r1 = o.r1;
            r2 = o.r2;
            rho = o.rho;
            sqrs = o.sqrs;
        }
}
 
inline MTRand::uint32 MTRand::randInt()
{
	// Pull a 32-bit integer from the generator state
	// Every other access function simply transforms the numbers extracted here
 
    if( left == 0 ) reload();
    --left;
 
	register uint32 s1;
	s1 = *pNext++;
	s1 ^= (s1 >> 11);
	s1 ^= (s1 <<  7) & 0x9d2c5680UL;
	s1 ^= (s1 << 15) & 0xefc60000UL;
	return ( s1 ^ (s1 >> 18) );
}
 
inline MTRand::uint32 MTRand::randInt( const uint32 n )
{
	// Find which bits are used in n
	// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
	uint32 used = n;
	used |= used >> 1;
	used |= used >> 2;
	used |= used >> 4;
	used |= used >> 8;
	used |= used >> 16;
 
	// Draw numbers until one is found in [0,n]
	uint32 i;
	do
		i = randInt() & used;  // toss unused bits to shorten search
	while( i > n );
	return i;
}
 
inline double MTRand::rand()
	{ return double(randInt()) * (1.0/4294967295.0); }
 
inline double MTRand::rand( const double n )
	{ return rand() * n; }
 
inline double MTRand::randExc()
	{ return double(randInt()) * (1.0/4294967296.0); }
 
inline double MTRand::randExc( const double n )
	{ return randExc() * n; }
 
inline double MTRand::randDblExc()
	{ return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
 
inline double MTRand::randDblExc( const double n )
	{ return randDblExc() * n; }
 
inline double MTRand::rand53()
{
	uint32 a = randInt() >> 5, b = randInt() >> 6;
	return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
}
 
inline double MTRand::randNorm( const double mean, const double stddev )
{
	// Return a real number from a normal (Gaussian) distribution with given
	// mean and standard deviation by polar form of Box-Muller transformation
    /*
	double x, y, r;
	do
	{
		x = 2.0 * rand() - 1.0;
		y = 2.0 * rand() - 1.0;
		r = x * x + y * y;
	}
	while ( r >= 1.0 || r == 0.0 );
	double s = sqrt( -2.0 * log(r) / r );
	return mean + x * s * stddev;
        */
    if (!valid)
    {
        do {
            r1 = 2.0 * rand() - 1.0;
            r2 = 2.0 * rand() - 1.0;
            sqrs = r1 * r1 + r2 * r2;
        } while (sqrs >= 1.0 || sqrs == 0.0);
        rho = sqrt(-2.0 * log(sqrs)/sqrs);
        valid = true;
    } else {
        valid = false;
    }
    return rho * (valid ? r1 : r2) * stddev + mean;
}
 
inline double MTRand::operator()()
{
	return rand();
}
 
inline void MTRand::save( uint32* saveArray ) const
{
	register const uint32 *s = state;
	register uint32 *sa = saveArray;
	register int i = N;
	for( ; i--; *sa++ = *s++ ) {}
	*sa = left;
}
 
inline void MTRand::load( uint32 *const loadArray )
{
	register uint32 *s = state;
	register uint32 *la = loadArray;
	register int i = N;
	for( ; i--; *s++ = *la++ ) {}
	left = *la;
	pNext = &state[N-left];
}
 
inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
{
	register const MTRand::uint32 *s = mtrand.state;
	register int i = mtrand.N;
	for( ; i--; os << *s++ << "\t" ) {}
	return os << mtrand.left;
}
 
inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
{
	register MTRand::uint32 *s = mtrand.state;
	register int i = mtrand.N;
	for( ; i--; is >> *s++ ) {}
	is >> mtrand.left;
	mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
	return is;
}
 
inline MTRand& MTRand::operator=( const MTRand& o )
{
	if( this == &o ) return (*this);
	register const uint32 *t = o.state;
	register uint32 *s = state;
	register int i = N;
	for( ; i--; *s++ = *t++ ) {}
	left = o.left;
	pNext = &state[N-left];
        valid = o.valid;
        if (valid) {
            sqrs = o.sqrs;
            r1 = o.r1;
            r2 = o.r2;
            rho = o.rho;
        }
	return (*this);
}
 
#endif  // MERSENNETWISTER_H
 
// Change log:
//
// v0.1 - First release on 15 May 2000
//      - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
//      - Translated from C to C++
//      - Made completely ANSI compliant
//      - Designed convenient interface for initialization, seeding, and
//        obtaining numbers in default or user-defined ranges
//      - Added automatic seeding from /dev/urandom or time() and clock()
//      - Provided functions for saving and loading generator state
//
// v0.2 - Fixed bug which reloaded generator one step too late
//
// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
//
// v0.4 - Removed trailing newline in saved generator format to be consistent
//        with output format of built-in types
//
// v0.5 - Improved portability by replacing static const int's with enum's and
//        clarifying return values in seed(); suggested by Eric Heimburg
//      - Removed MAXINT constant; use 0xffffffffUL instead
//
// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
//      - Changed integer [0,n] generator to give better uniformity
//
// v0.7 - Fixed operator precedence ambiguity in reload()
//      - Added access for real numbers in (0,1) and (0,n)
//
// v0.8 - Included time.h header to properly support time_t and clock_t
//
// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
//      - Allowed for seeding with arrays of any length
//      - Added access for real numbers in [0,1) with 53-bit resolution
//      - Added access for real numbers from normal (Gaussian) distributions
//      - Increased overall speed by optimizing twist()
//      - Doubled speed of integer [0,n] generation
//      - Fixed out-of-range number generation on 64-bit machines
//      - Improved portability by substituting literal constants for long enum's
//      - Changed license from GNU LGPL to BSD
//
// v1.1 - Corrected parameter label in randNorm from "variance" to "stddev"
//      - Changed randNorm algorithm from basic to polar form for efficiency
//      - Updated includes from deprecated <xxxx.h> to standard <cxxxx> forms
//      - Cleaned declarations and definitions to please Intel compiler
//      - Revised twist() operator to work on ones'-complement machines
//      - Fixed reload() function to work when N and M are unsigned
//      - Added copy constructor and copy operator from Salvador Espana
Поделиться
Обновлено: Ноябрь 3, 2018 — 09:53

Добавить комментарий