using System;
using System.Globalization;
namespace NamespaceName
{
	public class ClassName
	{
		const int maxiter = 1000;
		const double eps = 1e-12;
		public static void PrintMatrix(int m,int n,double[,] A)
		{
			NumberFormatInfo nfi = new NumberFormatInfo();
			nfi.NumberDecimalDigits = 12;
			double p = 0;
			for(int i=0;i<m;i++)
			{
				for(int j=0;j < n; j++)
					if(Math.Abs(A[i,j]) > eps)
						Console.Write("{0} , ",A[i,j].ToString("N",nfi));
					else
						Console.Write("{0} , ",p.ToString("N",nfi)); 
				Console.WriteLine();
			}
			Console.WriteLine();
		}
		public static double Hypot(double a,double b)
		{
			double absa,absb;
			absa = Math.Abs(a);
			absb = Math.Abs(b);
			if(absa > absb) return absa * Math.Sqrt(1 + ((double) absb / absa) * ((double) absb / absa));
			else return (absb == 0 ? 0 : absb * Math.Sqrt(1 + ((double) absa / absb) * ((double) absa / absb)));
		}
		public static void ToHessenberg(int n,double[,] A)
		{
			for(int i = 0;i < n-2;i++)
			{
				int k = i + 1;
				double maxA = Math.Abs(A[k,i]);
				for(int j= i + 1; j < n; j++)
				   if(Math.Abs(A[j,i]) > maxA)
				   {
					   k  = j;
					   maxA = Math.Abs(A[k,i]);
				   }
				if(maxA > 0)
				{
					for(int j = 0; j < n; j++)
					{
						double temp = A[k,j];
						A[k,j] = A[i + 1,j];
						A[i + 1,j] = temp;
					}
					for(int j=0;j < n;j++)
					{
						double temp = A[j,k];
						A[j,k] = A[j, i + 1];
						A[j,i + 1] = temp;
					}
					for(int j = i + 2;j < n;j++)
					{
						double m = (double)A[j,i]/A[i + 1,i];
						for(int l = 0; l < n;l++)
							A[j,l] -= m * A[i + 1,l];
						for(int l = 0;l < n;l++)
							A[l,i + 1] += m * A[l,j];
					}
				}
			}
		}
		public static void QRdecomp(int m,int n,double[,] A,double[,] Q)
		{
			for(int i = 0; i < m; i++)
				for(int j = 0;j < m;j++)
					Q[i,j] = (i == j ? 1 : 0);
				int min = (m < n ? m : n);
				for(int i = 0; i < min; i++)
					for(int j = i + 1; j < m; j++)
						if(A[j,i] != 0)
						{
							double r = Hypot(A[i,i],A[j,i]);
							double c = (double) A[i,i] / r;
							double s = (double) A[j,i] / r;
							for(int k = 0; k < n; k++)
							{
								double temp = A[i,k];
								A[i,k] = c * A[i,k] + s * A[j,k];
								A[j,k] = -s * temp + c * A[j,k];
							}
							for(int k = 0;k < m;k++)
							{
								double temp = Q[k,i];
								Q[k,i] = c * Q[k,i] + s * Q[k,j];
								Q[k,j] = -s * temp + c * Q[k,j];
							}
						}
				
		}
 
		public static void MultiplyMatrix(int m,int n,int p,double[,] A,double[,] B,double[,] C)
		{
			for(int i = 0;i < m;i++)
				for(int k=0;k < p;k++)
				{
					C[i,k] = 0;
					for(int j=0;j < n;j++)
						C[i,k] += A[i,j] * B[j,k];
				}
		}
		public static void CopyMatrix(int m,int n,double[,] A,double[,] B)
		{
			for(int i=0;i<m;i++)
				for(int j=0;j<n;j++)
					B[i,j] = A[i,j];
		}
		public static void Main(string[] args)
		{
			int n,m;
			double[,] A;
			double[,] Q;
			double[,] R;
			Console.WriteLine("Obliczanie przyblizonych wartosci wlasnych metoda rozkladu QR");
				Console.WriteLine("Podaj wymiar macierzy ktorej wartosci wlasne chcesz policzyc");
				int.TryParse(Console.ReadLine(),out n);
				Console.WriteLine("Podaj gorny zakres przedzialu do ktorego beda nalezec elementy macierzy");
				int.TryParse(Console.ReadLine(),out m);
				A = new double[n,n];
				Q = new double[n,n];
				R = new double[n,n];
				Random rand = new Random();
				for(int i=0;i<n;i++)
					for(int j = 0;j < n;j++)
						A[i,j] = (1 - 2 * rand.Next(2)) * rand.Next(m);
				Console.WriteLine("Wylosowana macierz");
				PrintMatrix(n,n,A);
				ToHessenberg(n,A);
				Console.WriteLine("Macierz po sprowadzeniu do postaci Hessenberga");
				PrintMatrix(n,n,A);
				for(int i = 0;i<maxiter;i++)
				{
					QRdecomp(n,n,A,Q);
					CopyMatrix(n,n,A,R);
					MultiplyMatrix(n,n,n,R,Q,A);
				}
				Console.WriteLine("Macierz po iteracyjnym sprowadzeniu do postaci Schura");
				PrintMatrix(n,n,A);
				int k = 0;
				NumberFormatInfo nfi = new NumberFormatInfo();
				nfi.NumberDecimalDigits = 12;
				while(k < n)
				{
					if(k + 1 < n && Math.Abs(A[k+1,k]) > eps)
					{
						/* Tutaj zalozylem ze w postaci Schura klatka 2x2 pojawia sie
						 na diagonali gdy macierz posiada zespolone wartosci wlasne
						 */
						double p = 0.5 * (A[k,k] + A[k + 1,k + 1]);
						double q = A[k,k] * A[k+1,k+1] - A[k,k+1] * A[k+1,k];
						double d = Math.Sqrt(Math.Abs(q - p * p));
						Console.WriteLine("x[{0}] = {1} - {2}*i",k,p.ToString("N",nfi),d.ToString("N",nfi));
						Console.WriteLine("x[{0}] = {1} + {2}*i",k+1,p.ToString("N",nfi),d.ToString("N",nfi));
						k += 2;
					}
					else
					{
						Console.WriteLine("x[{0}] = {1}",k,A[k,k].ToString("N",nfi));
						k++;
					}
				}		
		}
	}
}