//
// Copyright (c) 2002
// Ronald Kevin Burton
//
// Z poniższym kodem nie jest związana żadna gwarancja poprawności działania.
// Program został dołączony do książki ".NET CLR. Księga eksperta" w celu
// ilustracji koncepcji i zasad przedstawionych w tej książce. Program może być 
// używany na własne ryzyko.
//
// Przyznaje się prawo do użycia lub kopiowania tego oprogramowania do dowolnego celu
// bez konieczności ponoszenia żadnych opłat pod warunkiem, że powyższe uwagi zostaną 
// zachowane we wszystkich kopiach. Przyznaje się także prawo do modyfikacji kodu
// i dystrybucji zmodyfikowanego kodu pod warunkiem zachowania powyższych uwag
// oraz dołączenia informacji mówiącej o modyfikacji kodu.
//
// 
using System;
using System.Diagnostics;
using System.Runtime.InteropServices;

namespace Matrix
{
	/// <summary>
	/// Podsumowanie dla IMatrix.
	/// </summary>
	public interface IMatrix
	{
		void Add(object a, object b, ref object c);
		void Subtract(object a, object b, ref object c);
		void Multiply(object a, object b, ref object c);
		void Inverse(object a, ref object b);
	}
	[ClassInterface(ClassInterfaceType.None)]
	public class MatrixOperations : IMatrix
	{
		public void Add(object a, object b, ref object c)
		{
			object[,] oa = (object [,])a;
			object[,] ob = (object [,])b;
			object[,] oc = (object [,])c;
			if(oa.GetLength(0) != ob.GetLength(0) ||
			   oa.GetLength(0) != oc.GetLength(0) ||
			   ob.GetLength(0) != oc.GetLength(0) ||
			   oa.GetLength(1) != ob.GetLength(1) ||
			   oa.GetLength(1) != oc.GetLength(1) ||
			   ob.GetLength(1) != oc.GetLength(1))
			{
				throw new ArgumentException("Tablice muszą mieć taki sam kształt, aby możliwe było ich dodanie", "A,B,C");
			}
			if(oa[1,1].GetType() == typeof(double))
			{
				double [,] fa = new double [oa.GetLength(0), oa.GetLength(1)];
				double [,] fb = new double [ob.GetLength(0), ob.GetLength(1)];
				double [,] fc = new double [oc.GetLength(0), oc.GetLength(1)];
				for(int i = 0; i < oa.GetLength(0); i++)
					for(int j = 0; j < oa.GetLength(1); j++)
					{
						fa[i,j] = (double)oa[i+1,j+1];
					}
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
					{
						fb[i,j] = (double)ob[i+1,j+1];
					}
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
					{
						fc[i,j] = 0.0;
					}
				_Add(fa, fb, ref fc);
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
						oc[i+1,j+1] = fc[i,j];
				return;
			}
			if(oa[1,1].GetType() == typeof(int))
			{
				int [,] ia = new int [oa.GetLength(0), oa.GetLength(1)];
				int [,] ib = new int [ob.GetLength(0), ob.GetLength(1)];
				int [,] ic = new int [oc.GetLength(0), oc.GetLength(1)];
				for(int i = 0; i < oa.GetLength(0); i++)
					for(int j = 0; j < oa.GetLength(1); j++)
					{
						ia[i,j] = (int)oa[i+1,j+1];
					}
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
					{
						ib[i,j] = (int)ob[i+1,j+1];
					}
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
					{
						ic[i,j] = 0;
					}
				_Add(ia, ib, ref ic);
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
						oc[i+1,j+1] = ic[i,j];
				return;
			}
		}
		public void Subtract(object a, object b, ref object c)
		{
			object[,] oa = (object [,])a;
			object[,] ob = (object [,])b;
			object[,] oc = (object [,])c;
			if(oa.GetLength(0) != ob.GetLength(0) ||
			   oa.GetLength(0) != oc.GetLength(0) ||
			   ob.GetLength(0) != oc.GetLength(0) ||
			   oa.GetLength(1) != ob.GetLength(1) ||
			   oa.GetLength(1) != oc.GetLength(1) ||
			   ob.GetLength(1) != oc.GetLength(1))
			{
				throw new ArgumentException("Tablice muszą mieć taki sam kształt, aby możliwe było ich odjęcie", "A,B,C");
			}
			if(oa[1,1].GetType() == typeof(double))
			{
				double [,] fa = new double [oa.GetLength(0), oa.GetLength(1)];
				double [,] fb = new double [ob.GetLength(0), ob.GetLength(1)];
				double [,] fc = new double [oc.GetLength(0), oc.GetLength(1)];
				for(int i = 0; i < oa.GetLength(0); i++)
					for(int j = 0; j < oa.GetLength(1); j++)
					{
						fa[i,j] = (double)oa[i+1,j+1];
					}
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
					{
						fb[i,j] = (double)ob[i+1,j+1];
					}
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
					{
						fc[i,j] = 0.0;
					}
				_Subtract(fa, fb, ref fc);
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
						oc[i+1,j+1] = fc[i,j];
				return;
			}
			if(oa[1,1].GetType() == typeof(int))
			{
				int [,] ia = new int [oa.GetLength(0), oa.GetLength(1)];
				int [,] ib = new int [ob.GetLength(0), ob.GetLength(1)];
				int [,] ic = new int [oc.GetLength(0), oc.GetLength(1)];
				for(int i = 0; i < oa.GetLength(0); i++)
					for(int j = 0; j < oa.GetLength(1); j++)
					{
						ia[i,j] = (int)oa[i+1,j+1];
					}
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
					{
						ib[i,j] = (int)ob[i+1,j+1];
					}
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
					{
						ic[i,j] = 0;
					}
				_Subtract(ia, ib, ref ic);
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
						oc[i+1,j+1] = ic[i,j];
				return;
			}
		}
		public void Multiply(object a, object b, ref object c)
		{
			object[,] oa = (object [,])a;
			object[,] ob = (object [,])b;
			object[,] oc = (object [,])c;

			if(oa.GetLength(1) != ob.GetLength(0) ||
			   oc.GetLength(1) != ob.GetLength(1) ||
			   oa.GetLength(0) != oc.GetLength(0))
			{
				throw new ArgumentException("kolumny A = wiersze B kolumny C = kolumny B wiersze A = wiersze C", "A,B,C");
			}

			if(oa[1,1].GetType() == typeof(double))
			{
				double [,] fa = new double [oa.GetLength(0), oa.GetLength(1)];
				double [,] fb = new double [ob.GetLength(0), ob.GetLength(1)];
				double [,] fc = new double [oc.GetLength(0), oc.GetLength(1)];
				for(int i = 0; i < oa.GetLength(0); i++)
					for(int j = 0; j < oa.GetLength(1); j++)
					{
						fa[i,j] = (double)oa[i+1,j+1];
					}
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
					{
						fb[i,j] = (double)ob[i+1,j+1];
					}
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
					{
						fc[i,j] = 0.0;
					}
				_Multiply(fa, fb, ref fc);
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
						oc[i+1,j+1] = fc[i,j];
				return;
			}
			if(oa[1,1].GetType() == typeof(int))
			{
				int [,] ia = new int [oa.GetLength(0), oa.GetLength(1)];
				int [,] ib = new int [ob.GetLength(0), ob.GetLength(1)];
				int [,] ic = new int [oc.GetLength(0), oc.GetLength(1)];
				for(int i = 0; i < oa.GetLength(0); i++)
					for(int j = 0; j < oa.GetLength(1); j++)
					{
						ia[i,j] = (int)oa[i+1,j+1];
					}
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
					{
						ib[i,j] = (int)ob[i+1,j+1];
					}
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
					{
						ic[i,j] = 0;
					}
				_Multiply(ia, ib, ref ic);
				for(int i = 0; i < oc.GetLength(0); i++)
					for(int j = 0; j < oc.GetLength(1); j++)
						oc[i+1,j+1] = ic[i,j];
				return;
			}
		}

