////////////////////////  以下プログラム  ////////////////////////////
/*
 * 現状の安定バージョン
 * --->F2の結果に関する有効桁数が不明。
 * http://radhia.blog23.fc2.com/
 * 作成者：ラディア
 */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <math.h>

// MAX_KEISUU : 係数の最大個数
#define MAX_KEISUU (14+1)

// MAX_C : 染色体の最大値
#define MAX_C 7



// WHAT_NUMBER : 求める数の10進数表示
//#define WHAT_NUMBER 10000
#define WHAT_NUMBER 10000

// PENARTY_NUMBER : ペナルティ値
#define PENARTY_NUMBER 1000000

// GA -----------------------------------------------------------------

// POPULATION_COUNT : 集団の大きさ、個体数
// デフォルトは500
#define POPULATION_COUNT 500

// CUT_LINE : 足切りライン。この順位以下の固体は淘汰される。
#define CUT_LINE 50

// MUTATION_ODDS : 突然変異を起こす確率
#define MUTATION_ODDS 0.05

/*
INDIVIDUAL : 個体（染色体）。個体は係数を示す経路(pathway)で表現され、
             その経路から算出されるfitnessで評価される
			 今回に限ってはfitnessは小さいほどいいことにする
*/
typedef struct {
	double pathway[MAX_KEISUU - 1];
	double keisuu_t[MAX_KEISUU - 1];
	double fitness1;
	double fitness2;
} INDIVIDUAL;

// population : 個体の集合
INDIVIDUAL **population;



/* p[n]からt[n]を算出する
 * 
 * 関数として、t[n]とp[n]の比較を行い、
 * t[n] < p[n]であれば0を返す。
 */
void calculate_t(INDIVIDUAL *target) {
	int number_x = WHAT_NUMBER;
	int i;
	
	for(i=0; i<MAX_KEISUU-1; i++){
		if(target->pathway[i] > 0) {
			target->keisuu_t[i] = fmod(number_x, target->pathway[i]); // 剰余をtarget->keisuu_t[i]に入れる
			if(target->keisuu_t[i] < 0)
				//target->keisuu_t[i] = 0;
				printf("target->keisuu_t[%2d] = %d がマイナス\n", i, target->keisuu_t[i]);
		}
		else {
			//target->keisuu_t[i] = 0;
			//printf("target->pathway[%2d] = %d がマイナス\n", i, target->pathway[i]);
		//printf("target->keisuu_t[%2d] = %f\n", i, target->keisuu_t[i]);
		}
		if(target->pathway[i] > 0)
			number_x = number_x / target->pathway[i]; // 商を次のnumber_xとする
		else {
			number_x = number_x;
			printf("Error!target->pathway[%2d] = 0\n", i);
		}
	}
}

// 個体のF1の計算を行う
void calculate_fitness1(INDIVIDUAL *target) {
	double return_fitness1 = 0;
	int i;

	// fitnessを加算
	/* 	 
	 * return_fitness1 = Σ((2*p[i] + 3*t[i] -2)a)
	 */
	
	calculate_t(target);
	
	for(i=0; i< MAX_KEISUU-1; i++)
		return_fitness1 += (2 * target->pathway[i] + 3 * target->keisuu_t[i] -2);

	target->fitness1 = return_fitness1;

}

// calculate_fitness2で用いるd[i]の算出を行う
void calculate_d(int d[], INDIVIDUAL *target) {
	int i;
	
	for(i=0; i<MAX_KEISUU -1; i++){
		d[i] = (target->keisuu_t[i] > 0)? 1 : 0;
	}
}

// 個体のF2の計算を行う
void calculate_fitness2(INDIVIDUAL *target) {
	double return_fitness2 = 0;
	int i, j, d[MAX_KEISUU];

	// fitnessを加算
	/* 	 
	 * return_fitness2 = Σ((p[i] + d[i] - 1)(2b/p[i] + c/(p[i] + d[i]))/(Πp[j])
	 * 
	 * まずは b = c = 1 でやる
	 */
	calculate_d(d, target);
	
	for(i=0; i< MAX_KEISUU - 1; i++){
		return_fitness2 +=
			(target->pathway[i] + d[i] - 1) * (2 / target->pathway[i] + 1 / (target->pathway[i] + d[i]));

		//	ここで0になる
		for(j=0; j <= i - 1; j++) {
			if(target->pathway[j] > 0)
				return_fitness2 /= target->pathway[j];
			else {
				if(target->pathway[i] < 0) {
					printf("target->pathway[%2d] = %fにつき終了\n", i, target->pathway[i]);
					exit(1);
				}
			}
			//printf("return_fitness2 = %d, target->pathway[%2d] = %d\n", return_fitness2, j, target->pathway[j]);
		}
		
	}
	target->fitness2 = return_fitness2;
}

