#include <cstdio>
#include <complex>

int main()
{
	int j, k, l, t; 
	double p[20]; 
	std::complex<double> psi[20][20], old_psi[20][20], u[20][20]; 
	
	// ユニタリ行列の設定
	for (j = 0; j < 20; j++)
	{
		for (k = 0; k < 20; k++)
		{
			if (k == j)  u[j][k] = std::complex<double>(2.0 / 20 - 1.0, 0); 
			else  u[j][k] = std::complex<double>(2.0 / 20, 0); 
		}
	}
	
	// 初期確率振幅
	for (j = 0; j < 20; j++)
	{
		if (j == 0)
		{
			for (k = 0; k < 20; k++)
			{
				if (k == 0)  psi[j][k] = std::complex<double>(1.0, 0); 
				else  psi[j][k] = std::complex<double>(0, 0); 
			}
		}
		else
		{
			for (k = 0; k < 20; k++)	psi[j][k] = std::complex<double>(0, 0); 
		}
	}
	
	// 時刻100まで時間発展
	for (t = 0; t < 100; t++)
	{
		// 時刻tの確率振幅を配列old_psiにコピー
		for (j = 0; j < 20; j++)
		{
			for (k = 0; k < 20; k++)	old_psi[j][k] = psi[j][k]; 
			
			for (k = 0; k < 20; k++)	psi[j][k] = std::complex<double>(0, 0); 
		}
		
		// 時刻t+0.5の確率振幅（配列psiに格納）を計算
		for (k = 0; k < 20; k++)
		{
			for (j = 0; j < 20; j++)
			{
				for (l = 0; l < 20; l++)	psi[k][j] += u[j][l] * old_psi[k][l]; 
			}
		}
		
		// 時刻t+0.5の確率振幅を配列old_psiにコピー
		for (j = 0; j < 20; j++)
		{
			for (k = 0; k < 20; k++)	old_psi[j][k] = psi[j][k]; 
			
			for (k = 0; k < 20; k++)	psi[j][k] = std::complex<double>(0, 0); 
		}
		
		// 時刻t+1の確率振幅（配列psiに格納）を計算
		for (k = 0; k < 20; k++)
		{
			for (j = 0; j < 20; j++)	psi[j][k] += old_psi[k][j]; 
		}
	}
	
	// 時刻100における確率（配列pに格納）を計算
	for (j = 0; j < 20; j++)
	{
		p[j] = 0; 
		for (k = 0; k < 20; k++)
		{
			p[j] += pow(psi[j][k].real(), 2) + pow(psi[j][k].imag(), 2); 
		}
		
		std::printf("p(%d) = %f\n", j, p[j]); 
	}

	return 0; 
}