C#

using System;
using System.IO;
using System.Globalization;
namespace NamespaceName
{
	public class ClassName
	{
		const int maxiter = 1000;
		const double eps = 1e-12;
		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 printMatrix(int m,int n, double[,] A)
		{
			NumberFormatInfo nfi = new NumberFormatInfo();
			nfi.NumberDecimalDigits = 12;
			for(int i = 0;i < m;i++)
			{
				for(int j = 0; j < n; j++)
				{
					Console.Write("{0} , ",A[i,j].ToString("N",nfi));
				}
				Console.WriteLine();
			}
			Console.WriteLine();
		}
		public static void QR_Givens(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 copyMatrix(double[,] A,double[,] B,int m,int n)
		{
			for(int i = 0;i < m;i++)
				for(int j = 0;j < n;j++)
					B[i,j] = A[i,j];
		}
		public static void multiplyMatrix(double[,] A,double[,] B,double[,] C,int m,int n,int p)
		{
			for(int i = 0;i < m;i++)
				for(int j = 0;j < p;j++)
				{
					double sum = 0;
					for(int k = 0;k < n;k++)
						sum += A[i,k] * B[k,j];
					C[i,j] = sum;
				}
		}
		public static double rayleigh(int n, double[,] A,double[] q)
		{
			double norm = 0;
			double sum;
			double[] v = new double[n];
			for(int i = 0;i < n;i++)
				norm += q[i] * q[i];
			for(int i = 0;i < n;i++)
			{
				sum = 0;
				for(int j = 0;j < n;j++)
					sum += A[i,j]*q[j];
				v[i] = sum;
			}
			double r = 0;
			for(int i = 0;i < n;i++)
				r += q[i]*v[i];
			r /= (double)norm;
			return r;
		}
		public static void Main(string[] args)
		{
			char esc;
			int n;
			double[,] A,Q,R;
			double[] v;
			double[] a;
			double r;
			Random rnd = new Random();
			NumberFormatInfo nfi = new NumberFormatInfo();
			nfi.NumberDecimalDigits = 12;
			
			
			using(StreamWriter sw = new StreamWriter("polyroots.txt",true))
			{
				do
				{
				Console.WriteLine("Podaj stopien wielomianu");
				int.TryParse(Console.ReadLine(),out n);
				A = new double[n,n];
				Q = new double[n,n];
				R = new double[n,n];
				v = new double[n];
				a = new double[n + 1];
				for(int i = n;i >= 0;i--)
				{
					Console.Write("Podaj a[{0}]= ", i);
					double.TryParse(Console.ReadLine(),out a[i]);
				}
				for(int i = n;i >= 0;i--)
					if(a[i] < 0)
						sw.WriteLine("-{0}x^{1} ",(-a[i]).ToString("N",nfi),i);
					else
						sw.WriteLine("+{0}x^{1} ",a[i].ToString("N",nfi),i);
					sw.WriteLine();	
				for(int i=0;i<n;i++)
					A[0,i] = (double)(-a[n-i-1]/a[n]);
				for(int i = 1;i < n;i++)
					for(int j = 0;j<n;j++)
						A[i,j] = (i == j+1)?1:0;
				printMatrix(n,n,A);
				for(int i = 0;i < n;i++)
				{
					for(int j = 0; j < n; j++)
						sw.Write("{0} ",A[i,j].ToString("N",nfi));
					sw.WriteLine();
				}
				sw.WriteLine();
				for(int i = 0;i < n;i++)
					v[i] = i + rnd.NextDouble();
				for(int i = 0;i < maxiter;i++)
				{
					r = rayleigh(n,A,v);
					for(int j = 0;j < n;j++)
						A[j,j] -= r;
					QR_Givens(n,n,A,Q);
					copyMatrix(A,R,n,n);
					multiplyMatrix(R,Q,A,n,n,n);
					for(int j = 0;j < n;j++)
						A[j,j] += r;
					for(int j = 0;j < n;j++)
						v[j] = Q[n-1,j];
				}
				printMatrix(n,n,A);
				for(int i = 0;i < n;i++)
				{
					for(int j = 0; j < n; j++)
						sw.Write("{0} ",A[i,j].ToString("N",nfi));
					sw.WriteLine();
				}
				sw.WriteLine();
				Console.WriteLine("Pierwiastki danego rownania wielomianowego to: ");
				sw.WriteLine("Pierwiastki danego rownania wielomianowego to: ");
				int k = 0;
				while(k<n)
				{
					if(k + 1 < n && Math.Abs(A[k+1,k])>eps)
					{
						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 = q - p * p;
						Console.WriteLine("x[{0}]={1}-{2}i",k,p.ToString("N",nfi),Math.Sqrt(d).ToString("N",nfi));
						Console.WriteLine("x[{0}]={1}+{2}i",k+1,p.ToString("N",nfi),Math.Sqrt(d).ToString("N",nfi));
						sw.WriteLine("x[{0}]={1}-{2}i",k,p.ToString("N",nfi),Math.Sqrt(d).ToString("N",nfi));
						sw.WriteLine("x[{0}]={1}+{2}i",k+1,p.ToString("N",nfi),Math.Sqrt(d).ToString("N",nfi));
						k += 2;
					}
					else
					{
						Console.WriteLine("x[{0}]={1}",k,A[k,k].ToString("N",nfi));
						sw.WriteLine("x[{0}]={1}",k,A[k,k].ToString("N",nfi));
						k++;
					}
				}
				esc = (char) Console.ReadKey().Key;
			}
			while(esc != (char)ConsoleKey.Escape);
		}
		}
	}
}