/* 遺伝子の残りの部分を0で埋める */
void buries_generation(int num, INDIVIDUAL *target) {
	int i;

	for(i = num; i < MAX_KEISUU - 1; i++) {
		target->keisuu_t[i] = 0;
		target->pathway[i]  = 1;
	}
}

/* WHAT_NUMBERと一致しているかのチェック */
void check_generation_x(INDIVIDUAL *target) {
	int i, j;
	double sum = 0;
	
	for(i=MAX_KEISUU - 1; i > 0; i--) {
		if(target->pathway[i] < 0) {
			printf("target->pathway[%2d] = %fにつき終了\n", i, target->pathway[i]);
			exit(1);
		}
		sum += sum * target->pathway[i] + target->keisuu_t[i];
		if(sum == WHAT_NUMBER) {
			buries_generation(i, target);
			
			 // これが表示されないということは、WHAT_NUMBERを満たすような遺伝子が存在していないということになる！
			printf("WHAT_NUMBERを満たすような遺伝子が存在しています\n");
		}
	}
	if(sum > WHAT_NUMBER)
		target->fitness1 = PENARTY_NUMBER; // Error!遺伝子が条件を満たしていない(WHAT_NUMBERを超えている)
	if(sum < WHAT_NUMBER)
		target->fitness1 = PENARTY_NUMBER; // Error!遺伝子が条件を満たしていない(WHAT_NUMBERに達していない)
}


/*
 * 個体の突然変異体を作る
 * mutation(INDIVIDUAL *target)
 * mutation2(INDIVIDUAL *target)
 * mutation3(INDIVIDUAL *target)
 * increment(INDIVIDUAL *target)
 * decrement(INDIVIDUAL *target)
 */
// 染色体の2箇所をランダムに選んで入れ替え
void mutation(INDIVIDUAL *target) {
	int i, j, tmp;
	
	i = rand() % (MAX_KEISUU - 1);
	j = rand() % (MAX_KEISUU - 1);
	tmp = target->pathway[j];
	target->pathway[j] = target->pathway[i];
	target->pathway[i] = tmp;
	check_generation_x(target);

}

// 染色体の2箇所をランダムに選んでこの間の数を乱数化
void mutation2(INDIVIDUAL *target) {
	int i, j, k, tmp;
	
	i = rand() % (MAX_KEISUU - 1);
	j = rand() % (MAX_KEISUU - 1);
	if(i > j) {
		tmp = i;
		i = j;
		j = tmp;
	}
	//printf("i:%d, j:%d\n", i, j);
	for(k=i; k<=j; k++) {
		/*
		 * p[s]は全て2以上7以下の整数
		 */
		target->pathway[k] = rand() % (MAX_C - 1) + 2;
	}
	check_generation_x(target);
}

// 染色体をすべて乱数化
void mutation3(INDIVIDUAL *target) {
	int i, j, k, tmp;
	
	i = 0;
	j = MAX_KEISUU - 1;

	//printf("i:%d, j:%d\n", i, j);
	for(k=i; k<=j; k++) {
		//target->pathway[k] = rand() % (MAX_KEISUU - 1);
		/*
		 * p[s]は全て2以上7以下の整数
		 */
		target->pathway[k] = rand() % (MAX_C - 1) + 2;
		//printf("target->pathway[%d] = %d\n", k, target->pathway[k]);
	}
	check_generation_x(target);
}

// 染色体の2箇所をランダムに選んで１インクリメント
void increment(INDIVIDUAL *target) {
	int i, j, tmp;
	
	i = rand() % (MAX_KEISUU - 1);
	if(target->pathway[i] < MAX_C)
		target->pathway[i] = (target->pathway[i] + 1);
	check_generation_x(target);
}

// 染色体の2箇所をランダムに選んで2インクリメント
/*void double_increment(INDIVIDUAL *target) {
	int i, j, tmp;
	
	i = rand() % (MAX_KEISUU - 1);
	target->pathway[i] = (target->pathway[i] + 2);
	check_generation_x(target);
}*/

// 染色体の2箇所をランダムに選んで１デクリメント
void decrement(INDIVIDUAL *target) {
	int i, j, tmp;
	
	i = rand() % (MAX_KEISUU - 1);
	if(target->pathway[i] > 2 && target->pathway[i] <= MAX_C)
		target->pathway[i] = (target->pathway[i] - 1);
	check_generation_x(target);
}

