巡回セールスマン問題(分割法)

データ例,プログラム,MT.h

  以下,複数のファイル構成になっています.ファイル間の区切りを「---・・・」で示します.このプログラム例においては,ヘッダファイル MT.h に記述されたメルセンヌ・ツイスタを使用していまが,C++11 で記述可能であれば,標準ライブラリ内に含まれているメルセンヌ・ツイスタ法を使用した乱数生成関数を利用できます.

------------------------ケーススタディデータ------
問題の数 2
問題 data\pr439.tsp 乱数の数 10
 乱数 1 
 乱数 12
 乱数 123
 乱数 1234
 乱数 12345
 乱数 54321
 乱数 5432
 乱数 543
 乱数 54
 乱数 5
問題 data\kroA100.tsp 乱数の数 10
 乱数 1 
 乱数 12
 乱数 123
 乱数 1234
 乱数 12345
 乱数 54321
 乱数 5432
 乱数 543
 乱数 54
 乱数 5

------------------------データファイル------------
都市の数 439 選択方法(0:最良,1:最初) 1 近傍(2or3) 3 整数 1
出力(0:ディスプレイ,1:ファイル) -1 出力ファイル名 result\pr439
分割数 X 4 Y 3 最大試行回数 1000
7125 11300
7225 11050
         ・・・・・
1775 5375
2075 6475

------------------------partition.h---------------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "MT.h"

float kyori(int, int *, float **);

/*************************/
/* クラスPartitionの定義 */
/*************************/
class Partition {
		float **city;   //都市の位置データ
		float **city_i;   //都市の位置データ(作業領域)
		float *p_x;   // x軸の分割点
		float *p_y;   // y軸の分割点
		float **rg;   // 都市間の距離
		long seed;   // 乱数の初期値
		int fix;   // =1 : 近傍を固定
                   // =0 : 近傍を可変
		int max_try;   // 最大試行回数
		int n_city;   // 都市の数
		int **n_seq;   // 各領域の都市数
		int **n_seq1;   // 各領域の都市数(ワーク)
		int n_p_x;   // x軸方向の分割数
		int n_p_y;   // y軸方向の分割数
		int ***seq;   // 経路
		int ***seq1;   // 経路(ワーク)
		int *seq_w1;   // 作業領域
		int *seq_w2;   // 作業領域
		int neib;   // 近傍(2 or 3)
		int seisu;   // 位置データの表現方法
                     //   =1 : 整数
                     //   =-1 : 実数(距離を整数計算)
                     //   =-2 : 実数(距離を実数計算)
		int sel;   // エッジの選択方法
                   //   =0 : 最良のものを選択
                   //   =1 : 最初のものを選択
		int **state;   // 領域結合用ワーク
		char *i_file;   // 入力ファイル名
	public:
		int Max;   // 最適経路の長さ
		int out_m;   // 出力方法
                     //   =-1 : ディスプレイ(経路長だけ)
                     //   =0 : ディスプレイ
                     //   =1 : ファイル
                     //   =2 : ファイル(経路長だけ)
		char o_file[100];   // 出力ファイル名
		Partition(char *);   // コンストラクタ
		Partition();   // コンストラクタ
		~Partition();   // デストラクタ
		void Optimize(long);   // 最適化の実行
		void Output(int, int);   // 出力
		int Connect();   // 分割したものを一つにまとめる
};

------------------------iteration.h---------------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "MT.h"

float kyori(int, int *, float **);

/*************************/
/* クラスIterationの定義 */
/*************************/
class Iteration {
		float **city;   //都市の位置データ
		float **rg;   // 都市間の距離
		int fix;   // =1 : 近傍を固定
                   // =0 : 近傍を可変
		int max_try;   // 最大試行回数
		int n_city;   // 都市の数
		int out_d;   // 表示間隔
		int *seq_w1;   // 都市を訪れる順序(ワーク)
		int *seq_w2;   // 都市を訪れる順序(ワーク)
		int *seq_w3;   // 都市を訪れる順序(ワーク)
		int *seq_w4;   // 都市を訪れる順序(ワーク)
		int *seq_w5;   // 都市を訪れる順序(ワーク)
		int neib;   // 近傍(2 or 3)
		int out_lvl;   // 出力レベル
                       //   =0 : 最終出力だけ
                       //   n>0 : n世代毎に出力(負の時はファイル)
		int out_m;   // 出力方法
                     //   =-1 : 出力しない
                     //   =0 : すべてを出力
                     //   =1 : 評価値だけを出力(最終結果だけはすべてを出力)
		int seisu;   // 位置データの表現方法
                     //   =1 : 整数
                     //   =-1 : 実数(距離を整数計算)
                     //   =-2 : 実数(距離を実数計算)
		int sel;   // エッジの選択方法
                   //   =0 : 最良のものを選択
                   //   =1 : 最初のものを選択
		char o_file[100];   // 出力ファイル名
	public:
		int *seq;   // 都市を訪れる順序
		Iteration (int, int, int, int, int, int, int, int, int, char *, float **);// コンストラクタ
		~Iteration();   // デストラクタ
		int Optimize();   // 最適化の実行
		int Change(double *);   // 改善
		void Output(int, int, double);   // 出力
};

------------------------メンバー関数(Partition)---
/*********************************/
/* クラスPartitionのメンバー関数 */
/*********************************/
#include "partition.h"
#include "iteration.h"

