/****************************/
/* back propagation model */
/* coded by Y.Suganuma */
/****************************/
#include
#include
#include
#include
#include "MT.h"
/******************************************************/
/* バックプロパゲーションの制御（クラス BackControl） */
/******************************************************/
class BackControl {
protected:
double alpha, eata; //重み及びバイアス修正パラメータ
double eps; // 許容誤差
unsigned long seed; // 乱数の初期値;
int order; // 入力パターンの与え方（=0:順番，=1:ランダム）
int p_type; // 出力先・方法の指定
// =0 : 誤って認識した数だけ出力
// =1 : 認識結果を出力
// =2 : 認識結果と重みを出力
// （負の時は，認識結果と重みをファイルへも出力）
char o_file[100]; // 出力ファイル名
public:
BackControl(char *); // コンストラクタ
};
/*****************************************************/
/* バックプロパゲーションのデータ（クラス BackData） */
/*****************************************************/
class BackData {
int noiu; // 入力ユニットの数
int noou; // 出力ユニットの数
int noip; // 入力パターンの数
double **iptn; // iptn[i][j] : (i+1)番目の入力パターンの(j+1)番目の
// 入力ユニットの入力値
// i=0,noip-1 j=0,noiu-1
double **optn; // optn[i][j] : (i+1)番目の入力パターンに対する(j+1)
// 番目の出力ユニットの目標出力値
// i=0,noip-1 j=0,noou-1
public:
BackData(char *); // コンストラクタ
~BackData(); // デストラクタ
friend class Backpr;
};
/*******************************************/
/* バックプロパゲーション（クラス Backpr） */
/*******************************************/
class Backpr : BackControl {
char **con; // con[i][j] : 各ユニットに対するバイアスの与え方，及び，
// 接続方法
// [i][i] : ユニット(i+1)のバイアスの与え方
// =-3 : 入力値で固定
// =-2 : 入力値を初期値として，学習により変更
// =-1 : 乱数で初期値を設定し，学習により変更
// [i][j] : ユニット(i+1)と(j+1)の接続方法(j>i)
// =0 : 接続しない
// =1 : 接続する(重みの初期値を乱数で設定し，学習)
// =2 : 接続する(重みの初期値を入力で与え，学習)
// =3 : 接続する(重みを入力値に固定)
// i=0,nou-1 j=0,nou-1
int f_type; // シグモイド関数のタイプ
// 0 : [0,1]
// 1 : [-1,1]）
int noiu; // 入力ユニットの数
int noou; // 出力ユニットの数
int *nohu; // nohu[i] : レベル(i+1)の隠れ層のユニット数（隠れ層
// には下から順に番号が付けられ，出力層はレ
// ベル(nolvl+1)の隠れ層とも見做される）
// i=0,nolvl
int nolvl; // 隠れユニットの階層数
int nou; // 入力，出力，及び，隠れ層の各ユニットの数の和（各ユニ
// ットには最も上の出力ユニットから，隠れ層の各ユニット，
// 及び，入力ユニットに至る一連のユニット番号が付けられ
// る）
double *dp; // dp[i] : ユニット(i+1)の誤差 i=0,nou-1
double *op; // op[i] : ユニット(i+1)の出力 i=0,nou-1
double *theta; //theta[i] : ユニット(i+1)のバイアス i=0,nou
double **w; // w[i][j] : ユニット(i+1)から(j+1)への重み(j>i)
// w[j][i] : (i+1)から(j+1)への重みの前回修正量(j>i)
// w[i][i] : ユニット(i+1)のバイアスの前回修正量
// i=0,nou-1 j=0,nou-1
public:
Backpr(char *, char *); // コンストラクタ
~Backpr(); // デストラクタ
void Err_back(double *); // 重みの修正
void Forward(); // 誤差の計算
void Learn(BackData &, int); // 学習
int Recog(BackData &, int, int); // 認識と学習
};
/*************************************/
/* クラスBackControlのコンストラクタ */
/* name : 入力データファイル名 */
/*************************************/
BackControl::BackControl(char *name)
{
FILE *in;
in = fopen(name, "r");
fscanf(in, "%*s %lf %*s %d", &eps, &p_type);
if (p_type < 0)
fscanf(in, "%*s %s", o_file);
fscanf(in, "%*s %d %*s %lf %*s %lf", &order, &eata, &alpha);
fclose(in);
seed = (unsigned long)time(NULL);
}
/************************************/
/* クラスBackDataのコンストラクタ */
/* name : 入力データファイル名 */
/************************************/
BackData::BackData(char *name)
{
int i1, i2;
FILE *in;
in = fopen(name, "r");
/*
入力パターン数等
*/
fscanf(in, "%*s %d %*s %d %*s %d", &noip, &noiu, &noou);
/*
領域の確保
*/
iptn = new double * [noip];
for (i1 = 0; i1 < noip; i1++)
iptn[i1] = new double [noiu];
optn = new double * [noip];
for (i1 = 0; i1 < noip; i1++)
optn[i1] = new double [noou];
/*
入力パターン及び各入力パターンに対する出力パターンの入力
*/
for (i1 = 0; i1 < noip; i1++) {
fscanf(in, "%*s");
for (i2 = 0; i2 < noiu; i2++)
fscanf(in, "%lf", &(iptn[i1][i2]));
fscanf(in, "%*s");
for (i2 = 0; i2 < noou; i2++)
fscanf(in, "%lf", &(optn[i1][i2]));
}
fclose(in);
}
/**************************/
/* BackDataのデストラクタ */
/**************************/
BackData::~BackData()
{
int i1;
for (i1 = 0; i1 < noip; i1++)
delete [] iptn[i1];
delete [] iptn;
for (i1 = 0; i1 < noip; i1++)
delete [] optn[i1];
delete [] optn;
}
/************************************************/
/* クラスBackprのコンストラクタ */
/* name_c : 制御データ用入力ファイル名 */
/* name_s : ネットワーク記述入力ファイル名 */
/************************************************/
Backpr::Backpr(char *name_c, char *name_s) : BackControl(name_c)
{
double x1, x2;
int i0, i1, i2, id, k, k1, k2, k3, k4, l1, l2, n, sw;
FILE *in;
in = fopen(name_s, "r");
/*
乱数の初期化
*/
init_genrand(seed);
/*
入力ユニット，出力ユニットの数，関数タイプ
*/
fscanf(in, "%*s %d %*s %d %*s %d", &noiu, &noou, &f_type);
nou = noiu + noou; // 入力ユニットと出力ユニットの和
/*
隠れユニットの階層数と各階層のユニット数
*/
fscanf(in, "%*s %d", &nolvl);
nohu = new int [nolvl+1];
nohu[nolvl] = noou; // 出力ユニットの数
if (nolvl > 0) {
fscanf(in, "%*s");
for (i1 = 0; i1 < nolvl; i1++) {
fscanf(in, "%d", &(nohu[i1]));
nou += nohu[i1];
}
}
/*
領域の確保
*/
con = new char * [nou];
for (i1 = 0; i1 < nou; i1++) {
con[i1] = new char [nou];
for (i2 = 0; i2 < nou; i2++)
con[i1][i2] = (i1 == i2) ? -1 : 0;
}
w = new double * [nou];
for (i1 = 0; i1 < nou; i1++) {
w[i1] = new double [nou];
for (i2 = 0; i2 < nou; i2++)
w[i1][i2] = 0.0;
}
dp = new double [nou];
op = new double [nou];
theta = new double [nou];
for (i1 = 0; i1 < nou; i1++)
theta[i1] = 0.2 * genrand_real3() - 0.1;
/*
各ユニットのバイアスとユニット間の接続関係
*/
// バイアス
// バイアスデータの数
fscanf(in, "%*s %d", &n);
fscanf(in, "%*s %*s %*s");
if (n > 0) {
// バイアスデータの処理
for (i0 = 0; i0 < n; i0++) {
// ユニット番号
fscanf(in, "%d", &k1);
// 不適当なユニット番号のチェック
if (k1 < 1 || k1 > (nou-noiu)) {
printf("***error ユニット番号 %d が不適当\n", k1);
exit(1);
}
// バイアスの与え方
k1--;
fscanf(in, "%d", &id);
con[k1][k1] = id;
// バイアスの初期設定
switch (con[k1][k1]) {
case -1:
fscanf(in, "%lf %lf", &x1, &x2);
theta[k1] = (x2 - x1) * genrand_real3() + x1;
break;
case -2:
fscanf(in, "%lf", &(theta[k1]));
break;
case -3:
fscanf(in, "%lf", &(theta[k1]));
break;
default:
printf("***error バイアスの与え方が不適当\n");
exit(1);
}
}
}
// 接続方法
// 接続データの数
fscanf(in, "%*s %d", &n);
fscanf(in, "%*s %*s %*s");
if (n > 0) {
// 接続データの処理
for (i0 = 0; i0 < n; i0++) {
// 接続情報
fscanf(in, "%d %d %d %d", &k1, &k2, &k3, &k4);
// 不適切な接続のチェック
sw = 0;
if (k1 < 1 || k2 < 1 || k3 < 1 || k4 < 1)
sw = 1;
else {
if (k1 > nou || k2 > nou || k3 > nou || k4 > nou)
sw = 1;
else {
if (k1 > k2 || k3 > k4)
sw = 1;
else {
if (k4 >= k1)
sw = 1;
else {
l1 = -1;
k = 0;
for (i1 = nolvl; i1 >= 0 && l1 < 0; i1--) {
k += nohu[i1];
if (k1 <= k)
l1 = i1;
}
l2 = -1;
k = 0;
for (i1 = nolvl; i1 >= 0 && l2 < 0; i1--) {
k += nohu[i1];
if (k4 <= k)
l2 = i1;
}
if (l2 <= l1)
sw = 1;
}
}
}
}
if (sw > 0) {
printf("***error ユニット番号が不適当（%d %d %d %d）\n",
k1, k2, k3, k4);
exit(1);
}
// 重みの初期値の与え方
k1--;
k2--;
k3--;
k4--;
fscanf(in, "%d", &id);
if (id == 1)
fscanf(in, "%lf %lf", &x1, &x2);
else {
if (id > 1)
fscanf(in, "%lf", &x1);
else {
if (id != 0) {
printf("***error 接続方法が不適当\n");
exit(1);
}
}
}
// 重みの初期値の設定
for (i1 = k3; i1 <= k4; i1++) {
for (i2 = k1; i2 <= k2; i2++) {
con[i1][i2] = id;
switch (id) {
case 0:
w[i1][i2] = 0.0;
break;
case 1:
w[i1][i2] = (x2 - x1) * genrand_real3() + x1;
break;
case 2:
w[i1][i2] = x1;
break;
case 3:
w[i1][i2] = x1;
break;
}
}
}
}
}
fclose(in);
}
/************************/
/* Backprのデストラクタ */
/************************/
Backpr::~Backpr()
{
int i1;
for (i1 = 0; i1 < nou; i1++)
delete [] con[i1];
delete [] con;
for (i1 = 0; i1 < nou; i1++)
delete [] w[i1];
delete [] w;
delete [] nohu;
delete [] dp;
delete [] op;
delete [] theta;
}
/******************************************/
/* 誤差の計算，及び，重みとバイアスの修正 */
/* ptn[i1] : 出力パターン */
/******************************************/
void Backpr::Err_back(double *ptn)
{
double x1;
int i1, i2;
for (i1 = 0; i1 < nou-noiu; i1++) {
// 誤差の計算
if (i1 < noou) {
if (f_type == 0)
dp[i1] = (ptn[i1] - op[i1]) * op[i1] * (1.0 - op[i1]);
else
dp[i1] = 0.5 * (ptn[i1] - op[i1]) *
(op[i1] - 1.0) * (op[i1] + 1.0);
}
else {
x1 = 0.0;
for (i2 = 0; i2 < i1; i2++) {
if (con[i2][i1] > 0)
x1 += dp[i2] * w[i2][i1];
}
if (f_type == 0)
dp[i1] = op[i1] * (1.0 - op[i1]) * x1;
else
dp[i1] = 0.5 * (op[i1] - 1.0) * (op[i1] + 1.0) * x1;
}
// 重みの修正
for (i2 = i1+1; i2 < nou; i2++) {
if (con[i1][i2] == 1 || con[i1][i2] == 2) {
x1 = eata * dp[i1] * op[i2] + alpha * w[i2][i1];
w[i2][i1] = x1;
w[i1][i2] += x1;
}
}
// バイアスの修正
if (con[i1][i1] >= -2) {
x1 = eata * dp[i1] + alpha * w[i1][i1];
w[i1][i1] = x1;
theta[i1] += x1;
}
}
}
/********************************************************/
/* 与えられた入力パターンに対する各ユニットの出力の計算 */
/********************************************************/
void Backpr::Forward()
{
double sum;
int i1, i2;
for (i1 = nou-noiu-1; i1 >= 0; i1--) {
sum = -theta[i1];
for (i2 = i1+1; i2 < nou; i2++) {
if (con[i1][i2] > 0)
sum -= w[i1][i2] * op[i2];
}
op[i1] = (f_type == 0) ? 1.0 / (1.0 + exp(sum)) :
1.0 - 2.0 / (1.0 + exp(sum));
}
}
/*****************************/
/* 学習の実行 */
/* p : 認識パターン */
/* m_tri : 最大学習回数 */
/*****************************/
void Backpr::Learn(BackData &p, int m_tri)
{
int i1, i2, k0 = -1, k1;
/*
エラーチェック
*/
if (noiu != p.noiu || noou != p.noou) {
printf("***error 入力または出力ユニットの数が違います\n");
exit(1);
}
for (i1 = 0; i1 < m_tri; i1++) {
/*
パターンを与える順番の決定
*/
if (order == 0) { // 順番
k0++;
if (k0 >= p.noip)
k0 = 0;
}
else { // ランダム
k0 = (int)(genrand_real3() * p.noip);
if (k0 >= p.noip)
k0 = p.noip - 1;
}
/*
出力ユニットの結果を計算
*/
k1 = nou - noiu;
for (i2 = 0; i2 < noiu; i2++)
op[k1+i2] = p.iptn[k0][i2];
Forward();
/*
重みとバイアスの修正
*/
Err_back(&(p.optn[k0][0]));
}
}
/***********************************************/
/* 与えられた対象の認識と出力 */
/* p : 認識パターン */
/* pr : =0 : 出力を行わない */
/* =1 : 出力を行う */
/* =2 : 出力を行う（未学習パターン） */
/* tri : 現在の学習回数 */
/* return : 誤って認識したパターンの数 */
/***********************************************/
int Backpr::Recog(BackData &p, int pr, int tri)
{
int i1, i2, k1, No = 0, ln, sw;
FILE *out = NULL;
/*
ファイルのオープン
*/
if (p_type < 0 && pr > 0) {
if (pr == 1) {
out = fopen(o_file, "w");
fprintf(out, "***学習パターン***\n\n");
}
else {
out = fopen(o_file, "a");
fprintf(out, "\n***未学習パターン***\n\n");
}
}
/*
各パターンに対する出力
*/
for (i1 = 0; i1 < p.noip; i1++) {
// 入力パターンの設定
k1 = nou - noiu;
for (i2 = 0; i2 < noiu; i2++)
op[k1+i2] = p.iptn[i1][i2];
// 出力の計算
Forward();
// 結果の表示
if (p_type != 0 && pr > 0) {
printf("入力パターン%4d ", i1+1);
for (i2 = 0; i2 < noiu; i2++) {
printf("%5.2f", op[k1+i2]);
if (i2 == noiu-1)
printf("\n");
else {
if(((i2+1) % 10) == 0)
printf("\n ");
}
}
printf("\n 出力パターン(理想) ");
for (i2 = 0; i2 < noou; i2++) {
printf("%10.3f", p.optn[i1][i2]);
if (i2 == noou-1)
printf("\n");
else {
if(((i2+1) % 5) == 0)
printf("\n ");
}
}
}
sw = 0;
if (p_type != 0 && pr > 0)
printf(" (実際) ");
for (i2 = 0; i2 < noou; i2++) {
if (p_type != 0 && pr > 0) {
printf("%10.3f", op[i2]);
if (i2 == noou-1)
printf("\n");
else {
if(((i2+1) % 5) == 0)
printf("\n ");
}
}
if (fabs(op[i2]-p.optn[i1][i2]) > eps)
sw = 1;
}
if (sw > 0)
No++;
if (p_type < 0 && pr > 0) {
fprintf(out, "入力パターン%4d ", i1+1);
for (i2 = 0; i2 < noiu; i2++) {
fprintf(out, "%5.2f", op[k1+i2]);
if (i2 == noiu-1)
fprintf(out, "\n");
else {
if(((i2+1) % 10) == 0)
fprintf(out, "\n ");
}
}
fprintf(out, "\n 出力パターン(理想) ");
for (i2 = 0; i2 < noou; i2++) {
fprintf(out, "%10.3f", p.optn[i1][i2]);
if (i2 == noou-1)
fprintf(out, "\n");
else {
if(((i2+1) % 5) == 0)
fprintf(out, "\n ");
}
}
fprintf(out, " (実際) ");
for (i2 = 0; i2 < noou; i2++) {
fprintf(out, "%10.3f", op[i2]);
if (i2 == noou-1)
fprintf(out, "\n");
else {
if(((i2+1) % 5) == 0)
fprintf(out, "\n ");
}
}
}
if (p_type != 0 && pr > 0) {
getchar();
if (i1 == 0)
getchar();
}
}
/*
重みの出力
*/
if ((p_type < -1 || p_type > 1) && pr == 1) {
printf(" 重み\n");
for (i1 = 0; i1 < nou-noiu; i1++) {
printf(" to%4d from ", i1+1);
ln = -1;
for (i2 = 0; i2 < nou; i2++) {
if (con[i1][i2] > 0) {
if (ln <= 0) {
if (ln < 0)
ln = 0;
else
printf("\n ");
}
printf("%4d%11.3f", i2+1, w[i1][i2]);
ln += 1;
if (ln == 4)
ln = 0;
}
}
printf("\n");
}
printf("\n バイアス ");
ln = 0;
for (i1 = 0; i1 < nou-noiu; i1++) {
printf("%4d%11.3f", i1+1, theta[i1]);
ln += 1;
if (ln == 4 && i1 != nou-noiu-1) {
ln = 0;
printf("\n ");
}
}
if (ln != 0)
printf("\n");
getchar();
}
if (p_type < 0 && pr == 1) {
fprintf(out, " 重み\n");
for (i1 = 0; i1 < nou-noiu; i1++) {
fprintf(out, " to%4d from ", i1+1);
ln = -1;
for (i2 = 0; i2 < nou; i2++) {
if (con[i1][i2] > 0) {
if (ln <= 0) {
if (ln < 0)
ln = 0;
else
fprintf(out, "\n ");
}
fprintf(out, "%4d%11.3f", i2+1, w[i1][i2]);
ln += 1;
if (ln == 4)
ln = 0;
}
}
fprintf(out, "\n");
}
fprintf(out, "\n バイアス ");
ln = 0;
for (i1 = 0; i1 < nou-noiu; i1++) {
fprintf(out, "%4d%11.3f", i1+1, theta[i1]);
ln += 1;
if (ln == 4 && i1 != nou-noiu-1) {
ln = 0;
fprintf(out, "\n ");
}
}
if (ln != 0)
fprintf(out, "\n");
}
if (p_type < 0 && pr > 0)
fclose(out);
return No;
}
/****************/
/* main program */
/****************/
int main(int argc, char *argv[])
{
int conv; // 収束確認回数
int m_tri; // 最大学習回数
int max = 0, no = 1, sw, tri;
char f_name[100];
// エラー
if (argc != 3) {
printf("***error 入力データファイル名を指定して下さい\n");
exit(1);
}
else {
// ネットワークの定義
Backpr net(argv[1], argv[2]);
// 学習パターン等の入力
printf("学習回数は？ ");
scanf("%d", &m_tri);
printf("何回毎に収束を確認しますか？ ");
scanf("%d", &conv);
printf("学習パターンのファイル名は？ ");
scanf("%s", f_name);
BackData dt1(f_name);
// 学習
while (max < m_tri && no > 0) {
tri = ((max + conv) < m_tri) ? conv : m_tri - max;
max += tri;
net.Learn(dt1, tri); // 学習
no = net.Recog(dt1, 0, max); // 学習対象の認識
printf(" 回数 %d 誤って認識したパターン数 %d\n", max, no);
}
no = net.Recog(dt1, 1, max); // 学習対象の認識と出力
// 未学習パターンの認識
printf("未学習パターンの認識を行いますか？（=1:行う，=0:行わない） ");
scanf("%d", &sw);
if (sw > 0) {
printf("未学習パターンのファイル名は？ ");
scanf("%s", f_name);
BackData dt2(f_name);
no = net.Recog(dt2, 2, max); // 未学習対象の認識と出力
}
}
return 0;
}
/*
------------------------制御データ----------------
誤差 0.1 出力 -2 出力ファイル kekka
順番 0 η 0.5 α 0.8
------------------------構造データ----------------
入力ユニット数 2 出力ユニット数 1 関数タイプ 0
隠れ層の数 1 各隠れ層のユニット数（下から） 1
バイアス入力ユニット数 1
ユニット番号：出力ユニットから順に番号付け
入力方法：=-3:固定，=-2:入力後学習，=-1:乱数(default，[-0.1,0.1])）
値：バイアス値（ー２またはー３の時）または一様乱数の範囲（下限，上限）
1 -1 -0.05 0.05
接続方法の数 2
ユニット番号：ユニットk1からk2を，k3からk4に接続
接続方法：=0:接続なし，=1:乱数，=2:重み入力後学習，=3:重み固定
値：重み（２または３の時）または一様乱数の範囲（１の時：下限，上限）
3 4 1 2 1 -0.1 0.1
2 2 1 1 1 -0.1 0.1
------------------------学習データ----------------
パターンの数 4 入力ユニット数 2 出力ユニット数 1
入力１ 0 0
出力１ 0
入力２ 0 1
出力２ 1
入力３ 1 0
出力３ 1
入力４ 1 1
出力４ 0
------------------------認識データ----------------
パターンの数 4 入力ユニット数 2 出力ユニット数 1
入力１ 0 0
出力１ 0
入力２ 0 1
出力２ 1
入力３ 1 0
出力３ 1
入力４ 1 1
出力４ 0
//---------------------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"
// - 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> 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> 1) ^ mag01[y & 0x1UL];
}
for (;kk> 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
*/
