#include <cstdio>
#include <complex>

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