線性方程組a[][]x[]=b[] ： Gauss算法

/*================================================================================================*\

|                                   線性方程組a[][]x[]=b[] ： Gauss算法   

\*================================================================================================*/

#define eps 1e-10                                // 定義精度

#define fabs(x) ((x)>0?(x):-(x))                 // 取絕對值

#define zero(x) (fabs(x)<eps)                  // 判定是否為0

const int maxn = 100;                            // 未知數的個數

double a[maxn][maxn+1]; int n;                   // 增廣矩陣

void pmat() {                                    // 輸出增廣矩陣，默認為浮點型

      for (int i = 0, j; i < n; i++) {

             for (j = 0; j < n; j++)

                     printf("%lf ", a[i][j]);

             printf("%lf\n", a[i][j]);

      }

}

----------------------------------------------------------------------------------------

列主元消元法： 僅在一列當中選取絕對值最大的元素作為消去的主元素的Gauss算法。增廣矩陣表示。

----------------------------------------------------------------------------------------

double getMaxRow(int k, int &row) {            // 輸出矩陣第a[k..n-1][k]中最大元素和行

      double ret = 0;

      for (int i = k; i < n; i++)

             if (fabs(a[i][k]) > fabs(ret))

                     ret = a[row=i][k];

      return ret;

}

void swapRow(int k, int i) {                   // 交第k行和第i行： a[k] 和 a[i];

      for (int j = k; j <= n; j++)  

             swap(a[k][j], a[i][j]);

}

int gauss() {

      int i, j, k, row;

      double maxp;      // pmat();

      for (k = 0; k < n; k++) {

             maxp = getMaxRow(k, row);

             if (zero(maxp)) return 0;                 // 如果為奇異陣，則有無數解。

             if (row != k) swapRow(k, row);             // 需要交換兩行

             for (j = k + 1; j  <= n; j++) {            // 加減法消元。

                     a[k][j] /= maxp;

                     for (i = k+1; i < n; i++)

                             a[i][j] -= a[i][k]*a[k][j];

             }

      }//pmat();

      for (i = n-1; i >= 0; i--)                       // 結果保留在a[0..n-1][n]中

             for (j = i+1; j < n; j++)

                     a[i][n] -= a[i][j]*a[j][n];

      return 1;

}// gauss

----------------------------------------------------------------------------------------

全主元消元法： 選取剩余的最大元素作為消去的主元素，交換行與列的Gauss算法。增廣矩陣表示。

----------------------------------------------------------------------------------------

double getMax(int k, int &row, int &col) {          // 在a[k..n][k..n]中查找最大元素

      double ret = 0;                               // 并保留最大元所在的行與列

      for (int i = k; i < n; i++)

             for (int j = k; j < n; j++)

                     if (fabs(a[i][j]) > fabs(ret))

                             ret = a[row=i][col=i];

      return ret;

}

void swapCol(int k, int j) {                       // 交換第k列和第j列

      for (int i = 0; i < n; i++)

             swap(a[i][j], a[i][k]);

}

void swapRow(int k, int i) {                       // 交換第k行和第i行

      for (int j = k; j <= n; j++)

             swap(a[k][j], a[i][j]);

}

int index[maxn]; double t[maxn];                  // 兩個輔助數組。

bool gauss() {

      int i, j, k, row, col;

      double maxp;

      for (i = 0; i < n; i++) index[i] = i;

      for (k = 0; k < n; k++) {

             maxp = getMax(k, row, col);

             if (zero(maxp)) return 0;

             if (col != k) {

                     swapCol(k, col);

                     swap(index[col], index[k]); // 這里要交換索引。

             }

             if (row != k) swapRow(k, row);

             for (j = k+1; j <= n; j++) {

                     a[k][j] /= maxp;

                     for (i = k+1; i < n; i++)

                             a[i][j] -= a[i][k]*a[k][j];

             }

      }

      for (i = n-1; i >= 0; i--)

             for (j = i+1; j < n; j++)

                     a[i][n] -= a[i][j]*a[j][n];

      for (k = 0; k < n; k++) t[index[k]] = a[k][n];

      for (k = 0; k < n; k++) a[k][n] = t[k];

      return 1;

}

 

比較： 全主元消元法 和 列主元消元法比較起來，其精度更高，如果要求精度較高的情況，選擇全主元消元法，精度較低時，選擇列主元消元法。

----------------------------------------------------------------------------------------

高斯-諾爾當消元法：除主對角線元素為1之外其余元素都為零的消元法。增廣矩陣表示。

----------------------------------------------------------------------------------------

double getRow(int k, int &row) {            // 獲得當前列的a[k..n][k]的第一個非零元

      for (int i = k; i < n; i++)

             if (!zero(a[i][k])) return a[row=i][k];

      return 0;

}

void swapRow(int k, int i) {               // 交行第k行和第i行

      for (int j = k; j <= n; j++)

             swap(a[k][j], a[i][j]);

}

