using System;
using System.Collections.Generic;
using System.Text;
using System.IO;

namespace NumericRcipes
{
    public sealed class NRVec//向量類定義
    {
        private int n;
        private double[] data;

        public int Rows { get { return n; } }
        public double this[int i] { get { return data[i]; } set { data[i] = value; } }
        public NRVec() { }//構造空向量
        public NRVec(int nn) { n = nn; data = new double[n]; }//構造一個長度為nn的向量
        public NRVec(double a, int nn)//構造一個長度為nn，元素為a的向量
        {
            n = nn; data = new double[n];
            for (int i = 0; i < n; i++) data[i] = a;
        }
        public NRVec(double[] a)//由一維數組構造一個向量
        {
            n = a.GetLength(0);
            data = new double[n];
            a.CopyTo(data, 0);
        }
        public NRVec(NRVec rhs)//復制rhs向量
        {
            n = rhs.n;
            data = new double[n];
            rhs.data.CopyTo(data, 0);
        }
        public override string ToString()//向量的字符串輸出
        {
            StringBuilder sb = new StringBuilder();
            for (int i = 0; i < n; i++)
                sb.AppendFormat("{0:F4}\n", data[i]);
            return sb.ToString();
        }
        public void AssignTo(NRVec vec)//將向量賦值給vec
        {
            vec.n = n;
            vec.data = new double[vec.n];
            data.CopyTo(vec.data, 0);
        }
        public void AssignFrom(double a)//將a賦給向量的每個元素
        {
            for (int i = 0; i < n; i++) data[i] = a;
        }
        public static explicit operator double[](NRVec vec)//顯示轉換向量為一維數組
        {
            double[] result = new double[vec.n];
            vec.data.CopyTo(result, 0);
            return result;
        }
        public static explicit operator NRVec(double[] array)//顯示轉換一維數組為向量
        {
            NRVec result = new NRVec(array);
            return result;
        }
    }
    public sealed class NRMatrix : Object //矩陣類的定義
    {
        #region 矩陣的數據成員
        int rows;                            //記錄矩陣的行數
        int cols;                              //記錄矩陣的列數
        double[,] data;                     //矩陣元素的存儲數組
        public int Rows                 //獲取矩陣的行數
        {
            get { return rows; }
        }
        public int Cols                   //獲取矩陣的列數
        {
            get { return cols; }
        }
        public double this[int i, int j]//獲取和設置矩陣第i行第j列的元素
        {

            get
            {
                if (i >= rows || j >= cols || i < 0 || j < 0)
                    throw (new NRException("行標或列標越界，不能獲取元素!"));
                return data[i, j];
            }
            set
            {
                if (i >= rows || j >= cols || i < 0 || j < 0)
                    throw (new NRException("行標或列標越界，不能設置元素!"));
                data[i, j] = value;
            }
        }
        #endregion
        #region 矩陣的構造函數
        public NRMatrix() { }//構造空矩陣
        public NRMatrix(int rows, int cols)    //構造一個rows×cols的零矩陣
        {
            if (rows <= 0 || cols <= 0) throw (new NRException("不能生成少于1行1列的矩陣!"));
            this.rows = rows; this.cols = cols;
            data = new double[rows, cols];
        }
        public NRMatrix(double[,] a)//二維數組構造矩陣
        {
            rows = a.GetLength(0);
            cols = a.GetLength(1);
            data = new double[rows, cols];
            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    data[i, j] = a[i, j];
        }
        public NRMatrix(NRVec vec)//由向量vec生成對角矩陣
        {
            rows = vec.Rows;
            cols = rows;
            data = new double[rows, cols];
            for (int i = 0; i < rows; i++)
                data[i, i] = vec[i];
        }
        public NRMatrix(NRMatrix a)//復制a矩陣
        {
            rows = a.rows;
            cols = a.cols;
            data = new double[rows, cols];
            for (int i = 0; i < rows; i++)
                for (int j = 0; j < cols; j++)
                    data[i, j] = a[i, j];
        }
        #endregion
        public void AssignTo(NRMatrix mat)//將矩陣賦值給矩陣mat
        {
            mat.rows = rows;
            mat.cols = cols;
            mat.data = new double[mat.rows, mat.cols];
            for (int i = 0; i < mat.rows; i++)
                for (int j = 0; j < mat.cols; j++)
                    mat.data[i, j] = data[i, j];
        }
        public override string ToString()  //矩陣的字符串輸出
        {
            StringBuilder sb = new StringBuilder();
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                    sb.AppendFormat("{0:F4}\t", data[i, j]);
                sb.AppendFormat("\r\n");
            }
            return sb.ToString();
        }
        #region 矩陣基本運算符重載
        public override bool Equals(object o)  //重載Equals,矩陣相等的比較
        {
            // 在處理之前確認調用對象object是矩陣。
            if (o != null && o is NRMatrix)
            {

                if (rows == ((NRMatrix)o).rows && cols == ((NRMatrix)o).cols)
                {
                    for (int i = 0; i < rows; i++)
                        for (int j = 0; j < cols; j++)
                            if (this[i, j] != ((NRMatrix)o)[i, j]) return false;
                    return true;
                }
                else
                    return false;
            }
            return false;
        }
        public static bool operator ==(NRMatrix mat1, NRMatrix mat2)        //矩陣相等的比較
        {
            return mat1.Equals(mat2);
        }
        public override int GetHashCode()
        {
            return this.ToString().GetHashCode();
        }
        public static bool operator !=(NRMatrix mat1, NRMatrix mat2)        //矩陣不等的比較
        {
            return !mat1.Equals(mat2);
        }
        public static NRMatrix operator +(NRMatrix mat1, NRMatrix mat2)     //矩陣的加法運算
        {
            if (mat1.rows != mat2.rows || mat1.cols != mat2.cols)
                throw (new NRException("矩陣不同型，不能作加法！"));
            NRMatrix result = new NRMatrix(mat1.rows, mat1.cols);
            for (int i = 0; i < mat1.rows; i++)
                for (int j = 0; j < mat1.cols; j++)
                    result[i, j] = mat1[i, j] + mat2[i, j];
            return result;
        }
        public static NRMatrix operator -(NRMatrix mat1, NRMatrix mat2)     //矩陣的減法運算
        {
            if (mat1.rows != mat2.rows || mat1.cols != mat2.cols)
                throw (new NRException("矩陣不同型，不能作減法！"));
            NRMatrix result = new NRMatrix(mat1.rows, mat1.cols);
            for (int i = 0; i < mat1.rows; i++)
                for (int j = 0; j < mat1.cols; j++)
                    result[i, j] = mat1[i, j] - mat2[i, j];
            return result;
        }
        public static NRMatrix operator *(NRMatrix mat1, NRMatrix mat2)  //矩陣的乘積運算
        {
            if (mat1.cols != mat2.rows) throw (new NRException("矩陣不相容，不能作矩陣乘法！"));
            int M = mat1.rows, N = mat2.cols;
            NRMatrix result = new NRMatrix(M, N);
            for (int i = 0; i < M; i++)
                for (int j = 0; j < N; j++)
                {
                    double sum = 0.0;
                    for (int k = 0; k < mat1.cols; k++)
                        sum += mat1[i, k] * mat2[k, j];
                    result[i, j] = sum;
                }
            return result;
        }
        public static NRMatrix operator *(double k, NRMatrix mat)  //矩陣的數乘運算
        {
            NRMatrix result = new NRMatrix(mat.rows, mat.cols);
            for (int i = 0; i < mat.rows; i++)
                for (int j = 0; j < mat.cols; j++)
                    result[i, j] = k * mat[i, j];
            return result;
        }