/******************************/
/* 複雑な待ち行列問題 */
/* coded by Y.Suganuma */
/******************************/
#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <ctime>
#include <map>
#include <queue>
#include <vector>
#include <string>
#include "MT.h"
using namespace std;
/**********************/
/* 指数分布乱数の発生 */
/* m : 平均値 */
/* rerutn : 乱数 */
/**********************/
double Exp_b(double m)
{
return -m * log(genrand_real3());
}
/*********************/
/* クラスInletの定義 */
/*********************/
class Inlet {
string name; // 入り口名
string out; // 待ち行列名
int out_n; // 待ち行列番号
double arrive_time; // 客の到着時刻（負：客がない）
int a_t; // = -n : 到着する客の人数を負の値にしたもの
// =0 : 指数分布
// =1 : 一定時間間隔
double mean; // 到着時間間隔またはその平均
queue<double> que; // 客の到着時刻リスト
public:
/*************************************************************/
/* コンストラクタ */
/* name1 : 入り口名 */
/* a_t1; // = -n : 到着する客の人数を負の値にしたもの */
/* // =0 : 指数分布 */
/* // =1 : 一定時間間隔 */
/* m_a： 到着時間間隔またはその平均 */
/* que1 : 客の到着時刻リスト */
/* name2 : 待ち行列名 */
/*************************************************************/
Inlet (string name1, int a_t1, double m_a, queue<double> que1, string name2)
{
name = name1;
out = name2;
a_t = a_t1;
mean = m_a;
que = que1;
if (a_t == 0)
arrive_time = Exp_b(mean);
else if (a_t == 1)
arrive_time = 0;
else {
arrive_time = que.front();
que.pop();
}
}
friend class Q_base;
};
/*********************/
/* クラスQueueの定義 */
/*********************/
class Queue {
string name; // 待ち行列名
int nc; // 待ち行列への到着客数カウンタ
int max_q_l; // 最大待ち行列長
double c_ql; // 待ち行列長にその長さが継続した時間を乗じた値の累計
double ql_t; // 現在の待ち行列長になった時間
double max_wt; // 最大待ち時間
double c_wt; // 待ち時間の累計
int n; // =0 : 入り口から入る
// >0 : 複数の窓口から入る(窓口数)
vector<string> in; // 入り口名，または，窓口名
vector<int> in_n; // 入り口番号，または，窓口番号
int m; // 処理する窓口数
vector<string> out; // 窓口名
vector<int> out_n; // 窓口番号
queue<int> que; // 待ち行列
public:
/*********************************************/
/* コンストラクタ */
/* name1 : 待ち行列名 */
/* n1 : =0 : 入り口から入る */
/* >0 : 複数の窓口から入る(窓口数) */
/* in1 : 入り口名，または，窓口名 */
/* m1 : 処理する窓口数 */
/* out1 : 窓口名 */
/*********************************************/
Queue(string name1, int n1, vector<string> in1, int m1, vector<string> out1)
{
name = name1;
n = n1;
in = in1;
m = m1;
out = out1;
nc = 0;
max_q_l = 0;
c_ql = 0.0;
ql_t = 0.0;
max_wt = 0.0;
c_wt = 0.0;
}
friend class Q_base;
};
/**********************/
/* クラスEntityの定義 */
/**********************/
class Entity {
string name; // 窓口名
double end_time; // サービス終了時刻（負：何も処理していない）
int s_t; // =0 : 指数分布
// =1 : 一定時間
double mean; // 平均サービス時間
double busy_t; // 窓口がふさがった時刻
double c_busy; // 窓口がふさがっている時間の累計
int service; // サービス中の客番号（0のときは無し）
string in; // 待ち行列(入力)の名前
int in_n; // 待ち行列(入力)番号
int to; // =0 : システムの外に出る
// =1 : 待ち行列に入る
string out; // 待ち行列(出力)の名前
int out_n; // 待ち行列(出力)番号
public:
/****************************************/
/* コンストラクタ */
/* name1 : 窓口名 */
/* s_t1; // =0 : 指数分布 */
/* // =1 : 一定時間 */
/* m_s：サービス時間またはその平均 */
/* in : 待ち行列(入力)の名前 */
/* sw : =0 : システムの外に出る */
/* =1 : 待ち行列に入る */
/* out : 待ち行列(出力)の名前 */
/**********************************+++***/
Entity(string name1, int s_t1, double m_s, string name2, int sw, string name3)
{
name = name1;
to = sw;
in = name2;
if (to > 0)
out = name3;
end_time = -1.0;
s_t = s_t1;
mean = m_s;
service = 0;
busy_t = -1.0;
c_busy = 0.0;
}
friend class Q_base;
};
/************************/
/* クラスCustomerの定義 */
/************************/
class Customer
{
public :
double time; // 到着時刻
int state1; // 客の状態1
// =0 : 待ち行列に入っている
// =1 : サービスを受けている
int state2; // 客の状態2(待ち行列番号，または，窓口番号）
/*********************/
/* コンストラクタ */
/* s1,s2 : 状態 */
/* t : 到着時刻 */
/*******************************/
Customer::Customer(int s1, int s2, double t)
{
time = t;
state1 = s1;
state2 = s2;
}
};
/**********************/
/* クラスQ_baseの定義 */
/**********************/
class Q_base {
double p_time; // 現在時刻
int max_c; // 最大系内客数
int nc; // システムへの到着客数カウンタ
double now_c_t; // 現在の系内客数になった時間
double c_now_c; // 系内客数にその数が継続した時間を乗じた値の累計
double c_sys; // 滞在時間の累計
double max_sys; // 最大滞在時間
double end; // シミュレーション終了時間
map<int, Customer> cus; // 系内にいる客のリスト
vector<Inlet> inl; // Inletオブジェクトリスト
vector<Queue> que; // Queueオブジェクトリスト
vector<Entity> ent; // Entityオブジェクトリスト
public:
/***************************************/
/* コンストラクタ */
/* v_i : Inletオブジェクトリスト */
/* v_q : Queueオブジェクトリスト */
/* v_e : Entityオブジェクトリスト */
/* e : シミュレーション終了時刻 */
/***************************************/
Q_base(vector<Inlet> v_i, vector<Queue> v_q, vector<Entity> v_e, double e)
{
// 接続関係のチェック
cout << "\n";
// Inlet
inl = v_i;
que = v_q;
ent = v_e;
for (int i1 = 0; i1 < (int)inl.size()-1; i1++) {
for (int i2 = i1+1; i2 < (int)inl.size(); i2++) {
if (inl[i1].name == inl[i2].name) {
cout << "***error 同じ名前の入り口があります " + inl[i1].name + "\n";
exit(1);
}
}
}
int k;
for (int i1 = 0; i1 < (int)inl.size(); i1++) {
k = -1;
for (int i2 = 0; i2 < (int)que.size(); i2++) {
if (inl[i1].out == que[i2].name) {
k = i2;
break;
}
}
if (k >= 0)
inl[i1].out_n = k;
else {
cout << "***error 入り口から入る待ち行列名が不適当です " + inl[i1].out + "\n";
exit(1);
}
}
// Queue
for (int i1 = 0; i1 < (int)que.size()-1; i1++) {
for (int i2 = i1+1; i2 < (int)que.size(); i2++) {
if (que[i1].name == que[i2].name) {
cout << "***error 同じ名前の待ち行列があります " + que[i1].name + "\n";
exit(1);
}
}
}
for (int i1 = 0; i1 < (int)que.size(); i1++) {
if (que[i1].n == 0) {
k = -1;
for (int i2 = 0; i2 < (int)inl.size(); i2++) {
if (que[i1].in[0] == inl[i2].name) {
k = i2;
break;
}
}
if (k >= 0)
que[i1].in_n.push_back(k);
else {
cout << "***error 待ち行列に入る入り口名が不適当です " + que[i1].in[0] + "\n";
exit(1);
}
}
else {
for (int i2 = 0; i2 < (int)que[i1].n; i2++) {
k = -1;
for (int i3 = 0; i3 < (int)ent.size(); i3++) {
if (que[i1].in[i2] == ent[i3].name) {
k = i3;
break;
}
}
if (k >= 0)
que[i1].in_n.push_back(k);
else {
cout << "***error 待ち行列に入る窓口名が不適当です " + que[i1].in[i2] + "\n";
exit(1);
}
}
}
for (int i2 = 0; i2 < (int)que[i1].m; i2++) {
k = -1;
for (int i3 = 0; i3 < (int)ent.size(); i3++) {
if (que[i1].out[i2] == ent[i3].name) {
k = i3;
break;
}
}
if (k >= 0)
que[i1].out_n.push_back(k);
else {
cout << "***error 待ち行列を処理する窓口名が不適当です " + que[i1].out[i2] + "\n";
exit(1);
}
}
}
// Entity
for (int i1 = 0; i1 < (int)ent.size()-1; i1++) {
k = -1;
for (int i2 = i1+1; i2 < (int)ent.size(); i2++) {
if (ent[i1].name == ent[i2].name) {
cout << "***error 同じ名前の窓口があります " + ent[i1].name + "\n";
exit(1);
}
}
}
for (int i1 = 0; i1 < (int)ent.size(); i1++) {
k = -1;
for (int i2 = 0; i2 < (int)que.size(); i2++) {
if (ent[i1].in == que[i2].name) {
k = i2;
break;
}
}
if (k >= 0)
ent[i1].in_n = k;
else {
cout << "***error 窓口に入る待ち行列名が不適当です " + ent[i1].in + "\n";
exit(1);
}
if (ent[i1].to > 0) {
k = -1;
for (int i2 = 0; i2 < (int)que.size(); i2++) {
if (ent[i1].out == que[i2].name) {
k = i2;
break;
}
}
if (k >= 0)
ent[i1].out_n = k;
else {
cout << "***error 窓口の出口にある待ち行列名が不適当です " + ent[i1].out + "\n";
exit(1);
}
}
}
// 初期設定
p_time = 0.0;
max_c = 0;
nc = 0;
now_c_t = 0.0;
c_now_c = 0.0;
c_sys = 0.0;
max_sys = 0.0;
end = e;
// 乱数の初期設定
init_genrand((unsigned)time(NULL));
}
void Control(); // 全体の制御
void Next(int *); // 次の処理の決定
int End_o_s(); // 終了処理
void Arrive(int); // 客の到着
void Service(int); // サービス終了
void Output(); // 結果の出力
};
/**************/
/* 全体の制御 */
/**************/
void Q_base::Control()
{
int sw[2];
sw[0] = 0;
while (sw[0] > -2) {
Next(sw); // 次の処理の選択
if (sw[0] == -1)
sw[0] = End_o_s(); // シミュレーションの終了
else {
if (sw[0] == 0)
Arrive(sw[1]); // 客の到着処理
else
Service(sw[1]); // サービスの終了
}
}
}
/*********************************************/
/* 次の処理の決定 */
/* sw[0] : =-1 : シミュレーションの終了 */
/* =0 : 客の到着処理 */
/* =1 : 窓口のサービス終了 */
/* [1] : 入り口番号 or 窓口番号 */
/*********************************************/
void Q_base::Next(int *sw)
{
double tm = end; // 次の処理時刻
sw[0] = -1;
// 次の客の到着時刻
for (int i1 = 0; i1 < (int)inl.size(); i1++) {
Inlet x = inl[i1];
if (x.arrive_time >= 0.0 && x.arrive_time < tm) {
sw[0] = 0;
sw[1] = i1;
tm = x.arrive_time;
}
}
// サービス終了時刻
for (int i1 = 0; i1 < (int)ent.size(); i1++) {
Entity x = ent[i1];
if (x.service > 0 && x.end_time <= tm) {
sw[0] = 1;
sw[1] = i1;
tm = x.end_time;
}
}
if (sw[0] < 0)
end = p_time;
}
/**********************************/
/* 終了処理 */
/* return : =-1 : 終了前処理 */
/* =-2 : 実際の終了 */
/**********************************/
int Q_base::End_o_s()
{
int sw = -2;
p_time = end; // 現在時刻
for (int i1 = 0; i1 < (int)ent.size(); i1++) {
Entity x = ent[i1];
if (x.service > 0) { // サービス中の客はいるか？
if (sw == -2) {
sw = -1;
end = x.end_time; // 窓口i1のサービス終了時刻
}
else {
if (x.end_time > end)
end = x.end_time; // 窓口i1のサービス終了時刻
}
}
}
return sw;
}
/************************/
/* 客の到着処理 */
/* kk : 入り口番号 */
/************************/
void Q_base::Arrive(int kk)
{
/*
客数の増加と次の客の到着時刻の設定
*/
nc += 1; // 到着客数カウンタ
p_time = inl[kk].arrive_time; // 現在時刻
if (inl[kk].a_t == 0) // 次の客の到着時刻
inl[kk].arrive_time = p_time + Exp_b(inl[kk].mean);
else if (inl[kk].a_t == 1)
inl[kk].arrive_time = p_time + inl[kk].mean;
else {
if (inl[kk].que.empty())
inl[kk].arrive_time = -1.0;
else {
inl[kk].arrive_time = inl[kk].que.front();
inl[kk].que.pop();
}
}
if (inl[kk].arrive_time >= end)
inl[kk].arrive_time = -1.0;
/*
系内客数の処理
*/
c_now_c += cus.size() * (p_time - now_c_t); // 系内客数にその数が継続した時間を乗じた値の累計
now_c_t = p_time; // 現在の系内客数になった時刻
if ((int)cus.size()+1 > max_c)
max_c = cus.size() + 1; // 最大系内客数
/*
空いている窓口を探す
*/
int k1 = inl[kk].out_n;
que[k1].nc++;
int k = -1;
for (int i1 = 0; i1 < que[k1].m; i1++) {
int k2 = que[k1].out_n[i1]; // 処理窓口
if (ent[k2].service == 0) {
k = k2;
break;
}
}
/*
窓口に空きがない場合
*/
if (k < 0) {
Customer ct_p(0, k1, p_time);
cus.insert(pair<int, Customer>(nc, ct_p)); // 客の登録（系内客数）
que[k1].c_ql += que[k1].que.size() * (p_time - que[k1].ql_t); // 待ち行列長にその長さが継続した時間を乗じた値の累計
que[k1].que.push(nc); // 客の登録（待ち行列）
que[k1].ql_t = p_time; // 現在の待ち行列長になった時刻
if ((int)que[k1].que.size() > que[k1].max_q_l)
que[k1].max_q_l = que[k1].que.size(); // 最大待ち行列長
}
/*
すぐサービスをうけられる場合
*/
else {
Customer ct_p(1, k, p_time);
cus.insert(pair<int, Customer>(nc, ct_p)); // 客の登録（系内客数）
ent[k].service = nc; // サービスを受けている客番号
ent[k].busy_t = p_time; // 窓口がふさがった時刻
if (ent[k].s_t == 0) // 窓口のサービス終了時刻
ent[k].end_time = p_time + Exp_b(ent[k].mean);
else
ent[k].end_time = p_time + ent[k].mean;
}
}
/**********************************/
/* サービス終了時の処理 */
/* kk : サービス終了窓口番号 */
/**********************************/
void Q_base::Service(int kk)
{
map<int, Customer>::iterator it;
/*
時間の設定
*/
p_time = ent[kk].end_time; // 現在時刻
ent[kk].end_time = -1.0; // サービス終了時間
/*
システムの外へ出る場合
*/
if (ent[kk].to == 0) {
/*
系内客数の処理
*/
c_now_c += cus.size() * (p_time - now_c_t); // 系内客数にその数が継続した時間を乗じた値の累計
now_c_t = p_time; // 現在の系内客数になった時刻
/*
滞在時間の処理
*/
it = cus.find(ent[kk].service); // サービス中の客
double x1 = p_time - (it->second).time;
c_sys += x1; // 滞在時間の累計
if (x1 > max_sys)
max_sys = x1; // 最大滞在時間
cus.erase(ent[kk].service); // 客の削除（系内客数）
}
/*
他の窓口処理へ入る場合の処理
*/
else {
int k1 = ent[kk].out_n;
que[k1].nc++;
int sw = 1;
int k2 = 0;
if (que[k1].que.size() == 0) {
for (int i1 = 0; i1 < que[k1].m; i1++) {
k2 = que[k1].out_n[i1]; // 窓口
if (ent[k2].service == 0) {
sw = 0;
break;
}
}
}
/*
待ち行列がある場合
*/
if (sw > 0) {
que[k1].c_ql += que[k1].que.size() * (p_time - que[k1].ql_t); // 待ち行列長にその長さが継続した時間を乗じた値の累計
que[k1].que.push(ent[kk].service); // 客の登録（待ち行列）
que[k1].ql_t = p_time; // 現在の待ち行列長になった時刻
if ((int)que[k1].que.size() > que[k1].max_q_l)
que[k1].max_q_l = que[k1].que.size(); // 最大待ち行列長
it = cus.find(ent[kk].service);
(it->second).state1 = 0; // 客の状態変更(待ち行列)
(it->second).state2 = ent[kk].out_n; // 客の状態変更(待ち行列番号)
}
/*
すぐサービスをうけられる場合
*/
else {
ent[k2].service = ent[kk].service; // サービスを受けている客番号
ent[k2].busy_t = p_time; // 窓口がふさがった時刻
if (ent[k2].s_t == 0) // 窓口のサービス終了時刻
ent[k2].end_time = p_time + Exp_b(ent[k2].mean);
else
ent[k2].end_time = p_time + ent[k2].mean;
}
}
/*
窓口使用時間の処理
*/
ent[kk].service = 0; // 窓口を空き状態にする
ent[kk].c_busy += (p_time - ent[kk].busy_t); // 窓口がふさがっている時間の累計
/*
この窓口に対する待ち行列がある場合
*/
int k3 = ent[kk].in_n;
if (que[k3].que.size() > 0) {
que[k3].c_ql += que[k3].que.size() * (p_time - que[k3].ql_t); // 待ち行列長にその長さが継続した時間を乗じた値の累計
int n = que[k3].que.front(); // 待ち行列の先頭にいる客
que[k3].que.pop(); // 客の削除（待ち行列）
que[k3].ql_t = p_time; // 現在の待ち行列長になった時刻
it = cus.find(n);
double x1 = p_time - (it->second).time;
que[k3].c_wt += x1; // 待ち時間の累計
if (x1 > que[k3].max_wt)
que[k3].max_wt = x1; // 最大待ち時間
for (int i1 = 0; i1 < que[k3].m; i1++) {
int k4 = que[k3].out_n[i1]; // 窓口
if (ent[k4].service == 0) {
ent[k4].service = n; // 窓口の客番号
ent[k4].busy_t = p_time; // 窓口がふさがった時刻
if (ent[k4].s_t == 0) // 窓口のサービス終了時刻
ent[k4].end_time = p_time + Exp_b(ent[k4].mean);
else
ent[k4].end_time = p_time + ent[k4].mean;
(it->second).state1 = 1; // 客の状態変更(サービス中)
(it->second).state2 = k4; // 客の状態変更(窓口番号)
break;
}
}
}
}
/**************************/
/* 統計量の計算と最終出力 */
/**************************/
void Q_base::Output()
{
// System
printf("全客数 %d", nc);
printf(" 最大系内客数 %d 最大滞在時間 %.3f\n", max_c, max_sys);
printf("平均系内客数 %.3f", c_now_c / p_time);
printf(" 平均滞在時間 %.3f", c_sys / nc);
printf(" 終了時間 %.3f\n", p_time);
// Entity
for (int i1 = 0; i1 < (int)ent.size(); i1++) {
Entity e = ent[i1];
cout << "Entity " << e.name;
printf(" 稼働率 %.3f\n", e.c_busy / p_time);
}
// Queue
for (int i1 = 0; i1 < (int)que.size(); i1++) {
Queue q = que[i1];
cout << "Queue " << q.name;
printf(" 客数 %d", q.nc);
printf(" 最大待ち行列長 %d", q.max_q_l);
printf(" 最大待ち時間 %.3f\n", q.max_wt);
printf(" 平均待ち行列長 %.3f", q.c_ql / p_time);
printf(" 平均待ち時間 %.3f\n", q.c_wt / q.nc);
}
}
/****************/
/* main program */
/****************/
int main()
{
// 入り口
int n_i;
printf("入り口(Inlet)の数は？ ");
scanf("%d", &n_i);
vector<Inlet> inl;
for (int i1 = 0; i1 < n_i; i1++) {
double m_a;
string name1, name2;
printf("%d 番目の入り口(Inlet)\n", i1+1);
printf(" 名前は？ ");
cin >> name1;
int n;
queue<double> que;
printf(" 到着分布（=0:指数分布，=1:一定時間間隔，<0：指定，客数の負値）？ ");
scanf("%d", &n);
if (n == 0) {
printf(" 到着時間間隔の平均値は？ ");
scanf("%lf", &m_a);
}
else if (n == 1) {
printf(" 到着時間間隔は？ ");
scanf("%lf", &m_a);
}
else {
double x;
for (int i2 = 0; i2 < -n; i2++) {
printf(" 到着時間は？ ");
scanf("%lf", &x);
que.push(x);
}
}
printf(" 並ぶ待ち行列の名前は？ ");
cin >> name2;
Inlet inl_e(name1, n, m_a, que, name2);
inl.push_back(inl_e);
}
// 待ち行列
int n_q;
printf("待ち行列(Queue)の数は？ ");
scanf("%d", &n_q);
vector<Queue> que;
for (int i1 = 0; i1 < n_q; i1++) {
string name1;
printf("%d 番目の待ち行列(Queue)\n", i1+1);
printf(" 名前は？ ");
cin >> name1;
int n;
printf(" 入り口(0)，または，窓口(n>0，窓口の数)から？ ");
scanf("%d", &n);
vector<string> in;
if (n == 0) {
string name2;
printf(" 入り口の名前は？ ");
cin >> name2;
in.push_back(name2);
}
else {
for (int i2 = 0; i2 < n; i2++) {
string name3;
printf(" %d 番目の窓口の名前は？ ", i2+1);
cin >> name3;
in.push_back(name3);
}
}
int m;
printf(" 処理する窓口の数は？ ");
scanf("%d", &m);
vector<string> out;
for (int i2 = 0; i2 < m; i2++) {
string name4;
printf(" %d 番目の窓口の名前は？ ", i2+1);
cin >> name4;
out.push_back(name4);
}
Queue que_e(name1, n, in, m, out);
que.push_back(que_e);
}
// 窓口
int n_e;
printf("窓口(Entity)の数は？ ");
scanf("%d", &n_e);
vector<Entity> ent;
for (int i1 = 0; i1 < n_e; i1++) {
double m_s;
string name1;
printf("%d 番目の窓口(Entity)\n", i1+1);
printf(" 名前は？ ");
cin >> name1;
int s_t;
printf(" サービス分布（=0:指数分布，=1:一定時間）？ ");
scanf("%d", &s_t);
if (s_t == 0) {
printf(" サービス時間の平均値は？ ");
scanf("%lf", &m_s);
}
else {
printf(" サービス時間は？ ");
scanf("%lf", &m_s);
}
printf(" 待ち行列(入力)の名前は？ ");
string name2;
cin >> name2;
int sw;
printf(" 終了後，外部(0)，または，待ち行列(1)？ ");
scanf("%d", &sw);
string name3;
if (sw > 0) {
printf(" 待ち行列(出力)の名前は？ ");
cin >> name3;
}
Entity ent_e(name1, s_t, m_s, name2, sw, name3);
ent.push_back(ent_e);
}
// 全体の制御を行うクラス
double end;
printf("シミュレーション終了時間？ ");
scanf("%lf", &end);
Q_base base(inl, que, ent, end); // 全体の制御を行うクラス
// 実行
base.Control();
// 出力
base.Output();
return 0;
}
/*
------------入力例(簡単な場合)-----------
1
Inlet
0
5
Queue
1
Queue
0
Inlet
2
Entity1
Entity2
2
Entity1
0
4
Queue
0
Entity2
0
4
Queue
0
10000
------------入力例(複雑な場合)-----------
2
Inlet1
0
5
Queue1
Inlet2
0
5
Queue2
3
Queue1
0
Inlet1
1
Entity1
Queue2
0
Inlet2
1
Entity2
Queue3
2
Entity1
Entity2
2
Entity3
Entity4
4
Entity1
0
4
Queue1
1
Queue3
Entity2
0
4
Queue2
1
Queue3
Entity3
0
3
Queue3
0
Entity4
0
3
Queue3
0
10000
*/
-----------------------MT.h--------------------
/*
A C-program for MT19937, with initialization improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.
Before using, initialize the state by using init_genrand(seed)
or init_by_array(init_key, key_length).
Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
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.
Any feedback is very welcome.
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
*/
/*
The original version of http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c was modified by Takahiro Omi as
- delete line 47 "#include<stdio.h>"
- delete line 174 int main(void){...}
- change N -> MT_N
- change N -> MT_N
- change the file name "mt19937ar.c" -> "MT.h"
*/
/* Period parameters */
#define MT_N 624
#define MT_M 397
#define MATRIX_A 0x9908b0dfUL /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
static unsigned long mt[MT_N]; /* the array for the state vector */
static int mti=MT_N+1; /* mti==MT_N+1 means mt[MT_N] is not initialized */
/* initializes mt[MT_N] with a seed */
void init_genrand(unsigned long s)
{
mt[0]= s & 0xffffffffUL;
for (mti=1; mti<MT_N; mti++) {
mt[mti] =
(1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
/* In the previous versions, MSBs of the seed affect */
/* only MSBs of the array mt[]. */
/* 2002/01/09 modified by Makoto Matsumoto */
mt[mti] &= 0xffffffffUL;
/* for >32 bit machines */
}
}
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(unsigned long init_key[], int key_length)
{
int i, j, k;
init_genrand(19650218UL);
i=1; j=0;
k = (MT_N>key_length ? MT_N : key_length);
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ init_key[j] + j; /* non linear */
mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
i++; j++;
if (i>=MT_N) { mt[0] = mt[MT_N-1]; i=1; }
if (j>=key_length) j=0;
}
for (k=MT_N-1; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
- i; /* non linear */
mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
i++;
if (i>=MT_N) { mt[0] = mt[MT_N-1]; i=1; }
}
mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
}
/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void)
{
unsigned long y;
static unsigned long mag01[2]={0x0UL, MATRIX_A};
/* mag01[x] = x * MATRIX_A for x=0,1 */
if (mti >= MT_N) { /* generate N words at one time */
int kk;
if (mti == MT_N+1) /* if init_genrand() has not been called, */
init_genrand(5489UL); /* a default initial seed is used */
for (kk=0;kk<MT_N-MT_M;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
for (;kk<MT_N-1;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+(MT_M-MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
y = (mt[MT_N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[MT_N-1] = mt[MT_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
mti = 0;
}
y = mt[mti++];
/* Tempering */
y ^= (y >> 11);
y ^= (y << 7) & 0x9d2c5680UL;
y ^= (y << 15) & 0xefc60000UL;
y ^= (y >> 18);
return y;
}
/* generates a random number on [0,0x7fffffff]-interval */
long genrand_int31(void)
{
return (long)(genrand_int32()>>1);
}
/* generates a random number on [0,1]-real-interval */
double genrand_real1(void)
{
return genrand_int32()*(1.0/4294967295.0);
/* divided by 2^32-1 */
}
/* generates a random number on [0,1)-real-interval */
double genrand_real2(void)
{
return genrand_int32()*(1.0/4294967296.0);
/* divided by 2^32 */
}
/* generates a random number on (0,1)-real-interval */
double genrand_real3(void)
{
return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
/* divided by 2^32 */
}
/* generates a random number on [0,1) with 53-bit resolution*/
double genrand_res53(void)
{
unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
return(a*67108864.0+b)*(1.0/9007199254740992.0);
}
/* These real versions are due to Isaku Wada, 2002/01/09 added */