/**************************/
/* コンストラクタ         */
/*      name : ファイル名 */
/**************************/
Partition::Partition(char *name)
{
	double x, y;
	float max_x, max_y, min_x, min_y, s_x, s_y;
	int i1, i2, i3, max = 0, n;
	FILE *in;
					// ファイルのオープン
	i_file = name;
	in     = fopen(name, "r");
	if (in == NULL) {
		printf("***error  データファイル名が不適当\n");
		exit(1);
	}
					// 基本データ
	fscanf(in, "%*s %d %*s %d %*s %d %*s %d", &n_city, &sel, &neib, &seisu);
	fscanf(in, "%*s %d %*s %s", &out_m, o_file);
	fscanf(in, "%*s %*s %d %*s %d %*s %d", &n_p_x, &n_p_y, &max_try);

	if (neib < 0) {
		neib = -neib;
		fix  = 0;
	}
	else
		fix = 1;
					// 都市の位置データ
	city = new float * [n_city];
	for (i1 = 0; i1 < n_city; i1++) {
		city[i1] = new float [2];
		fscanf(in, "%f %f", &city[i1][0], &city[i1][1]);
	}
					// ファイルのクローズ
	fclose(in);
					// 距離テーブルの作成
	rg = new float * [n_city];

	for (i1 = 0; i1 < n_city; i1++) {
		rg[i1] = new float [n_city];
		for (i2 = i1+1; i2 < n_city; i2++) {
			x          = city[i2][0] - city[i1][0];
			y          = city[i2][1] - city[i1][1];
			rg[i1][i2] = (float)sqrt(x * x + y * y);
			if (seisu > -2)
				rg[i1][i2] = (int)(rg[i1][i2] + 0.5);
		}
	}

	for (i1 = 1; i1 < n_city; i1++) {
		for (i2 = 0; i2 < i1; i2++)
			rg[i1][i2] = rg[i2][i1];
	}
					// 作業領域
	state  = new int * [n_p_y];
	n_seq  = new int * [n_p_y];
	n_seq1 = new int * [n_p_y];
	for (i1 = 0; i1 < n_p_y; i1++) {
		n_seq[i1]  = new int [n_p_x];
		n_seq1[i1] = new int [n_p_x];
		state[i1]  = new int [n_p_x];
	}

	seq  = new int ** [n_p_y];
	seq1 = new int ** [n_p_y];
	for (i1 = 0; i1 < n_p_y; i1++) {
		seq[i1]  = new int * [n_p_x];
		seq1[i1] = new int * [n_p_x];
	}

	seq_w1 = new int [n_city];
	seq_w2 = new int [n_city];
	p_x    = new float [n_p_x];
	p_y    = new float [n_p_y];
					// 都市の分割
	for (i1 = 0; i1 < n_city; i1++)
		seq_w1[i1] = 0;

	min_x = city[0][0];
	max_x = city[0][0];
	min_y = city[0][1];
	max_y = city[0][1];

	for (i1 = 1; i1 < n_city; i1++) {
		if (city[i1][0] < min_x)
			min_x = city[i1][0];
		else {
			if (city[i1][0] > max_x)
				max_x = city[i1][0];
		}
		if (city[i1][1] < min_y)
			min_y = city[i1][1];
		else {
			if (city[i1][1] > max_y)
				max_y = city[i1][1];
		}
	}

	s_x          = (max_x - min_x) / n_p_x;
	p_x[0]       = min_x + s_x;
	p_x[n_p_x-1] = max_x;
	for (i1 = 1; i1 < n_p_x-1; i1++)
		p_x[i1] = p_x[0] + i1 * s_x;

	s_y = (max_y - min_y) / n_p_y;
	p_y[0]       = min_y + s_y;
	p_y[n_p_y-1] = max_y;
	for (i1 = 1; i1 < n_p_y-1; i1++)
		p_y[i1] = p_y[0] + i1 * s_y;

	for (i1 = 0; i1 < n_p_y; i1++) {
		for (i2 = 0; i2 < n_p_x; i2++) {
			n = 0;
			for (i3 = 0; i3 < n_city; i3++) {
				if (seq_w1[i3] == 0) {
					if (city[i3][0] <= p_x[i2] && city[i3][1] <= p_y[i1]) {
						seq_w1[i3] = 1;
						seq_w2[n]  = i3;
						n++;
					}
				}
			}
			n_seq1[i1][i2] = n;
			if (n > 0) {
				seq[i1][i2]  = new int [n_city];
				seq1[i1][i2] = new int [n_city];
				for (i3 = 0; i3 < n; i3++)
					seq1[i1][i2][i3] = seq_w2[i3];
				if (n > max)
					max = n;
			}
		}
	}
					// 作業領域
	printf("max %d\n", max);
	city_i = new float * [max];
	for (i1 = 0; i1 < max; i1++)
		city_i[i1] = new float [2];
}

/******************/
/* コンストラクタ */
/******************/
Partition::Partition()
{
	n_city = 0;
}

