------------------------制御データ----------------
誤差 0.1 出力 -2 出力ファイル kekka
順番 0 η 0.5 α 0.8 乱数 123
------------------------構造データ----------------
入力ユニット数 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
------------------------backpr.h------------------
/******************************************************/
/* バックプロパゲーションの制御（クラス BackControl） */
/******************************************************/
class BackControl {
protected:
double alpha, eata; //重み及びバイアス修正パラメータ
double eps; // 許容誤差
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のconstructor--
/*************************************/
/* クラスBackControlのコンストラクタ */
/* name : 入力データファイル名 */
/*************************************/
#include <stdio.h>
#include "backpr.h"
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 %*s %ld", &order, &eata, &alpha, &seed);
fclose(in);
}
------------------------BackDataのconstructor-----
/************************************/
/* クラスBackDataのコンストラクタ */
/* name : 入力データファイル名 */
/************************************/
#include <stdio.h>
#include "backpr.h"
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のdestructor------
/**************************/
/* BackDataのデストラクタ */
/**************************/
#include "backpr.h"
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のconstructor-------
/************************************************/
/* クラスBackprのコンストラクタ */
/* name_c : 制御データ用入力ファイル名 */
/* name_s : ネットワーク記述入力ファイル名 */
/************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "backpr.h"
#include "MT.h"
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;
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のdestructor--------
/************************/
/* Backprのデストラクタ */
/************************/
#include "backpr.h"
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;
}
------------------------Err_back.cpp--------------
/******************************************/
/* 誤差の計算，及び，重みとバイアスの修正 */
/* ptn[i1] : 出力パターン */
/******************************************/
#include "backpr.h"
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;
}
}
}
------------------------Forward.cpp---------------
/********************************************************/
/* 与えられた入力パターンに対する各ユニットの出力の計算 */
/********************************************************/
#include <math.h>
#include "backpr.h"
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));
}
}
------------------------Learn.cpp-----------------
/*****************************/
/* 学習の実行 */
/* p : 認識パターン */
/* m_tri : 最大学習回数 */
/*****************************/
#include <stdio.h>
#include <stdlib.h>
#include "backpr.h"
#include "MT.h"
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]));
}
}
------------------------Recog.cpp-----------------
/***********************************************/
/* 与えられた対象の認識と出力 */
/* p : 認識パターン */
/* pr : =0 : 出力を行わない */
/* =1 : 出力を行う */
/* =2 : 出力を行う（未学習パターン） */
/* tri : 現在の学習回数 */
/* return : 誤って認識したパターンの数 */
/***********************************************/
#include <stdio.h>
#include <math.h>
#include "backpr.h"
int Backpr::Recog(BackData &p, int pr, int tri)
{
int i1, i2, k1, No = 0, ln, sw;
FILE *out;
/*
ファイルのオープン
*/
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----------------------
/****************************/
/* back propagation model */
/* coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <stdlib.h>
#include "backpr.h"
/****************/
/* 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;
}
-----------------------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 */