using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace NumMethods
{
    
    class Mtr
    {
        double[,] data;
        int n, k;

        public Mtr(int n, int k)
        {
            this.n = n;
            this.k = k;
            data = new double[n, k];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                    data[i, j] = 0;
        }

        public Mtr(Mtr c)
        {
            this.n = c.n;
            this.k = c.k;
            data = new double[n, k];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                    data[i, j] = c.data[i,j];
        }

        public int getColumns()
        {
            return k;
        }
        public int getRows()
        {
            return n;
        }

        public void makeIdentity()
        {
            if (n != k)
            {
                Console.WriteLine("MAKE IDENTITY ERROR: N != K");
                return;
            }
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                {
                    data[i, j] = ((i == j) ? 1 : 0);
                }        
        }


        public Mtr(double[,] d)
        {
            n = d.GetLength(0);
            k = d.GetLength(1);
            data = new double[n, k];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < n; ++j)
                    data[i, j] = d[i, j];
        }

        public void setValue(double v, int i, int j)
        {
            data[i, j] = v;
        }


        public double getValue(int i, int j)
        {
            return data[i, j];
        }

        static public Mtr operator *(Mtr inp1, Mtr inp2)
        {
            if (inp1.k != inp2.n)
            {
                Console.WriteLine("MATRIX MULT ERROR");
                return new Mtr(1, 1);
            }
            Mtr rs = new Mtr(inp1.n, inp2.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp2.k; ++j)
                {
                    double tmp = 0;
                    for (int c = 0; c < inp1.k; ++c)
                        tmp += inp1.data[i, c] * inp2.data[c, j];
                    rs.setValue(tmp, i, j);
                }
            }

            return rs;
        }

        static public Mtr operator *(Mtr inp1, double alpha)
        {

            Mtr rs = new Mtr(inp1.n, inp1.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp1.k; ++j)
                {
                    rs.setValue(inp1.data[i, j] * alpha, i, j);
                }
            }

            return rs;
        }

       static  public Mtr operator -(Mtr inp1, double alpha)
        {

            Mtr rs = new Mtr(inp1.n, inp1.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp1.k; ++j)
                {
                    rs.setValue(inp1.data[i, j] - alpha, i, j);
                }
            }

            return rs;
        }

       static public Mtr operator -(Mtr inp1, Mtr inp2)
       {
           if (inp1.n != inp2.n || inp1.k != inp2.k)
           {
               Console.WriteLine("Operator- ERROR: Matr size not equal");
               return new Mtr(1, 1);
           }

           Mtr rs = new Mtr(inp1.n, inp1.k);
           for (int i = 0; i < inp1.n; ++i)
           {
               for (int j = 0; j < inp1.k; ++j)
               {
                   rs.setValue(inp1.data[i, j] - inp2.data[i,j], i, j);
               }
           }

           return rs;
       }

       static public Mtr operator +(Mtr inp1, Mtr inp2)
       {
           if (inp1.n != inp2.n || inp1.k != inp2.k)
           {
               Console.WriteLine("Operator+ ERROR: Matr size not equal");
               return new Mtr(1, 1);
           }

           Mtr rs = new Mtr(inp1.n, inp1.k);
           for (int i = 0; i < inp1.n; ++i)
           {
               for (int j = 0; j < inp1.k; ++j)
               {
                   rs.setValue(inp1.data[i, j] + inp2.data[i, j], i, j);
               }
           }

           return rs;
       }

        /// Tensor multiplication of vectors
       static public Mtr operator ^(Mtr inp1, Mtr inp2)
       {
           if (inp1.k != inp2.n || inp1.n != inp2.k)
           {
               Console.WriteLine("Operator^ ERROR: Matr size not equal");
               return new Mtr(1, 1);
           }

           Mtr rs = new Mtr(inp1.n, inp2.k);
           for (int i = 0; i < inp1.n; ++i)
           {
               for (int j = 0; j < inp2.k; ++j)
               {
                   rs.setValue(inp1.data[i, 0] * inp2.data[0, j], i, j);
               }
           }

           return rs;
       }

        public Mtr getTranspozed()
        {
            Mtr ret = new Mtr(k, n);
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < k; ++j)
                {
                    ret.setValue(data[i, j], j, i);
                }
            }
            return ret;
        }

        public void transpoze()
        {
            double[,] ret = new double[k, n];
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < k; ++j)
                {
                    ret[j, i] = data[i, j];
                }
            }
            data = ret;
        }

        public Mtr getSubMatrix(int si, int sj, int ei, int ej)
        {
            if (si >= n || sj >= k || ei >= n || ej >= k ||
                si < 0 || sj < 0 || ei < 0 || ej < 0)
            {
                Console.WriteLine("getSubMatrix RANGE ERROR");
                return new Mtr(1, 1);
            }

            Mtr retMtr = new Mtr(ei - si + 1, ej - sj + 1);
            int fi = 0;
            int fj = 0;
            for (int i = si; i <= ei; ++i)
            {
                fj = 0;
                for (int j = sj; j <= ej; ++j)
                {
                    retMtr.setValue(data[i, j], fi, fj);
                    ++fj;
                }
                ++fi;
            }
            return retMtr;
        }

        public bool isZeroMatrix()
        {
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                    if (data[i, j] != 0)
                        return false;
            return true;
        }

        public void printMatrix()
        {
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < k; ++j)
                {
                   
                    Console.Write("{0, 3:f} ", data[i, j]);
                }
                Console.WriteLine();

            }
            Console.WriteLine(); Console.WriteLine(); Console.WriteLine();
        
        }

        public double getMaximumAboveDiag(ref int ii, ref int jj)
        {
            double val = -1000 * 1000 * 1000;

            for (int i = 0; i < n; ++i)
            {
                for (int j = i + 1; j < k; ++j)
                {
                    if (val < Math.Abs(data[i, j]))
                    {
                        val = Math.Abs(data[i, j]);
                        ii = i;
                        jj = j;
                    }
                }
            }
            return val;
        }
    }






    class Program
    {
     static   int mSign(double vl)
        {
            if (vl >= 0)
                return 1;
            else 
                return -1;
        }

     static bool isNull(Mtr col)
     {
         for (int i = 1; i < col.getRows(); ++i)
             if (col.getValue(i, 0) != 0)
                 return false;

         return true;
     }
     static Mtr uniteMtr(Mtr A, Mtr B)
     {
         int bi = A.getRows() - 1;
        
         for (int i = B.getRows() - 1; i >= 0; --i)
         {
             int bj = A.getColumns() - 1;
             for (int j = B.getColumns() - 1; j >= 0; --j)
             {
                 A.setValue(B.getValue(i, j), bi, bj);
                 --bj;
             }
             --bi;
         }
         return A;
     }
        


     static void planarRotationMethod(Mtr A)
     {
         const double EPS = 0.001;

         A.printMatrix();

         int ii = 0, jj = 0;

         Mtr B = A;
         int n = A.getRows();
         int k = A.getColumns();
         Mtr globalT = new Mtr(n, k);
         globalT.makeIdentity();
         for (int iteration = 0; iteration < 300; ++iteration)
         {
             Mtr locT = new Mtr(n, k);
             locT.makeIdentity();           

             double aii = B.getMaximumAboveDiag(ref ii, ref jj);
             if (Math.Abs(aii) < EPS)
                 break;

             double p = 2 * aii;
             double q = B.getValue(ii, ii) - B.getValue(jj, jj);
             double d = Math.Sqrt(p * p + q * q);
             double r = 0, c = 0, s = 0;
             
             if (Math.Abs(q) > EPS)
             {
                 r = Math.Abs(q) / (2 * d);
                 c = Math.Sqrt(0.5f + r);                 
                 s = Math.Sqrt(0.5f - r) * Math.Sign(p * q);                 
             }
             else 
             {
                 c = s = Math.Sqrt(2.0f) / 2.0;
             }
             locT.setValue(c, ii, ii);
             locT.setValue(c, jj, jj);
             locT.setValue(-s, ii, jj);
             locT.setValue(s, jj, ii);

             globalT = globalT * locT;            
             
             B = locT.getTranspozed() * B * locT;
         }

         Console.WriteLine("Eigen numbers");
         B.printMatrix();
         Console.WriteLine("Eigen vectors");
         globalT.printMatrix();
     }

     static Mtr QR_Decomposition(ref Mtr A)
     {
         int si, sj;
         si = A.getRows();
         sj = A.getColumns();

         Mtr y;
         Mtr U = new Mtr(si, sj);
         U.makeIdentity();
         for (int i = 0; i < A.getColumns() - 1; ++i)
         {
             y = A.getSubMatrix(i, i, A.getRows() - 1, i);
             if (isNull(y))
                 continue;
             Mtr locV = new Mtr(si - i, sj - i);
             Mtr locU = new Mtr(si, sj);
             locU.makeIdentity();
             locV.makeIdentity();

             Mtr omega = new Mtr(si - i, 1);
             Mtr vecx = new Mtr(si - i, 1);
             vecx.setValue(1, 0, 0);
             double yNorm = Math.Sqrt((y.getTranspozed() * y).getValue(0, 0));
             int sgn = mSign(y.getValue(0, 0));

             omega = y + vecx * (sgn * yNorm);
             double omegaNorm = Math.Sqrt((omega.getTranspozed() * omega).getValue(0, 0));
             omega = omega * (1 / omegaNorm);
             locV = locV - (omega ^ omega.getTranspozed()) * 2;


             locU = uniteMtr(locU, locV);

             
             A = locU * A;
             U = locU * U;        
         }
         return U;
     }
        
     static void LuDecomposition(Mtr A, ref Mtr L, ref Mtr U)
     {
         if (A.getRows() != A.getColumns())
             return;
         int n = A.getRows();

         L.makeIdentity();

         for (int i = 0; i < n; ++i)
         {
            
             for (int j = i; j < n; ++j)
             {
                 double temp = A.getValue(i, j);
                 
                 for (int k = 0; k < i; ++k)
                 {
                     temp = temp - U.getValue(k, j) * L.getValue(i, k);
                 }           
                 U.setValue(temp, i, j);
             }

             for (int j = i + 1; j < n; ++j)
             {
                 double temp = A.getValue(j, i);

                 for (int k = 0; k < i; ++k)
                 {
                     temp = temp - U.getValue(k, i) * L.getValue(j, k);
                 }
                 temp = temp / U.getValue(i, i);
                 L.setValue(temp, j, i);
             }         
         }     
     }

     static void inversedLuSolver(Mtr A, Mtr b)
     {
         int n = A.getColumns();

         Mtr L = new Mtr(n, n);
         Mtr U = new Mtr(n, n);
         LuDecomposition(A,ref L, ref U);



         Mtr inversed = new Mtr(n, n);

         for (int j = n - 1; j >= 0; --j)
         {
             double tmp = 1;
             for (int k = j + 1; k < n; ++k)
                 tmp -= U.getValue(j, k) * inversed.getValue(k, j);
             tmp /= U.getValue(j, j);
             inversed.setValue(tmp, j, j);

             for (int i = j - 1; i >= 0; --i)
             {
                 tmp = 0;
                 for (int k = i + 1; k < n; ++k)
                     tmp += U.getValue(i, k) * inversed.getValue(k, j);
                 tmp = -tmp / U.getValue(i, i);
                 inversed.setValue(tmp, i, j);             
             }

             for (int i = j - 1; i >= 0; --i)
             {
                 tmp = 0;
                 for (int k = i + 1; k < n; ++k)
                     tmp += inversed.getValue(j, k) * L.getValue(k, i);
                 inversed.setValue(-tmp, j, i);
             }        
         }

         Mtr answer = inversed * b;

         Console.WriteLine("LU DECOMPOSITION SOLVER: ");
         answer.printMatrix();

         Console.WriteLine("-------------------------------------");
     
     }

     static void QRDecomposeSolver(Mtr A, Mtr b)
     {
      
         int SIZE = b.getRows();
         Mtr U;
         U = QR_Decomposition(ref A);
         b = U * b;
         
         double[] results = new double[SIZE];

         for (int x = SIZE - 1; x >= 0; --x)
         {
             double vl = b.getValue(x, 0);
             for (int j = x + 1; j < SIZE; ++j)
             {
                 vl -= A.getValue(x, j) * results[j];
             }
             double dividor = A.getValue(x, x);
             results[x] = vl / dividor;
         }

         Console.WriteLine("QR decomposition solver:");

         for (int i = 0; i < SIZE; ++i)
         {
             Console.WriteLine("X: {0, 3:f} ", results[i]);
         }
         Console.WriteLine("-------------------------------------");
     
     }

     static int CubeDecompositionSolver(Mtr A, Mtr b)
     {
         int n = A.getRows();
         Mtr U = new Mtr(n, n);
         double uii;
         for (int i = 0; i < n; ++i)
         {
             double tmp = A.getValue(i, i);
             for (int k = 0; k < i; ++k)
                 tmp -= U.getValue(k, i) * U.getValue(k, i);
             if (tmp < 0)
             {
                 Console.WriteLine("Negative tmp in square root!");
                 return -1;
             }
             
             tmp = Math.Sqrt(tmp);
             uii = tmp;
             if (uii == 0)
             {
                 Console.WriteLine("Zero division in uii!");
                 return -1;
             }
             U.setValue(uii, i, i);
             for (int j = i + 1; j < n; ++j)
             {
                 tmp = A.getValue(i, j);
                 for (int k = 0; k < i; ++k)
                     tmp -= U.getValue(k, i) * U.getValue(k, j);
                 tmp /= uii;
                 U.setValue(tmp, i, j);
             }  

         }

         double[] y = new double[n];

         for (int i = 0; i < n; ++i)
         {
             double tmp = b.getValue(i, 0);
             for (int j = 0; j < i; ++j)
             {
                 tmp -= U.getValue(j, i) * y[j];
             }
             tmp /= U.getValue(i, i);
             y[i] = tmp;         
         }

         Mtr ans = new Mtr(n, 1);

         for (int i = n - 1; i >= 0; --i)
         {
             double tmp = y[i];
             for (int k = i + 1; k < n; ++k)
             {
                 tmp -= U.getValue(i, k) * ans.getValue(k, 0);
             }
             tmp /= U.getValue(i, i);
             ans.setValue(tmp, i, 0);         
         }

         Console.WriteLine("Cube decomposition solver:");
         
         ans.printMatrix();

         Console.WriteLine("-------------------------------------");

         return 1;
         

     }

     static void SimpleIteration(Mtr A, Mtr b)
     { 
         int n = A.getRows();
         Mtr B = new Mtr(n, n);
         Mtr C = new Mtr(n, 1);
         Mtr X = new Mtr(n, 1);
         
         for (int i = 0; i < n; ++i)
         {
             for (int j = 0; j < n; ++j)
             {
                 if (i == j)
                     continue;
                 B.setValue(-A.getValue(i, j) / A.getValue(i, i), i, j);             
             }
             C.setValue(b.getValue(i, 0) / A.getValue(i, i), i, 0);
             X.setValue(1, i, 0);
         }

         for (int i = 0; i < 100; ++i)
         {
             X = B * X + C;
         }

         Console.WriteLine("Simple iteration solver:");

         X.printMatrix();


         Console.WriteLine("-------------------------------------");
     
     }


     static double getNorm(Mtr a)
     { 
        
        return Math.Sqrt((a.getTranspozed() * a).getValue(0,0));
     }

     static int pmAlgorithm(Mtr A, ref double bEigen, ref Mtr vecEigen)
     {
         int n = A.getRows();
         const double pmEPS = 0.0000001;

         Mtr y = new Mtr(n, 1);
         double[] diffs = new double[n];

         
         for (int i = 0; i < n; ++i)
         {    y.setValue(1, i, 0);
                diffs[i] = 100;
         }
         double lastSum = 0;
         for (int i = 0; i < 1000; ++i)
         {
             lastSum = 0;
             Mtr x = y * (1 / getNorm(y));
             Mtr ny = A * x;
             bool stop = true;
             for (int j = 0; j < n; ++j)
             {
                 if (Math.Abs(x.getValue(j, 0)) > pmEPS)
                 {
                     if (Math.Abs(ny.getValue(j, 0) / x.getValue(j, 0) - diffs[j]) > pmEPS)
                         stop = false;
                     diffs[j] = ny.getValue(j, 0) / x.getValue(j, 0);
                     lastSum += diffs[j];
                 }
             }
             y = ny;
             if (stop)
                 break;
         }

         
         
         vecEigen = y * (1 / getNorm(y)); ;
         bEigen = lastSum / n;
         return 1;
     }

     static void findMinimalByPm(Mtr A)
     {
         double maxEigen = 0;
         Mtr res = new Mtr(1, 1);
         pmAlgorithm(A, ref maxEigen, ref res);

         Mtr substr = new Mtr(A.getRows(), A.getColumns());
         substr.makeIdentity();
         A = A - substr * maxEigen;
         double nEigen = 0;
         pmAlgorithm(A, ref nEigen, ref res);

         Console.WriteLine(nEigen + maxEigen);
         res.printMatrix();

     
     
     }


     static void Main(string[] args)
        {
            int SIZE = 2;
            Mtr A = new Mtr(SIZE, SIZE);
            Mtr b = new Mtr(SIZE, 1);

            Mtr L = new Mtr(SIZE, SIZE);
            Mtr P = new Mtr(SIZE, SIZE);


            Mtr symmetrical = new Mtr(SIZE, SIZE);
           

            Random rndd = new Random(13335);
            for (int i = 0; i < A.getRows(); ++i)
            {
                for (int j = 0; j < A.getColumns(); ++j)
                {
                    int vl = rndd.Next(100);
                    A.setValue(vl, i, j);
                    if (i == j)
                        A.setValue(vl * 30, i, j);
                    symmetrical.setValue(vl, i, j);
                    symmetrical.setValue(vl, j, i);

                }
                b.setValue(rndd.Next(100), i, 0);
            }
            inversedLuSolver(A, b);

            QRDecomposeSolver(A, b);
            SimpleIteration(A, b);


            symmetrical.setValue(5, 0, 0);
            symmetrical.setValue(1, 0, 1);
           // symmetrical.setValue(2, 0, 2);

            symmetrical.setValue(1, 1, 0);
            symmetrical.setValue(4, 1, 1);
          //  symmetrical.setValue(1, 1, 2);

           // symmetrical.setValue(2, 2, 0);
          //  symmetrical.setValue(1, 2, 1);
          //  symmetrical.setValue(3, 2, 2);
         

            planarRotationMethod(symmetrical);
            findMinimalByPm(symmetrical);
        }
    }
}


/*

 *  Matrices for cube decompose solver
 * 5x5
694,00 690,00 140,00 611,00 466,00
690,00 962,00 529,00 693,00 815,00
140,00 529,00 911,00 243,00 638,00
611,00 693,00 243,00 674,00 527,00
466,00 815,00 638,00 527,00 983,00

 * 4x4
377,00 207,00 262,00 250,00
207,00 883,00 178,00 1,00
262,00 178,00 823,00 358,00
250,00 1,00 358,00 412,00

*/