/****************/
/* デストラクタ */
/****************/
Partition::~Partition()
{
	int i1, i2;

	if (n_city > 0) {

		for (i1 = 0; i1 < n_city; i1++) {
			delete [] rg[i1];
			delete [] city[i1];
			delete [] city_i[i1];
		}
		delete [] rg;
		delete [] city;
		delete [] city_i;

		for (i1 = 0; i1 < n_p_y; i1++)
			delete [] state[i1];
		delete [] state;

		delete [] seq_w1;
		delete [] seq_w2;
		delete [] p_x;
		delete [] p_y;

		for (i1 = 0; i1 < n_p_y; i1++) {
			for (i2 = 0; i2 < n_p_x; i2++) {
				if (n_seq1[i1][i2] > 0) {
					delete [] seq[i1][i2];
					delete [] seq1[i1][i2];
				}
			}
			delete [] seq[i1];
			delete [] seq1[i1];
		}
		delete [] seq;
		delete [] seq1;

		for (i1 = 0; i1 < n_p_y; i1++) {
			delete [] n_seq[i1];
			delete [] n_seq1[i1];
		}
		delete [] n_seq;
		delete [] n_seq1;
	}
}

/******************************/
/* 最適化の実行               */
/*      seed_i : 乱数の初期値 */
/******************************/
void Partition::Optimize(long seed_i)
{
	int i1, i2, i3, k, max, nb, r = 0;
	Iteration *it;
					// 初期設定
	seed = seed_i;
	init_genrand(seed);
					// 分割数と開始時間の出力
	if (out_m > 0)
		Output(0, r);

	for (i1 = 0; i1 < n_p_y; i1++) {
		for (i2 = 0; i2 < n_p_x; i2++) {
			n_seq[i1][i2] = n_seq1[i1][i2];
			for (i3 = 0; i3 < n_seq1[i1][i2]; i3++)
				seq[i1][i2][i3] = seq1[i1][i2][i3];
		}
	}
					// 分割毎の最適化
	for (i1 = 0; i1 < n_p_y; i1++) {
		for (i2 = 0; i2 < n_p_x; i2++) {
			if (n_seq[i1][i2] > 3) {
						// 近傍の大きさ
				nb = (n_seq[i1][i2] > 3) ? neib : 2;
						// 都市位置データの設定
				for (i3 = 0; i3 < n_seq[i1][i2]; i3++) {
					k             = seq[i1][i2][i3];
					city_i[i3][0] = city[k][0];
					city_i[i3][1] = city[k][1];
				}
						// 最適化
				it  = new Iteration (n_seq[i1][i2], max_try, seisu, sel, nb, fix,
                                     0, -1, 0, o_file, city_i);
				max = it->Optimize();
						// 結果の保存
				for (i3 = 0; i3 < n_seq[i1][i2]; i3++) {
					k          = it->seq[i3];
					seq_w1[i3] = seq[i1][i2][k];
				}
				for (i3 = 0; i3 < n_seq[i1][i2]; i3++)
					seq[i1][i2][i3] = seq_w1[i3];
						// 出力
				r = (seisu > -2) ? (int)kyori(n_seq[i1][i2], seq[i1][i2], rg) :
                                   (int)(kyori(n_seq[i1][i2], seq[i1][i2], rg) + 0.5);
				printf("   y %d x %d n_city %d range %d (trial %d)\n",
                       i1+1, i2+1, n_seq[i1][i2], r, max);
			}
		}
	}
					// 経路の接続
	r = Connect();
					// 出力
	Output(n_city, r);
}

/***********************/
/* 出力                */
/*      n_c : 都市の数 */
/*      r : 距離       */
/***********************/
void Partition::Output(int n_c, int r)
{
	int i1, k = 0, n;
	char *now;
	time_t aclock;
	FILE *out;

	if (out_m <= 0) {
		out = stdout;
		fprintf(out, "距離 %d\n", r);
		getchar();
	}
	else {
		time(&aclock);
		now = ctime(&aclock);
		out = fopen(o_file, "a");
		if (n_c > 0) {
			printf("距離 %d\n", r);
			fprintf(out, "   距離 %d 時間 %s\n", r, now);
		}
		else
			fprintf(out, "問題 %s 乱数 %ld 分割 %d %d 時間 %s",
                    i_file, seed, n_p_x, n_p_y, now);
	}

	if (n_c > 0 && (out_m == 0 || out_m == 1)) {
		for (i1 = 0; i1 < n_c; i1++) {
			n = seq_w1[i1];
			if (seisu > 0)
				fprintf(out, "  %d %d %d\n", n, (int)city[n][0], (int)city[n][1]);
			else
				fprintf(out, "  %d %f %f\n", n, city[n][0], city[n][1]);
			if (out_m == 0) {
				k++;
				if (k == 10) {
					getchar();
					k = 0;
				}
			}
		}
	}

	if (out_m > 0)
		fclose(out);
}