bool gauss() {

      int i, j, k, row;

      double ret;

      for (k = 0; k < n; k++) {

             ret = getRow(k, row);//pmat();

             if (zero(ret)) return 0;

             if (row != k) swapRow(k, row);

             for (j = k; j <= n; j++)

                     a[k][j] /= ret;

             for (j = k+1; j <= n; j++) {

                     for (i = 0; i < k; i++)

                             a[i][j] -= a[i][k]*a[k][j];

                     for (i = k+1; i < n; i++)

                             a[i][j] -= a[i][k]*a[k][j];

             }

      }

      return 1;

}

 

 

 

 

 

 

 

/*================================================================================================*\

|                                   Gauss消元算法求解開關燈問題   

\*================================================================================================*/

----------------------------------------------------------------------------------------

開燈問題告訴N(N<=63)盞燈和M(M<=N)個開關,每個開關可以控制K(K<=N)盞燈，給定N盞燈的初始狀態S和要求通過開關控制得到的目標狀態E，求可以達到目標狀態的方案數。

----------------------------------------------------------------------------------------

簡單分析：對于每個開關，有按和不按兩種選擇（記為0/1）; 對于每盞燈有變和不變兩種情況（0/1）,如果初態和終態不一樣，那么這盞燈是一定要變化的。由此我們就可以得到一個0/1的矩陣：讓N盞燈燈作為列向量，開關作為橫向量，把每盞燈是否變化作為第M列（由0開始）這樣就得到一個N*(M+1)的矩陣,該矩陣有如下性質：

1. 如果N = M ，那么矩陣為增廣矩陣。

2. 該矩陣相當于方程組A * X = B,因此可以求其解。

   1. 若方程組有唯一解，那么，N = M (逆命題：如果M = N ,那么方程組有唯一解 不成立)

   2. 若方程組無實數解，那么，該方程不可以化成嚴格上三角形式（具體的證明見相關資料，這里不再證明）

   3. 若方程組有多接，即存在自由變元，因為每個自由變元可以取0/1兩種情況，那么總共有2^m(m為變元數)解。

下面是經過驗證的代碼：

int getRow(int p, int q, int &row) {

      for (int i = p; i < n; i++)

             if (!zero(a[i][q])) return a[row=i][q];

      return row=0;

}

void swapRow(int p, int row, int q) {

      for (int k = q; k <= m; k++)

             swap(a[p][k], a[row][k]);

}

i64 gauss() {

      int i = -1, j = -1, k, p, q, ret, row;

      while(++i < n && ++j < m) {

             ret = getRow(i, j, row);

             if (zero(ret)) { i--; continue;}

             if (row != i) swapRow(i, row, j);

             for (p = i+1; p < n; p++) if (a[p][j])

                     for (q = j; q <= m; q++)

                             a[p][q] ^= a[i][q];

      }

      for (k = i; k < n; k++) if(a[k][m]) return -1;

      return (i64)1 << (m-i);

}    //link: hdu3364 http://acm.hdu.edu.cn/showproblem.php?pid=3364

----------------------------------------------------------------------------------------

開關問題：有N個相同的開關，每個開關都與某些開關有著聯系，每當你打開或者關閉某個開關的時候，其他的與此開關相關聯的開關也會相應地發生變化，即這些相聯系的開關的狀態如果原來為開就變為關，如果為關就變為開。

求： 1. 方案數（自由變元的數目）  2. 給定一個最少的開關方案

----------------------------------------------------------------------------------------

這類問題是上面問題的一種簡化：

對于問題一、可以直接套用上面的公式（N=M）

對于問題二、如果構造得到的方程組只有一個解，那么問題解決，這里主要討論一下多解，存在自由變元的情況。如果存在自由變元，我們就要枚舉每個自由變元（0/1）然后比較選擇最小。枚舉時間復雜度為2^m(m為自由變元的個數)

下面是簡單的枚舉自由元的算法。

int gans(int a[][maxn+1]) {

      int i, j, ret = a[n-1][n];

      for (i = n-2; i >= 0; i--) {

             for (j = i+1; j < n; j++)

                     a[i][n] ^= a[i][j] && a[j][n];

             ret += b[i][n];

      }

      return ret;

}

void dfs(int p, int k) {

      if (p == k) {

             memcpy(b, a, sizeof(b));

             int ret = gans(b);

             if (ret < ans) ans = ret;

             return;

      }

      a[p][n] = 1; dfs(p-1, k);

      a[p][n] = 0; dfs(p-1, k);

}

int gauss() { //……代碼見上（n=m）……//

      dfs(n-1, i-1);

      return ans;

}

Link: pku_1222 http://162.105.81.212/JudgeOnline/problem?id=1222

      pku_1681 http://162.105.81.212/JudgeOnline/problem?id=1681

pku_1753 http://162.105.81.212/JudgeOnline/problem?id=1753

pku_1830 http://162.105.81.212/JudgeOnline/problem?id=1830

      pku_3185 http://162.105.81.212/JudgeOnline/problem?id=3185

 

 

 

 

posted on 2010-08-06 22:57 小志 閱讀(411) 評論(0)  編輯 收藏 引用