/*================================================================================================*\
| 线性方程组a[][]x[]=b[] : Gauss算法
\*================================================================================================*/
定理: 有方程组AX = b, 只要A非奇异,则可通过逐次消元及行的变换,将方程组化为三角形方程组,且求出唯一解。
#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个相同的开关,每个开关都与某些开关有着联系,每当你打开或者关闭某个开关的时候,其他的与此开关相关联的开关也会相应地发生变化,即这些相联系的开关的状态如果原来为开就变为关,如果为关就变为开。
对于这类问题,巧妙的运用位运算和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
|
/*================================================================================================*\
| Gauss消元算法求解开关灯问题
\*================================================================================================*/
开关问题:有N个相同的开关,每个开关都与某些开关有着联系,每当你打开或者关闭某个开关的时候,其他的与此开关相关联的开关也会相应地发生变化,即这些相联系的开关的状态如果原来为开就变为关,如果为关就变为开。
对于这类问题,巧妙的运用位运算和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
|
|