MyMathLib系列(行列式计算2)

/// <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, j] * 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;         }     } }