// 血液型相性シミュレーター (2007/01/11)

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// GSL乱数の設定

#define SEED 777
const gsl_rng_type *T;
gsl_rng *r;
int seed;
// #define FRAND() (rand()/(RAND_MAX+1.0))


void marriage(int m, int f, int na, int *M, int *F, double *g)
{
	// 結婚相談所
	//    引数: m-男性ナンバー, f-女性ナンバー, na-血液型の種類, M-男性血液型配列, F-女性血液型配列, g-相性配列

	int i;
	double tmp;
	int tm[2],tf[2],child[2];
	int ttm,ttf,t1,t2,t;

	// 相性により結婚するかしないかの判断

	if (gsl_rng_uniform_pos(r)>g[M[m]*na+F[f]])
	{
		// 結婚しない

		return;
	}

	// 血液型遺伝要素の抽出

	switch(M[m])
	{
		case(0):
		{
			tm[0]=0;  tm[1]=0;
			break;
		}
		case(1):
		{
			tm[0]=0;  tm[1]=1;
			break;
		}
		case(2):
		{
			tm[0]=1;  tm[1]=1;
			break;
		}
		case(3):
		{
			tm[0]=1;  tm[1]=2;
			break;
		}
		case(4):
		{
			tm[0]=2;  tm[1]=2;
			break;
		}
		case(5):
		{
			tm[0]=0;  tm[1]=2;
			break;
		}
	}
	switch(F[f])
	{
		case(0):
		{
			tf[0]=0;  tf[1]=0;
			break;
		}
		case(1):
		{
			tf[0]=0;  tf[1]=1;
			break;
		}
		case(2):
		{
			tf[0]=1;  tf[1]=1;
			break;
		}
		case(3):
		{
			tf[0]=1;  tf[1]=2;
			break;
		}
		case(4):
		{
			tf[0]=2;  tf[1]=2;
			break;
		}
		case(5):
		{
			tf[0]=0;  tf[1]=2;
			break;
		}
	}

	// 子供を作る

	for (i=0 ; i<2 ; i++)
	{
		// 減数分裂

		if (gsl_rng_uniform_pos(r)<0.5) ttm=0; else ttm=1;
		if (gsl_rng_uniform_pos(r)<0.5) ttf=0; else ttf=1;

		// 子供の血液型の設定

		if (tm[ttm]<=tf[ttf])
		{
			t1=tm[ttm];  t2=tf[ttf];
		}
		else
		{
			t1=tf[ttf];  t2=tm[ttm];
		}
		t=t1*3+t2;
		switch(t)
		{
			case(0):
			{
				// AA型

				child[i]=0;
				break;
			}
			case(1):
			{
				// AO型

				child[i]=1;
				break;
			}
			case(2):
			{
				// AB型

				child[i]=5;
				break;
			}
			case(4):
			{
				// OO型

				child[i]=2;
				break;
			}
			case(5):
			{
				// BO型

				child[i]=3;
				break;
			}
			case(8):
			{
				// BB型

				child[i]=4;
				break;
			}
		}
	}
	M[m]=child[0];
	F[f]=child[1];
}


void quicksort(double* a, int* list, int left,int right)
{
	// 任意精度配列のクイックソート(小さいもん順)
	//    a:配列, list:index, left:左値index, right:右値index

	int i,j,p,q;
	double s,t;

	if (left<right)
	{
		s=a[(left+right)/2];
		p=list[(left+right)/2];
		i=left-1;
		j=right+1;
		while (true)
		{
			while(a[++i]<s);
			while(a[--j]>s);
			if (i>=j) break;
			t=a[i];  a[i]=a[j];  a[j]=t;
			q=list[i];  list[i]=list[j];  list[j]=q;
		}
		quicksort(a,list,left,i-1);
		quicksort(a,list,j+1,right);
	}
}


