/// <summary>/// 行列式计算,本程序属于MyMathLib的一部分。欢迎使用，參考，提意见。/// 有时间用函数语言改写，做自己得MathLib,里面的算法经过验证，但没经过/// 严格測试，如需參考，请谨慎./// </summary>public static partial class LinearAlgebra{ /// <summary>/// 获取指定i,j的余子式/// </summary>/// <param name="Determinants">N阶行列式</param>/// <param name="i">第i行</param>/// <param name="j">第j列</param>/// <returns>计算结果</returns>public static T[,] GetDeterminantMij<T>(T[,] Determinants, int i, int j){var theN = Determinants.GetLength(0);var theNewDeter = new T[theN - 1, theN - 1];int theI = -1;for (int k = 0; k < theN; k++){if (k == i - 1){continue;}theI++;int theJ = -1;for (int l = 0; l < theN; l++){if (l == j - 1){continue;}theJ++;theNewDeter[theI, theJ] = Determinants[k, l];}}return theNewDeter;}/// <summary>/// 获取指定i,j的余子式/// </summary>/// <param name="Determinants">N阶行列式</param>/// <param name="Rows">要取得行</param>/// <param name="Cols">要取得列</param>/// <returns>计算结果</returns>public static T[,] GetDeterminantMij<T>(T[,] Determinants, int[] Rows, int[] Cols){if (Rows.Length != Cols.Length){throw new Exception("所取行数和列数必须相等!");}var theN = Determinants.GetLength(0);var theNewN = theN - Rows.Length;var theNewDeter = new T[theNewN, theNewN];int theI = -1;for (int k = 0; k < theN; k++){if (Rows.Contains(k + 1)){continue;}theI++;int theJ = -1;for (int l = 0; l < theN; l++){if (Cols.Contains(l + 1)){continue;}theJ++;theNewDeter[theI, theJ] = Determinants[k, l];}}return theNewDeter;}/// <summary>/// 获取指定k阶子式N/// </summary>/// <param name="Determinants">N阶行列式</param>/// <param name="Rows">要取得行</param>/// <param name="Cols">要取得列</param>/// <returns>计算结果</returns>public static T[,] GetDeterminantKN<T>(T[,] Determinants, int[] Rows, int[] Cols){if (Rows.Length != Cols.Length){throw new Exception("所取行数和列数必须相等!");}var theNewN = Rows.Length;var theNewDeter = new T[theNewN, theNewN];for (int k = 0; k < Rows.Length; k++){for (int l = 0; l < Cols.Length; l++){theNewDeter[k, l] = Determinants[Rows[k] - 1, Cols[l] - 1];}}return theNewDeter;}/// <summary>/// 计算余子式的符号。

/// </summary> /// <param name="i"></param> /// <param name="j"></param> /// <returns></returns> public static int CalcDeterMijSign(int i, int j) { int theSign = 1; if ((i + j) % 2 == 1) { theSign = -1; } return theSign; } /// <summary> /// 计算余子式的符号。 /// </summary> /// <param name="i"></param> /// <param name="j"></param> /// <returns></returns> public static int CalcDeterMijSign(int[] Rows, int[] Cols) { int theSign = 1; var theSum = Rows.Sum() + Cols.Sum(); if (theSum % 2 == 1) { theSign = -1; } return theSign; } /// <summary> /// 降阶法计算行列式 /// </summary> /// <param name="Determinants">N阶行列式</param> /// <param name="ZeroOptimization">是否0优化</param> /// <returns>计算结果</returns> public static decimal CalcDeterminantAij(decimal[,] Determinants, bool ZeroOptimization = false) { var theN = Determinants.GetLength(0); //假设为2阶，直接计算 if (theN == 2) { return Determinants[0, 0] * Determinants[1, 1] - Determinants[0, 1] * Determinants[1, 0]; } if (ZeroOptimization) { //找0最多的行 int theRowIndex = 0; int theMaxZeroCountR = -1; for (int i = 0; i < theN; i++) { int theZeroNum = 0; for (int j = 0; j < theN; j++) { if (Determinants[i, j] == 0) { theZeroNum++; } } if (theZeroNum > theMaxZeroCountR) { theRowIndex = i; theMaxZeroCountR = theZeroNum; } } //找0最多的列 int theColIndex = 0; int theMaxZeroCountC = -1; for (int i = 0; i < theN; i++) { int theZeroNum = 0; for (int j = 0; j < theN; j++) { if (Determinants[j, i] == 0) { theZeroNum++; } } if (theZeroNum > theMaxZeroCountC) { theColIndex = i; theMaxZeroCountC = theZeroNum; } } if (theMaxZeroCountR >= theMaxZeroCountC) { decimal theRetDec = 0; //第i=theRowIndex+1行展开 int i = theRowIndex + 1; for (int j = 1; j <= theN; j++) { var theSign = CalcDeterMijSign(i, j); var theNewMij = GetDeterminantMij(Determinants, i, j); theRetDec += theSign * Determinants[i - 1, j - 1] * CalcDeterminantAij(theNewMij, ZeroOptimization); } return theRetDec; } else { decimal theRetDec = 0; //第j=theColIndex+1列展开 int j = theColIndex + 1; for (int i = 1; i <= theN; i++) { var theSign = CalcDeterMijSign(i, j); var theNewMij = GetDeterminantMij(Determinants, i, j); theRetDec += theSign * Determinants[i, j] * CalcDeterminantAij(theNewMij, ZeroOptimization); } return theRetDec; } } else { //採用随机法展开一行 var i = new Random().Next(1, theN); decimal theRetDec = 0; for (int j = 1; j <= theN; j++) { var theSign = CalcDeterMijSign(i, j); var theNewMij = GetDeterminantMij(Determinants, i, j); theRetDec += theSign * Determinants[i-1, j-1] * CalcDeterminantAij(theNewMij, ZeroOptimization); } return theRetDec; } } /// <summary> /// 计算范德蒙行列式 /// </summary> /// <param name="Determinants">范德蒙行列式简记序列</param> /// <returns>计算结果</returns> public static decimal CalcVanDerModeDeter(decimal[] VanDerModeDeter) { var theN = VanDerModeDeter.Length; if (theN == 1) { return 1; } decimal theRetDec = 1; for (int i = 0; i < theN; i++) { for (int j = i + 1; j < theN; j++) { theRetDec *= (VanDerModeDeter[j] - VanDerModeDeter[i]); } } return theRetDec; } /// <summary> /// 获取奇数序列 /// </summary> /// <param name="N"></param> /// <returns></returns> private static int[] GetLaplaceRowsOdd(int N) { var theRet = new List<int>(); for (int i = 0; i < N; i = i + 2) { theRet.Add(i + 1); } return theRet.ToArray(); } /// <summary> /// 依据拉普拉斯定理计算行列式值。