/************************/
/* 分割された領域の接続 */
/************************/
int Partition::Connect()
{
	double wd, wd1, wa1, wa2, min = 0;
	int i1, i2, i3, i4, k, k1 = 0, k2 = 0, k3 = 0, k4 = 0, min_c = 0, n, r,
        r1 = 0, r2 = 0, r3 = 0, r4 = 0, s1 = 0, s2 = 0, sw = 1;
/*
     領域が1つの場合
*/
	if (n_p_x == 1 && n_p_y == 1) {
		for (i1 = 0; i1 < n_seq[0][0]; i1++)
			seq_w1[i1] = seq[0][0][i1];
	}
/*
     初期設定
*/
	else {

		for (i1 = 0; i1 < n_p_y; i1++) {
			for (i2 = 0; i2 < n_p_x; i2++)
				state[i1][i2] = (n_seq[i1][i2] > 0) ? 0 : 1;
		}
/*
     実行
*/
		while (sw > 0) {
					// 最小節点領域
			min_c = n_city;
			sw    = 0;
			for (i1 = 0; i1 < n_p_y; i1++) {
				for (i2 = 0; i2 < n_p_x; i2++) {
					if (state[i1][i2] == 0 && n_seq[i1][i2] < min_c) {
						sw    = 1;
						r1    = i1;
						r2    = i2;
						min_c = n_seq[i1][i2];
					}
				}
			}
					// 結合する対象領域の決定
			if (sw > 0) {
				sw = 0;
				for (i1 = 0; i1 < n_p_y; i1++) {
					for (i2 = 0; i2 < n_p_x; i2++) {
						if (state[i1][i2] == 0 && (i1 != r1 || i2 !=r2)) {
								// 節点の数>2
							if (n_seq[r1][r2] > 1) {
								for (i3 = 0; i3 < n_seq[r1][r2]; i3++) {
									k1  = seq[r1][r2][i3];
									k2  = (i3 == n_seq[r1][r2]-1) ? seq[r1][r2][0] :
                                                                    seq[r1][r2][i3+1];
									wd1 = rg[k1][k2];
									for (i4 = 0; i4 < n_seq[i1][i2]; i4++) {
										k3  = seq[i1][i2][i4];
										k4  = (i4 == n_seq[i1][i2]-1) ? seq[i1][i2][0] :
                                                                        seq[i1][i2][i4+1];
										wd  = wd1 + rg[k3][k4];
										wa1 = rg[k1][k3] + rg[k2][k4];
										wa2 = rg[k1][k4] + rg[k2][k3];
										if (sw == 0 || wa1-wd < min) {
											min = wa1 - wd;
											r3  = i1;
											r4  = i2;
											s1  = (i3 == n_seq[r1][r2]-1) ? 0 : i3 + 1;
											s2  = (i4 == n_seq[i1][i2]-1) ? 0 : i4 + 1;
											sw  = -1;
										}
										if (sw == 0 || wa2-wd < min) {
											min = wa2 - wd;
											r3  = i1;
											r4  = i2;
											s1  = i3;
											s2  = (i4 == n_seq[i1][i2]-1) ? 0 : i4 + 1;
											sw  = 1;
										}
									}
								}
							}
								// 節点の数=1
							else {
								k1 = seq[r1][r2][0];
								if (n_seq[i1][i2] > 1) {
									for (i4 = 0; i4 < n_seq[i1][i2]; i4++) {
										k3  = seq[i1][i2][i4];
										k4  = (i4 == n_seq[i1][i2]-1) ? seq[i1][i2][0] :
                                                                        seq[i1][i2][i4+1];
										wd  = rg[k3][k4];
										wa1 = rg[k1][k3] + rg[k1][k4];
										if (sw == 0 || wa1-wd < min) {
											min = wa1 - wd;
											r3  = i1;
											r4  = i2;
											s1  = 0;
											s2  = (i4 == n_seq[i1][i2]-1) ? 0 : i4 + 1;
											sw  = 1;
										}
									}
								}
								else {
									k3  = seq[i1][i2][0];
									wa1 = rg[k1][k3];
									if (sw == 0 || wa1 < min) {
										min = wa1;
										r3  = i1;
										r4  = i2;
										s1  = 0;
										s2  = 0;
										sw  = 1;
									}
								}
							}
						}
					}
				}
					// 領域の結合
				seq_w1[0] = seq[r1][r2][s1];
				k         = 1;
				n         = s2;
				for (i1 = 0; i1 < n_seq[r3][r4]; i1++) {
					seq_w1[k] = seq[r3][r4][n];
					k++;
					n++;
					if (n > n_seq[r3][r4]-1)
						n = 0;
				}
				if (sw > 0) {
					n = s1 + 1;
					for (i1 = 0; i1 < n_seq[r1][r2]-1; i1++) {
						if (n > n_seq[r1][r2]-1)
							n = 0;
						seq_w1[k] = seq[r1][r2][n];
						k++;
						n++;
					}
				}
				else {
					n = s1 - 1;
					for (i1 = 0; i1 < n_seq[r1][r2]-1; i1++) {
						if (n < 0)
							n = n_seq[r1][r2] - 1;
						seq_w1[k] = seq[r1][r2][n];
						k++;
						n--;
					}
				}
					// 状態の変更
				n_seq[r1][r2] += n_seq[r3][r4];
				state[r3][r4]  = 1;
				for (i1 = 0; i1 < n_seq[r1][r2]; i1++)
					seq[r1][r2][i1] = seq_w1[i1];
				sw  = 1;
			}
		}
	}

	r   = (seisu > -2) ? (int)kyori(n_city, seq_w1, rg) :
                         (int)(kyori(n_city, seq_w1, rg) + 0.5);
	Max = r;

	return r;
}

------------------------メンバー関数(Iteration)---
/*********************************/
/* クラスIterationのメンバー関数 */
/*********************************/
#include "iteration.h"