int main(int argc, char *argv[])
{
	int i,j,k,is;

	FILE *fno;
	char buff[512];
	int buff_size=sizeof(buff);

	int nstep;         // 交配のステップ数
	int na=6;          // 血液型の種類
	int tn;            // 男性、女性それぞれの総人口
	int *n;            // AA,AO,OO,BO,BB,ABの人口
	double *g;         // AA-AA,AA-AO,...の相性
	int *M;            // 男性の血液型
	int *F;            // 女性の血液型
	int *ind;          // インデックス配列
	double *ran;       // ソート用乱数導入配列

	// GSL乱数の設定

	T=gsl_rng_default;
	r=gsl_rng_alloc(T);
//	seed = time(NULL);
	seed=SEED;
	gsl_rng_set(r,seed);

	// 最初のメモリ確保

	n=new int[na];
	g=new double[na*na];

	// インプット(ファイル:input.txt)
	//
	//    交配のステップ数
	//
	//    初期AA型人口 (男性、女性片方の数を指定, 男女は同数になる)
	//    初期AO型人口
	//    初期OO型人口
	//    初期BO型人口
	//    初期BB型人口
	//    初期AB型人口
	// 
	//    AA-AAの相性 (0:絶対結婚しません,1:絶対結婚します, 0-1の範囲の数値)
	//    AA-AOの相性
	//    AA-OOの相性
	//    AA-BOの相性
	//    AA-BBの相性
	//    AA-ABの相性
	//
	//    AO-AOの相性
	//    AO-OOの相性
	//    AO-BOの相性
	//    AO-BBの相性
	//    AO-ABの相性
	//
	//    OO-OOの相性
	//    OO-BOの相性
	//    OO-BBの相性
	//    OO-ABの相性
	//
	//    BO-BOの相性
	//    BO-BBの相性
	//    BO-ABの相性
	//
	//    BB-BBの相性
	//    BB-ABの相性
	//
	//    AB-ABの相性

	fno=fopen("input.txt","rt");
	fgets(buff,buff_size,fno);  nstep=atoi(buff);
	for (i=0 ; i<na ; i++)
	{
		fgets(buff,buff_size,fno);  n[i]=atoi(buff);
	}
	for (i=0 ; i<na ; i++)
		for (j=i ; j<na ; j++)
		{
			fgets(buff,buff_size,fno);  g[i*na+j]=atof(buff);  g[j*na+i]=g[i*na+j];
		}
	fclose(fno);

	// 総人口の設定

	tn=0;
	for (i=0 ; i<na ; i++)
		tn+=n[i];

	// 血液型配列のメモリ確保

	M=new int[tn];
	F=new int[tn];
	ind=new int[tn];
	ran=new double[tn];

	// 血液型情報の初期化

	k=0;
	for (i=0 ; i<na ; i++)
		for (j=0 ; j<n[i] ; j++)
		{
			M[k]=i;
			F[k]=i;
			k++;
		}

	// アウトプットファイルの初期化

	fno=fopen("output.txt","wt");
	fclose(fno);

	// シミュレーション開始

	for (is=0 ; is<nstep ; is++)
	{
		// インデックス配列を初期化

		for (i=0 ; i<tn ; i++)
			ind[i]=i;

		// ソート用乱数配列に乱数を代入

		for (i=0 ; i<tn ; i++)
			ran[i]=gsl_rng_uniform_pos(r);

		// 結婚相手を選ぶためのソーティング

		quicksort(ran,ind,0,tn-1);

		// 交配開始

		for (i=0 ; i<tn ; i++)
			marriage(ind[i],i,na,M,F,g);

		// 血液型のチェック(男性、女性の総数)

		for (i=0 ; i<na; i++)
			n[i]=0;
		for (i=0 ; i<tn ; i++)
		{
			n[M[i]]++;
			n[F[i]]++;
		}

		// ファイルに出力

		fno=fopen("output.txt","at");
		for (i=0 ; i<na ; i++)
			fprintf(fno,"%10d",n[i]);
		fprintf(fno,"\n");
		fclose(fno);
	}

	// メモリ開放

	delete[] n;
	delete[] g;
	delete[] M;
	delete[] F;
	delete[] ind;
	delete[] ran;

	return 0;
}