		public void Inverse(object a, ref object b)
		{
			object[,] oa = (object [,])a;
			object[,] ob = (object [,])b;
			if(oa.GetLength(0) != oa.GetLength(1) ||
				oa.GetLength(0) != ob.GetLength(0) ||
				ob.GetLength(0) != ob.GetLength(1))
			{
				throw new ArgumentException("Ta procedura może odwrócić tylko macierz kwadratową", "A");
			}

			if(oa[1,1].GetType() == typeof(double))
			{
				// Ustalono, iż jest to 
				// macierz N x N.
				int n = oa.GetLength(0);
				// Utworzenie macierzy rozszerzonej
				double [,] fa = new double [n, n * 2];
				double [,] x = new double [ n, n];
				int [] nrow = new int[n];
				for(int i = 0; i < n; i++)
				{
					nrow[i] = i;
					for(int j = 0; j < n; j++)
					{
						fa[i,j] = (double)oa[i+1,j+1];
					}
					for(int j = n; j < n * 2; j++)
					{
						if(j == i + n)
							fa[i,j] = 1.0;
						else
							fa[i,j] = 0.0;
					}
				}
				// Pivot
				int p;
				int ncopy;
				double max;
				double mij;
				double sum;
				for(int i = 0; i < n - 1; i++)
				{
					p = i;
					max = fa[nrow[p],i];
					for(int j = i; j < n; j++)
					{
						if(fa[nrow[j],i] > max)
						{
							max = fa[nrow[j],i];
							p = j;
						}
					}
					if(fa[nrow[p],i] == 0.0)
						throw new ArgumentException("Brak jednoznacznego rozwiązania", "A");
					// Symulowana zamiana wierszy
					if(nrow[i] != nrow[p])
					{
						ncopy = nrow[i];
						nrow[i] = nrow[p];
						nrow[p] = ncopy;
					}
					for(int j = i + 1; j < n; j++)
					{
						mij = fa[nrow[j],i]/fa[nrow[i],i];
						for(int k = 0; k < n * 2; k++)
						{
							fa[nrow[j],k] = fa[nrow[j],k] - mij*fa[nrow[i],k];
						}
					}
					if(fa[nrow[n-1],n-1] == 0.0)
						throw new ArgumentException("Brak jednoznacznego rozwiązania", "A");
				}
				// Rozpoczęcie podstawiania
				for(int k = 0; k < n; k++)
				{
					x[n-1,k] = fa[nrow[n-1],n + k]/fa[nrow[n-1],n-1];
					for(int i = n - 2; i >= 0; i--)
					{
						sum = 0.0;
						for(int j = i + 1;j < n;j++)
							sum += fa[nrow[i],j] * x[j,k];
						x[i, k] = (fa[nrow[i], n + k] - sum)/fa[nrow[i],i];
					}
				}
				// Przekazanie wyników
				for(int i = 0; i < ob.GetLength(0); i++)
					for(int j = 0; j < ob.GetLength(1); j++)
						ob[i+1,j+1] = x[i,j];
			}
		}
		private void _Add(double [,] a, double [,] b, ref double [,] c)
		{
			// Pętla dla wierszy A
			for(int i = 0; i < a.GetLength(0); i++)
				// Pętla dla kolumn A i wierszy B
				for(int j = 0; j < a.GetLength(1); j++)
				{
					c[i,j] = a[i,j] + b[i,j];
				}
		}
		private void _Add(int [,] a, int [,] b, ref int [,] c)
		{
			// Pętla dla wierszy A
			for(int i = 0; i < a.GetLength(0); i++)
				// Pętla dla kolumn A i wierszy B
				for(int j = 0; j < a.GetLength(1); j++)
				{
					c[i,j] = a[i,j] + b[i,j];
				}
		}
		private void _Subtract(double [,] a, double [,] b, ref double [,] c)
		{
			// Pętla dla wierszy A
			for(int i = 0; i < a.GetLength(0); i++)
				// Pętla dla kolumn A i wierszy B
				for(int j = 0; j < a.GetLength(1); j++)
				{
					c[i,j] = a[i,j] - b[i,j];
				}
		}
		private void _Subtract(int [,] a, int [,] b, ref int [,] c)
		{
			// Pętla dla wierszy A
			for(int i = 0; i < a.GetLength(0); i++)
				// Pętla dla kolumn A i wierszy B
				for(int j = 0; j < a.GetLength(1); j++)
				{
					c[i,j] = a[i,j] - b[i,j];
				}
		}
		private void _Multiply(double [,] a, double [,] b, ref double [,] c)
		{
			// Pętla dla wierszy A
			for(int i = 0; i < a.GetLength(0); i++)
				// Pętla dla kolumn A i wierszy B
				for(int j = 0; j < b.GetLength(1); j++)
				{
					// Pętla dla kolumn B
					for(int k = 0; k < b.GetLength(1); k++)
						c[i,k] += a[i,j] * b[j,k];
				}
		}
		private void _Multiply(int [,] a, int [,] b, ref int [,] c)
		{
			// Pętla dla wierszy A
			for(int i = 0; i < a.GetLength(0); i++)
				// Pętla dla kolumn A i wierszy B
				for(int j = 0; j < b.GetLength(1); j++)
				{
					// Pętla dla kolumn B
					for(int k = 0; k < b.GetLength(1); k++)
						c[i,k] += a[i,j] * b[j,k];
				}
		}
	}
}