/**********************************/
/* コンストラクタ                 */
/*      n_city_i : 都市の数       */
/*      max_try_i : 最大試行回数  */
/*      sei_i : 整数 or 実数      */
/*      sel_i : エッジの選択方法  */
/*      neib_i : 近傍             */
/*      fix_i : 近傍の扱い方      */
/*      out_lvl_i : 出力レベル    */
/*      out_m_i : 出力方法        */
/*      out_d_i : 表示間隔        */
/*      o_file_i : 出力ファイル名 */
/*      city_i : 都市の位置データ */
/**********************************/
Iteration::Iteration (int n_city_i, int max_tri_i, int sei_i, int sel_i, int neib_i, int fix_i, 
                      int out_lvl_i, int out_m_i, int out_d_i, char *o_file_i, float **city_i)
{
	double x, y;
	int ct, i1, i2, sw;
					// 値の設定
	n_city  = n_city_i;
	max_try = max_tri_i;
	seisu   = sei_i;
	sel     = sel_i;
	neib    = neib_i;
	fix     = fix_i;
	out_lvl = out_lvl_i;
	out_m   = out_m_i;
	out_d   = out_d_i;

	strcpy(o_file, o_file_i);
					// 都市の位置データ
	city = new float * [n_city];
	for (i1 = 0; i1 < n_city; i1++) {
		city[i1]    = new float [2];
		city[i1][0] = city_i[i1][0];
		city[i1][1] = city_i[i1][1];
	}
					// 距離テーブルの作成
	rg = new float * [n_city];

	for (i1 = 0; i1 < n_city; i1++) {
		rg[i1] = new float [n_city];
		for (i2 = i1+1; i2 < n_city; i2++) {
			x          = city[i2][0] - city[i1][0];
			y          = city[i2][1] - city[i1][1];
			rg[i1][i2] = (float)sqrt(x * x + y * y);
			if (seisu > -2)
				rg[i1][i2] = (int)(rg[i1][i2] + 0.5);
		}
	}

	for (i1 = 1; i1 < n_city; i1++) {
		for (i2 = 0; i2 < i1; i2++)
			rg[i1][i2] = rg[i2][i1];
	}
					// 都市を訪れる順序(初期設定)
	seq    = new int [n_city];
	seq_w1 = new int [n_city];
	seq_w2 = new int [n_city];
	seq_w3 = new int [n_city];
	seq_w4 = new int [n_city];
	seq_w5 = new int [n_city];

	for (i1 = 0; i1 < n_city; i1++) {
		sw = 0;
		while (sw == 0) {
			ct = (int)(genrand_real3() * n_city);
			if (ct >= n_city)
				ct = n_city - 1;
			seq[i1] = ct;
			sw      = 1;
			for (i2 = 0; i2 < i1 && sw > 0; i2++) {
				if (ct == seq[i2])
					sw = 0;
			}
		}
	}
}

/****************/
/* デストラクタ */
/****************/
Iteration::~Iteration ()
{
	int i1;

	if (n_city > 0) {

		for (i1 = 0; i1 < n_city; i1++) {
			delete [] city[i1];
			delete [] rg[i1];
		}
		delete [] city;
		delete [] rg;

		delete[] seq;
		delete [] seq_w1;
		delete [] seq_w2;
		delete [] seq_w3;
		delete [] seq_w4;
		delete [] seq_w5;
	}
}

/****************/
/* 最適化の実行 */
/****************/
int Iteration::Optimize ()
{
	double max;
	int n_tri, sw;
					// 初期設定
	n_tri = 0;
	max   = kyori(n_city, seq, rg);

	if (out_m >= 0 && abs(out_lvl) > 0) {
		if (seisu > -2)
			printf("***試行回数 %d 距離 %d\n", n_tri, (int)max);
		else
			printf("***試行回数 %d 距離 %f\n", n_tri, max);
		Output(out_lvl, n_tri, max);
	}
					// 実行
	sw = 1;
	for (n_tri = 1; n_tri <= max_try && sw > 0; n_tri++) {
						// 改善
		sw = Change(&max);
						// 出力
		if (out_d > 0 && n_tri%out_d == 0) {
			if (seisu > -2)
				printf("***試行回数 %d 距離 %d\n", n_tri, (int)max);
			else
				printf("***試行回数 %d 距離 %f\n", n_tri, max);
		}

		if (out_m >= 0 && abs(out_lvl) > 0) {
			if (n_tri%abs(out_lvl) == 0)
				Output(out_lvl, n_tri, max);
		}
	}
					// 最終出力
	if (out_m >= 0) {
		n_tri--;
		if (seisu > -2)
			printf("***試行回数 %d 距離 %d\n", n_tri, (int)max);
		else
			printf("***試行回数 %d 距離 %f\n", n_tri, max);
		Output(out_lvl, n_tri, max);
	}

	return n_tri;
}