/// </summary> /// <param name="Determinants">N阶行列式</param> /// <param name="Rows">初始展开行,里面採用奇数行展开</param> /// <returns>计算结果</returns> public static decimal CalcDeterByLaplaceLaw(decimal[,] Determinants, int[] Rows) { var n = Determinants.GetLength(0); var k = Rows.Length; //假设阶数小于3，则不是必需採用拉普拉斯展开 if (n <= 3) { return CalcDeterminantAij(Determinants, false); } //从P（theN,theK) var theRetList = GetCombination(n, k); decimal theRetDec = 0; foreach (var theCols in theRetList) { var theSign = CalcDeterMijSign(Rows, theCols.ToArray()); var theKN = GetDeterminantKN(Determinants, Rows, theCols.ToArray()); var theN = GetDeterminantMij(Determinants, Rows, theCols.ToArray()); decimal theRetKN = 0; //假设剩余阶数>4则採用随机半数处理. if (n - k >= 4) { var theRows = GetLaplaceRowsOdd(n - k); theRetKN = CalcDeterByLaplaceLaw(theKN, theRows); } else { theRetKN = CalcDeterminantAij(theKN); } decimal theRetAk = 0; if (k >= 4) { var theRows = GetLaplaceRowsOdd(k); theRetAk = CalcDeterByLaplaceLaw(theN, theRows); } else { theRetAk = CalcDeterminantAij(theN); } theRetDec += theSign * theRetKN * theRetAk; } return theRetDec; } /// <summary> /// 从N个数中取k个数的组合结果。考虑到组合数没有顺序区分，因此仅仅要考虑从小 /// 到大的排列下的组合情况就可以，另外，假设组合也不用考虑元素反复的 /// 问题。假设有反复数，仅仅要除重就可以。

/// </summary> /// <param name="N">N个数1-N</param> /// <param name="k">取K个</param> /// <returns></returns> public static List<List<int>> GetCombination(int N, int k) { var theList = new List<int>(); for (int i = 1; i <= N; i++) { theList.Add(i); } return GetCombination(theList, k); } /// <summary> /// 从N个中取k个数,算法原理C(N,k)=C(N-1,k)+ (a + C(Na-1,k-1));当中Na是N中去掉a后的集合. /// </summary> /// <param name="N">元素总个数</param> /// <param name="k">取k个</param> /// <returns></returns> public static List<List<int>> GetCombination(List<int> N, int k) { if (k==0) { return null; } if (N.Count < k) { return null; } if (k == 1) { var theResultsList = new List<List<int>>(); foreach (var theN in N) { var theList = new List<int>(); theList.Add(theN); theResultsList.Add(theList); } return theResultsList; } if (N.Count == k) { var theResultsList = new List<List<int>>(); var theList = new List<int>(); theList.AddRange(N); theResultsList.Add(theList); return theResultsList; } var theRet3 = new List<List<int>>(); int theLeft = N[0]; var theRight = new List<int>(); theRight.AddRange(N); theRight.Remove(N[0]); var theRet2 = GetCombination(theRight, k); theRet3.AddRange(theRet2); theRet2 = GetCombination(theRight, k - 1); for (int n = 0; n < theRet2.Count; n++) { var theList = new List<int>(); theList.Add(theLeft); theList.AddRange(theRet2[n]); theRet3.Add(theList); } return theRet3; } } }