// 染色体の2箇所をランダムに選んで2デクリメント
/*void double_decrement(INDIVIDUAL *target) {
	int i, j, tmp;
	
	i = rand() % (MAX_KEISUU - 1);
	target->pathway[i] = (target->pathway[i] - 2);
	check_generation_x(target);
}*/

/*
// 2つの個体から一点交叉によって新しい個体を得る
void crossover(INDIVIDUAL *target, INDIVIDUAL *source1, INDIVIDUAL *source2) {
	int used_flag[MAX_KEISUU - 1], change_pos[(MAX_KEISUU -1) / 2], point,
		i, j, k;
	// 交叉位置を乱数で決定
	point = MAX_KEISUU / 2 + rand() % ((MAX_KEISUU -1) / 2);
	// 染色体の前半をsource1からコピー
	memcpy(target->pathway, source1->pathway, sizeof(int) * point);
	// 染色体の後半をsource2からコピー
	memcpy(target->pathway + point, source2->pathway+point,
		sizeof(int) * (MAX_KEISUU - 1 - point));
}*/

// 2つの個体から2点交叉によって新しい個体を得る
void crossover(INDIVIDUAL *target, INDIVIDUAL *source1, INDIVIDUAL *source2) {
	int change_pos1[(MAX_KEISUU -1) / 2], change_pos2[(MAX_KEISUU -1) / 2],
	point1, point2, i, j, k, tmp;
	// 交叉位置を乱数で決定
	//point1 = MAX_KEISUU / 2 + rand() % ((MAX_KEISUU -1) / 2);
	//point2 = MAX_KEISUU / 2 + rand() % ((MAX_KEISUU -1) / 2);
	point1 = rand() % (MAX_KEISUU -1) + 1;
	point2 = rand() % (MAX_KEISUU -1) + 1;
	//printf("point1:%d, point2:%d\n", point1, point2);
	if(point1 < point2) {
		tmp = point1;
		point1 = point2;
		point2 = tmp;
	}
	
	// 染色体の前半をsource1からコピー
	memcpy(target->pathway, source1->pathway, sizeof(int) * point1);
	// 染色体の真ん中をsource2からコピー
	memcpy(target->pathway + point1, source2->pathway+point2,
		sizeof(int) * (MAX_KEISUU - 1 - point1));
	// 染色体の後半をsource1からコピー
	memcpy(target->pathway + point2, source1->pathway+point2,
		sizeof(int) * (MAX_KEISUU - 1 - point2));
	check_generation_x(target);
}

// GAに関するすべての操作の前に初期化を行う
void ga_init(void) {
	int i;
	INDIVIDUAL tmp;
	// 集団のメモリ領域を確保
	population = (INDIVIDUAL **)malloc(sizeof(INDIVIDUAL *) * POPULATION_COUNT);
	for(i = 0; i < POPULATION_COUNT; i++)
		population[i] = (INDIVIDUAL *)malloc(sizeof(INDIVIDUAL));
	// 集団の初期状態を生成する
	
	srand((unsigned) time(NULL));
	
	for(i = 0; i < MAX_KEISUU - 1; i++) {
		tmp.keisuu_t[i] = 0; // 初期化
		
		//		tmp.pathway[i] = i + 1;
		/*
		 * p[s]は全て2以上7(= MAX_C)以下の整数
		 */
		tmp.pathway[i] = rand() % (MAX_C - 1) + 2;
		//printf("tmp.pathway[%d] = %d\n", i, tmp.pathway[i]);

		//printf("tmp.pathway[%2d] = %f\n", i, tmp.pathway[i]);
	}
	check_generation_x(&tmp);
	for(i = 0; i < POPULATION_COUNT; i++) {
		memcpy(population[i], &tmp, sizeof(INDIVIDUAL));
		mutation(&tmp);
	}
}

// 各個体のfitnessを計算し昇順に並べ替える
void ga_ranking() {
	int i, j;
	INDIVIDUAL *tmp;
	// 各個体のfitnessを計算
	for(i = 0; i < POPULATION_COUNT; i++){
		
		calculate_fitness1(population[i]);
		calculate_fitness2(population[i]);
	}
	/* fitnessに基づき並べ替える
	 * fitness1及びfitness2の小さい個体ほど上位にくる
	 */

	for(i = 0; i < POPULATION_COUNT - 1; i++) {
		for(j = i + 1; j < POPULATION_COUNT; j++) {
			//if(population[i]->fitness1 < population[j]->fitness1) {
			//if(population[i]->fitness1 + population[i]->fitness2 > population[j]->fitness1 + population[j]->fitness2)
			if(population[i]->fitness1 > population[j]->fitness1 && population[i]->fitness2 > population[j]->fitness2) {
				//printf("population[i] = %d, population[j] = %d\n", population[i], population[j]);
				tmp = population[j];
				population[j] = population[i];
				population[i] = tmp;
				//printf("population changed.\n");
			}
		}
	}
	//	for(i = 0; i < POPULATION_COUNT; i++)
	//	  printf("fit[%d]. %g,%d\n", i, population[i]->fitness, calculate_gt(population[i]));
}