/*******************************/
/* 出力                        */
/*      sw : >=0 : 出力先未定  */
/*           < 0 : ファイル    */
/*      n_tri : 現在の試行回数 */
/*      r : 距離               */
/*******************************/
void Iteration::Output(int sw, int n_tri, double r)
{
	int i1, k = 0, n, pr;
	char *now;
	time_t aclock;
	FILE *out;

	if (sw >= 0) {
		printf("   出力先は(0:出力なし,n:画面にn個づつ,-1:ファイル)? ");
		scanf("%d", &pr);
	}
	else
		pr = -1;

	if (pr != 0) {

		if (pr > 0) {
			out = stdout;
			getchar();
		}
		else {
			time(&aclock);
			now = ctime(&aclock);
			out = fopen(o_file, "a");
			if (seisu > -2)
				fprintf(out, "***試行回数 %d 距離 %d 時間 %s\n", n_tri, (int)r, now);
			else
				fprintf(out, "***試行回数 %d 距離 %d 時間 %s\n", n_tri, (int)(r+0.5), now);
		}

		if (out_m == 0) {
			for (i1 = 0; i1 < n_city; i1++) {
				n = seq[i1];
				if (seisu > 0)
					fprintf(out, "  %d %d %d\n", n, (int)city[n][0], (int)city[n][1]);
				else
					fprintf(out, "  %d %f %f\n", n, city[n][0], city[n][1]);
				if (pr > 0) {
					k++;
					if (k == pr) {
						getchar();
						k = 0;
					}
				}
			}
		}

		if (pr <= 0)
			fclose(out);
	}
}

