#include <cstdio>
#include <complex>

int main()
{
	int x, t; 
	double p[101]; 
	std::complex<double> psi[101][2], old_psi[101][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 (x = 0; x <= 100; x++)
	{
		if (x == 0)
		{
			psi[x][0] = alpha; 
			psi[x][1] = beta; 
		}
		else
		{
			psi[x][0] = std::complex<double>(0, 0); 
			psi[x][1] = std::complex<double>(0, 0); 
		}
	}
	
	// 時刻100まで時間発展
	for (t = 0; t < 100; t++)
	{
		// 時刻tの確率振幅を配列old_psiにコピー
		for (x = 0; x <= 100; x++)
		{
			old_psi[x][0] = psi[x][0]; 
			old_psi[x][1] = psi[x][1]; 
			
			psi[x][0] = std::complex<double>(0, 0); // psi[x][0]の値を0にしておく
			psi[x][1] = std::complex<double>(0, 0); // psi[x][1]の値を0にしておく
		}
		
		// 時刻t+1の確率振幅（配列psiに格納）を計算
		for (x = 0; x <= 100; x++)
		{
			if (x == 0)
			{
				// 場所xのベクトルに行列P0を掛けて同じ場所に流出
				psi[x][1] += a * old_psi[x][0]; 
				psi[x][1] += b * old_psi[x][1]; 
				
				// 場所xのベクトルに行列Qを掛けて右に流出
				psi[x + 1][1] += c * old_psi[x][0]; 
				psi[x + 1][1] += d * old_psi[x][1]; 
			}
			else if (x == 100)
			{
				// 場所xのベクトルに行列Pを掛けて左に流出
				psi[x - 1][0] += a * old_psi[x][0]; 
				psi[x - 1][0] += b * old_psi[x][1]; 
			}
			else
			{
				// 場所xのベクトルに行列Pを掛けて左に流出
				psi[x - 1][0] += a * old_psi[x][0]; 
				psi[x - 1][0] += b * old_psi[x][1]; 
				
				// 場所xのベクトルに行列Qを掛けて右に流出
				psi[x + 1][1] += c * old_psi[x][0]; 
				psi[x + 1][1] += d * old_psi[x][1]; 
			}
		}
	}
	
	// 時刻100における確率（配列pに格納）を計算
	for (x = 0; x <= 100; x++)
	{
		p[x] = 0; // 初期化
		
		// 第0成分の大きさの2乗を加算
		p[x] += pow(psi[x][0].real(), 2) + pow(psi[x][0].imag(), 2); 
		
		// 第1成分の大きさの2乗を加算
		p[x] += pow(psi[x][1].real(), 2) + pow(psi[x][1].imag(), 2); 
		
		std::printf("p(%d) = %f\n", x, p[x]); 
	}
	
	return 0; 
}