// 次世代の集団を形成する
void ga_next_generation(void) {
	int i;
	// 集団は順位づけられているものとし、CUT_LINE以下の個体を
	// CUT_LINE以上の個体から一点交叉によって生成される個体に置き換える
	for(i = CUT_LINE; i < POPULATION_COUNT; i++)
		crossover(population[i], population[rand() % CUT_LINE],
		population[rand() % CUT_LINE]);
	//printf("ga_next_generation()\n");
	// すべての個体はMUTATION_ODDSの確率で突然変異を起こす
	// 一部の固体はMUTATION_ODDS/2の確率で突然変異（特に変化の激しいもの）
	for(i = 0; i < POPULATION_COUNT; i++) {
		if((rand() % (int)(1.0 / MUTATION_ODDS)) == 0)
			mutation(population[i]);
		else if((rand() % (int)(1.0 / MUTATION_ODDS)) == 0)
			increment(population[i]);
		/*else if((rand() % (int)(1.0 / MUTATION_ODDS)) == 0)
			double_increment(population[i]);*/
		else if((rand() % (int)(1.0 / MUTATION_ODDS)) == 0)
			decrement(population[i]);
		/*else if((rand() % (int)(1.0 / MUTATION_ODDS)) == 0)
			double_decrement(population[i]);*/
		/*else if((rand() % (int)(1.0 / MUTATION_ODDS)) == 0)
			mutation(population[i]);*/
		else if((rand() % (int)(1.0 / MUTATION_ODDS / 2)) == 0)
			mutation2(population[i]);
		else if((rand() % (int)(1.0 / MUTATION_ODDS / 2)) == 0)
			mutation3(population[i]);
	}
}

// メイン関数
int main(void) {
	int i, j, k, not_change_count;
	time_t now_t, old_t;
	INDIVIDUAL ace;
	
	// 処理時間の計測
		time(&old_t);
		
	// srand()を導入すると、解が毎回異なるようになる。
//	srand((unsigned) time(NULL));



	//printf("%d進数表示問題を解く \n", MAX_KEISUU);
	// GAによる近似解の探索
	
	
	
	// 初期化
	ga_init();
	// aceのfitness1としてありえない大きさの値を初期値とする
	ace.fitness1 = PENARTY_NUMBER+1;
	ace.fitness2 = PENARTY_NUMBER+1;
	// aceが10世代交代しない場合に終了
	for(i = 0, not_change_count = 0; not_change_count < 10;
		ga_next_generation(), i++, not_change_count++) {
			ga_ranking();
			/* 
			 * 過去のaceより優秀なら、新たなaceとして採用する。
			 * ただしfitnessは0以上が原則
			 */
			//if(ace.fitness < population[0]->fitness) {
			if(ace.fitness1 > population[0]->fitness1 && population[0]->fitness1 >= 0
				&& ace.fitness2 > population[0]->fitness2 && population[0]->fitness2 >= 0) {
				memcpy(&ace, population[0], sizeof(INDIVIDUAL));
				not_change_count = 0;

				printf("ace.fitness1 = %3f, ace.fitness2 = %.30f\n", ace.fitness1, ace.fitness2);
				for(j=0; j<MAX_KEISUU-1;j++)
					//printf("ace.keisuu_t[%2d] = %d\n", j, ace.keisuu_t[j]);
				for(j=0; j<MAX_KEISUU-1;j++)
					//printf("ace.pathway[%2d] = %d\n", j, ace.pathway[j]);
				//printf("%3d\n", i);
				//printf("ace changed.\n");
				;
			}
	}
	// 結果表示
	/*puts("GA");
	for(i = 0; i < MAX_KEISUU - 1; i++) {
		//printf("a[%d] = %2d, ", i, population[0]->pathway[i]);
		printf("%2dt^%d+ ", population[0]->pathway[i], i);
	}
	printf("\n");*/
	//printf("fitness(k=%2d): %f\n", k, population[0]->fitness);
	//printf("%d\n", population[0]->fitness1);
	//printf("examination %d\n", calculate_gt(population[0]));

	// 処理時間の計測
	time(&now_t);
	
	printf("%d sec\n", now_t - old_t);
	return 0;
}