/**************************************/
/* エッジの入れ替え                   */
/*      r_m : 距離                    */
/*      return : =0 : 改善がなかった  */
/*               =1 : 改善があった    */
/**************************************/
int Iteration::Change(double *r_m)
{
	double max, max1 = 0.0, r;
	int ch = 0, i0, i1, i2, i3, i4, k, k1 = 0, k2 = 0, k3, k4,
        n, nn, n1 = 0, n2 = 0, n3, n4, sw = 0, sw1 = 0, sw2;

	max = *r_m;

/*
     近傍を可変
*/
	if (fix == 0) {
					// 初期設定(k=2)
		k = 2;
		for (i1 = 0; i1 < n_city; i1++) {
			seq_w4[i1] = seq[i1];
			seq_w3[i1] = 0;
		}
						// 評価
		sw2 = 0;
		for (i0 = 0; i0 < n_city-2 && sw2 < 2; i0++) {

			n = (i0 == 0) ? n_city-1 : n_city;

			for (i1 = i0+2; i1 < n && sw2 < 2; i1++) {
							// 相手の場所
				k3 = i1;
				k4 = k3 + 1;
				if (k4 > n_city-1)
					k4 = 0;
							// 順番の入れ替え
				n3 = -1;
				for (i2 = 0; i2 < n_city && n3 < 0; i2++) {
					if (seq_w4[i2] == seq[i0+1])
						n3 = i2 + 1;
				}
				nn = n3;
				n4 = -1;
				for (i2 = 0; i2 < n_city && n4 < 0; i2++) {
					if (nn > n_city-1)
						nn = 0;
					if (seq_w4[nn] == seq[k3] || seq_w4[nn] == seq[k4])
						n4 = seq_w4[nn];
					else
						nn++;
				}
				if (n4 == seq[k4]) {
					n4 = k3;
					k3 = k4;
					k4 = n4;
				}
							// 評価
				seq_w1[0] = seq[k4];
				seq_w1[1] = seq[i0+1];
				n4        = -1;
				nn        = 2;
				while (n4 < 0) {
					if (n3 > n_city-1)
						n3 = 0;
					seq_w1[nn] = seq_w4[n3];
					if (seq_w4[n3] == seq[k3])
						n4 = 1;
					nn++;
					n3++;
				}
				seq_w1[nn] = seq[i0];
				nn++;
				n3 = -1;
				n4 = -1;
				for (i2 = 0; i2 < n_city && n3 < 0; i2++) {
					if (seq_w4[i2] == seq[i0]) {
						n3 = i2 - 1;
						if (n3 < 0)
							n3 = n_city - 1;
					}
				}
				while (n4 < 0) {
					if (seq_w4[n3] == seq[k4])
						n4 = 1;
					else {
						seq_w1[nn] = seq_w4[n3];
						nn++;
						n3--;
						if (n3 < 0)
							n3 = n_city - 1;
					}
				}

				r = kyori(n_city, seq_w1, rg);
							// 最適値の保存
				if (sw2 == 0 || r < max1) {
					sw2  = 1;
					max1 = r;
					n1   = k3;
					n2   = k4;
					k1   = i0;
					k2   = i0 + 1;
					for (i2 = 0; i2 < n_city; i2++)
						seq_w5[i2] = seq_w1[i2];
					if (sel > 0 && max1 < max)
						sw2 = 2;
				}
			}
		}
						// 最適値の保存と近傍の増加
		if (sw2 > 0) {
			if (max1 < max) {
				sw  = 1;
				max = max1;
				for (i1 = 0; i1 < n_city; i1++)
					seq_w2[i1] = seq_w5[i1];
			}
			if (k < neib) {
				for (i1 = 0; i1 < n_city; i1++)
					seq_w4[i1] = seq_w5[i1];
				seq_w3[k1] = 1;
				seq_w3[k2] = 1;
				seq_w3[n1] = 1;
				seq_w3[n2] = 1;
				k1         = n2;
				k++;
			}
			else
				sw1 = 1;
		}
		else
			sw1 = 1;
					// 実行(k>2)
		while (sw1 == 0) {
						// 評価
			sw2 = 0;
			for (i1 = 0; i1 < n_city; i1++) {
							// 相手の場所
				k3 = i1;
				k4 = k3 + 1;
				if (k4 > n_city-1)
					k4 = 0;

				if (seq_w3[k3] == 0 && seq_w3[k4] == 0) {
							// 順番の入れ替え
					n3 = -1;
					for (i2 = 0; i2 < n_city && n3 < 0; i2++) {
						if (seq_w4[i2] == seq[k2])
							n3 = i2 + 1;
					}
					nn = n3;
					n4 = -1;
					for (i2 = 0; i2 < n_city && n4 < 0; i2++) {
						if (nn > n_city-1)
							nn = 0;
						if (seq_w4[nn] == seq[k3] || seq_w4[nn] == seq[k4])
							n4 = seq_w4[nn];
						else
							nn++;
					}
					if (n4 == seq[k4]) {
						n4 = k3;
						k3 = k4;
						k4 = n4;
					}
							// 評価
					seq_w1[0] = seq[k4];
					seq_w1[1] = seq[k2];
					n4        = -1;
					nn        = 2;
					while (n4 < 0) {
						if (n3 > n_city-1)
							n3 = 0;
						seq_w1[nn] = seq_w4[n3];
						if (seq_w4[n3] == seq[k3])
							n4 = 1;
						nn++;
						n3++;
					}
					seq_w1[nn] = seq[k1];
					nn++;
					n3 = -1;
					n4 = -1;
					for (i2 = 0; i2 < n_city && n3 < 0; i2++) {
						if (seq_w4[i2] == seq[k1]) {
							n3 = i2 - 1;
							if (n3 < 0)
								n3 = n_city - 1;
						}
					}
					while (n4 < 0) {
						if (seq_w4[n3] == seq[k4])
							n4 = 1;
						else {
							seq_w1[nn] = seq_w4[n3];
							nn++;
							n3--;
							if (n3 < 0)
								n3 = n_city - 1;
						}
					}

					r = kyori(n_city, seq_w1, rg);
							// 最適値の保存
					if (sw2 == 0 || r < max1) {
						sw2  = 1;
						max1 = r;
						n1   = k3;
						n2   = k4;
						for (i2 = 0; i2 < n_city; i2++)
							seq_w5[i2] = seq_w1[i2];
					}
				}
			}
						// 最適値の保存と近傍の増加
			if (sw2 > 0) {
				if (max1 < max) {
					sw  = 1;
					max = max1;
					for (i1 = 0; i1 < n_city; i1++)
						seq_w2[i1] = seq_w5[i1];
				}
				if (k < neib) {
					for (i1 = 0; i1 < n_city; i1++)
						seq_w4[i1] = seq_w5[i1];
					seq_w3[n1] = 1;
					seq_w3[n2] = 1;
					k1         = n2;
					k++;
				}
				else
					sw1 = 1;
			}
			else
				sw1 = 1;
		}
	}
/*
     近傍を固定
*/
	else {
		n3  = (int)(genrand_real3() * (n_city - 2));
		if (n3 > n_city-3)
			n3 = n_city - 3;
                         // 2近傍
		for (i1 = 0; i1 <= n_city-3 && ch == 0; i1++) {

			if (n3 == 0)
				n1 = n_city - 2;
			else
				n1 = n_city - 1;

			for (i2 = n3+2; i2 <= n1 && ch == 0; i2++) {
                              // 枝の場所((n3,n3+1), (k1,k2))
				k1 = i2;
				if (i2 == n_city-1)
					k2 = 0;
				else
					k2 = i2 + 1;
                              // 枝の入れ替え
				seq_w1[0] = seq[n3];
				k         = 1;
				for (i3 = k1; i3 >= n3+1; i3--) {
					seq_w1[k] = seq[i3];
					k++;
				}

				nn = k2;
				while (nn != n3) {
					seq_w1[k] = seq[nn];
					k++;
					nn++;
					if (nn > n_city-1)
						nn = 0;
				}
                              // 評価
				r = kyori(n_city, seq_w1, rg);

				if (r < max) {
					max = r;
					sw  = 1;
					for (i3 = 0; i3 < n_city; i3++)
						seq_w2[i3] = seq_w1[i3];
					if (sel > 0)
						ch = 1;
				}
			}

			n3++;
			if (n3 > n_city-3)
				n3 = 0;
		}
                         // 3近傍
		if (neib == 3 && ch == 0) {

			for (i1 = 0; i1 <= n_city-3 && ch == 0; i1++) {

				n1 = n_city - 2;
				n2 = n_city - 1;

				for (i2 = n3+1; i2 <= n1 && ch == 0; i2++) {

					for (i3 = i2+1; i3 <= n2 && ch == 0; i3++) {
                              // 枝の場所((n3,n3+1), (i2,i2+1), (k1,k2))
						k1 = i3;
						if (i3 == n_city-1)
							k2 = 0;
						else
							k2 = i3 + 1;
                              // 枝の入れ替えと評価
                                   // 入れ替え(その1)
						seq_w1[0] = seq[n3];
						k         = 1;
						for (i4 = i2; i4 >= n3+1; i4--) {
							seq_w1[k] = seq[i4];
							k++;
						}

						for (i4 = k1; i4 >= i2+1; i4--) {
							seq_w1[k] = seq[i4];
							k++;
						}

						nn = k2;
						while (nn != n3) {
							seq_w1[k] = seq[nn];
							k++;
							nn++;
							if (nn > n_city-1)
								nn = 0;
						}
                                   // 評価(その1)
						r = kyori(n_city, seq_w1, rg);

						if (r < max) {
							max = r;
							sw  = 1;
							for (i3 = 0; i3 < n_city; i3++)
								seq_w2[i3] = seq_w1[i3];
							if (sel > 0)
								ch = 1;
						}
                                   // 入れ替え(その2)
						seq_w1[0] = seq[n3];
						k         = 1;
						for (i4 = k1; i4 >= i2+1; i4--) {
							seq_w1[k] = seq[i4];
							k++;
						}

						for (i4 = n3+1; i4 <= i2; i4++) {
							seq_w1[k] = seq[i4];
							k++;
						}

						nn = k2;
						while (nn != n3) {
							seq_w1[k] = seq[nn];
							k++;
							nn++;
							if (nn > n_city-1)
								nn = 0;
						}
                                   // 評価(その2)
						r = kyori(n_city, seq_w1, rg);

						if (r < max) {
							max = r;
							sw  = 1;
							for (i3 = 0; i3 < n_city; i3++)
								seq_w2[i3] = seq_w1[i3];
							if (sel > 0)
								ch = 1;
						}
                                   // 入れ替え(その3)
						seq_w1[0] = seq[n3];
						k         = 1;
						for (i4 = i2+1; i4 <= k1; i4++) {
							seq_w1[k] = seq[i4];
							k++;
						}

						for (i4 = i2; i4 >= n3+1; i4--) {
							seq_w1[k] = seq[i4];
							k++;
						}

						nn = k2;
						while (nn != n3) {
							seq_w1[k] = seq[nn];
							k++;
							nn++;
							if (nn > n_city-1)
								nn = 0;
						}
                                   // 評価(その3)
						r = kyori(n_city, seq_w1, rg);

						if (r < max) {
							max = r;
							sw  = 1;
							for (i3 = 0; i3 < n_city; i3++)
								seq_w2[i3] = seq_w1[i3];
							if (sel > 0)
								ch = 1;
						}
                                   // 入れ替え(その4)
						seq_w1[0] = seq[n3];
						k         = 1;
						for (i4 = i2+1; i4 <= k1; i4++) {
							seq_w1[k] = seq[i4];
							k++;
						}

						for (i4 = n3+1; i4 <= i2; i4++) {
							seq_w1[k] = seq[i4];
							k++;
						}

						nn = k2;
						while (nn != n3) {
							seq_w1[k] = seq[nn];
							k++;
							nn++;
							if (nn > n_city-1)
								nn = 0;
						}
                                   // 評価(その4)
						r = kyori(n_city, seq_w1, rg);

						if (r < max) {
							max = r;
							sw  = 1;
							for (i3 = 0; i3 < n_city; i3++)
								seq_w2[i3] = seq_w1[i3];
							if (sel > 0)
								ch = 1;
						}
					}
				}

				n3++;
				if (n3 > n_city-3)
					n3 = 0;
			}
		}
	}
                         // 設定
	if (sw > 0) {
		*r_m = max;
		for (i1 = 0; i1 < n_city; i1++)
			seq[i1] = seq_w2[i1];
	}

	return sw;
}

------------------------kyori.cpp-----------------
/*********************************/
/* 距離の計算                    */
/*      n_c : 都市の数           */
/*      p : 都市番号             */
/*      rg : 都市間の距離        */
/*      return : 距離            */
/*********************************/
#include <math.h>
#include <stdio.h>

float kyori(int n_c, int *p, float **rg)
{
	float range = 0;
	int i1, n1, n2;

	n1 = p[0];

	for (i1 = 1; i1 < n_c; i1++) {
		n2     = p[i1];
		range += rg[n1][n2];
		n1     = n2;
	}

	n2     = p[0];
	range += rg[n1][n2];

	return range;
}

------------------------main----------------------
/****************************/
/* 巡回セールスマン問題     */
/* (分割法)               */
/*      coded by Y.Suganuma */
/****************************/
#include "partition.h"

/****************/
/* main program */
/****************/
int main(int argc, char *argv[])
{
	double mean;
	long seed;
	int i0, i1, n, nm, max;
	char i_file[100];
	FILE *in, *out;
	Partition *pt;
					// 入力ミス
	if (argc <= 1) {
		printf("***error  ファイル名を入力して下さい\n");
		exit(1);
	}
					// 入力OK
	else {
						// ファイルのオープン
		in = fopen(argv[1], "r");
		if (in == NULL) {
			printf("***error  ファイル名が不適当です\n");
			exit(1);
		}
						// 入力データファイル名と問題数
		fscanf(in, "%*s %d", &nm);

		for (i0 = 0; i0 < nm; i0++) {
							// 各問題の実行
			fscanf(in, "%*s %s %*s %d", i_file, &n);
			pt   = new Partition(i_file);
			mean = 0.0;
			max  = -1;
								// 乱数の初期値を変える
			for (i1 = 0; i1 < n; i1++) {
									// 乱数の初期値
				fscanf(in, "%*s %ld", &seed);
				printf("\n+++++問題 %s 乱数 %ld+++++\n", i_file, seed);
									// 最適化
				pt->Optimize(seed);
									// 最適値とその平均の計算
				mean += pt->Max;
				if (max < 0 || pt->Max < max)
					max = pt->Max;
			}
							// 結果
			if (pt->out_m <= 0)
				printf("     -----最小 %d 平均 %f-----\n", max, mean/n);
			else {
				out = fopen(pt->o_file, "a");
				fprintf(out, "     -----最小 %d 平均 %f-----\n", max, mean/n);
				fclose(out);
			}
		}

		fclose(in);
	}

	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 */