首页 > 其他分享 >高斯消元

高斯消元

时间：2024-07-22 19:51:20浏览次数：6  
标签：int 矩阵 mxpos band inline 高斯消 Mod

高斯消元

高斯消元法通常用于求解如下的 \(n\) 元线性方程组

\[\begin{cases} a_{1, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n, 1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases} \]

高斯消元解线性方程组

考虑求解如下方程组：

\[\begin{cases} x_1 + 3x_2 + 4x_3 = 5 \\ x_1 + 4x_2 + 7x_3 = 3 \\ 9x_1 + 3x_2 + 2x_3 = 2 \end{cases} \]

把系数写到矩阵里：

\[\begin{bmatrix} 1 & 3 & 4 & 5 \\ 1 & 4 & 7 & 3 \\ 9 & 3 & 2 & 2 \\ \end{bmatrix} \]

把未知项的系数单独拿出来考虑：

\[\begin{bmatrix} 1 & 3 & 4 \\ 1 & 4 & 7 \\ 9 & 3 & 2 \\ \end{bmatrix} \]

如果能够将其化为：

\[\begin{bmatrix} ? & ? & ? \\ 0 & ? & ? \\ 0 & 0 & ? \\ \end{bmatrix} \]

就可以从最后一行逆推上去。设第一行第一列为主元，利用加减消元法配系数消掉后两行的第一列：

\[\begin{bmatrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & -24 & -34 \\ \end{bmatrix} \]

再以第二行第二列为主元，加减消元得到：

\[\begin{bmatrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & 0 & 38 \\ \end{bmatrix} \]

即得目标状态，实现时把 \(b\) 的那一列一起带入消元即可。

每次选取系数最大的为主元可以有效降低精度误差。

P3389 【模板】高斯消元法

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-9;
const int N = 1e2 + 7;

double g[N][N], ans[N];

int n;

inline bool Gauss() {
	for (int i = 1; i <= n; ++i) {
		int mxpos = i;

		for (int j = i + 1; j <= n; ++j)
			if (fabs(g[j][i]) > fabs(g[mxpos][i]))
				mxpos = j;

		swap(g[i], g[mxpos]);

		if (fabs(g[i][i]) < eps)
			return false;

		for (int j = i + 1; j <= n; ++j) {
			double div = g[j][i] / g[i][i];

			for (int k = i + 1; k <= n + 1; ++k)
				g[j][k] -= g[i][k] * div;
		}
	}

	return true;
}

signed main() {
	scanf("%d", &n);

	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= n + 1; ++j)
			scanf("%lf", g[i] + j);

	if (!Gauss())
		return puts("No Solution"), 0;

	for (int i = n; i; --i) {
		ans[i] = g[i][n + 1];

		for (int j = i + 1; j <= n; ++j)
			ans[i] -= g[i][j] * ans[j];

		ans[i] /= g[i][i];
	}

	for (int i = 1; i <= n; ++i)
		printf("%.2lf\n", ans[i]);

	return 0;
}

约旦消元法解线性方程组

其为高斯消元的进阶版，能避免回带来求出答案，也就是要把矩阵变为：

\[\begin{bmatrix} ? & 0 & 0 & ? \\ 0 & ? & 0 & ? \\ 0 & 0 & ? & ? \\ \end{bmatrix} \]

首先，我们依照朴素的高斯消元不难得到：

\[\begin{bmatrix} ? & ? & ? & ? \\ 0 & ? & ? & ? \\ 0 & 0 & ? & ? \\ \end{bmatrix} \]

观察上下两个矩阵，不难得到，我们在消元时不仅仅消去主元所在式子下方的式子，而对于上方的式子也应当予以处理。

注意一下特殊情况：

  • 无解：当左边系数为 \(0\) 而右边系数不为 \(0\) 时即为无解。

  • 多解：当左边系数为 \(0\) 而右边系数也为 \(0\) 时即为多解。

inline bool Gauss() {
	for (int i = 1; i <= n; ++i) {
		int mxpos = i;
		
		for (int j = i + 1; j <= n; ++j)
			if (fabs(g[j][i]) > fabs(g[mxpos][i]))
				mxpos = j;
       	
        if (i != mxpos)
            swap(g[i], g[mxpos]);
		
		if (fabs(g[i][i]) < eps)
			return false;
		
		for (int j = 1; j <= n; ++j)
			if (j != i) {
				double div = g[j][i] / g[i][i];
				
				for (int k = i + 1; k <= n + 1; ++k)
					g[j][k] -= g[i][k] * div;
			}
	}
	
	for (int i = 1; i <= n; ++i)
		ans[i] = g[i][n + 1] / g[i][i];
	
	return true;
}

约旦消元法解异或方程组

异或方程组是指形如：

\[\begin{cases} a_{1, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{2, n} x_n = b_2 \\ \cdots \\ a_{n, 1} x_1 \oplus a_{n, 2} x_2 \oplus \cdots \oplus a_{n, n} x_n = b_n \end{cases} \]

的方程组，其中 \(a_{i, j}, b_i\) 的值为 \(0\) 或 \(1\) 。

由于 \(\oplus\) 符合交换律与结合律，故可以按照高斯消元法逐步消元求解。值得注意的是，我们在消元的时候应使用异或消元而非加减消元，且不需要进行乘除改变系数（因为系数均为 \(0\) 和 \(1\) ）。

注意到异或方程组的增广矩阵是 \(01\) 矩阵，可以用 bitset 优化到 \(O(\dfrac{n^2 m}{\omega})\) 。

inline bool Gauss() {
	for (int i = 1; i <= n; ++i) {
		int mxpos = i;

		while (mxpos <= n && !g[mxpos].test(i))
			++mxpos;

		if (mxpos > n)
			return false;
		
        if (i != mxpos)
			swap(g[i], g[mxpos]);

		for (int j = 1; j <= m; ++j)
			if (j != i && g[j].test(i))
				g[j] ^= g[i];
	}

	return true;
}

矩阵求逆

对于矩阵 \(A\) ，若存在矩阵 \(A^{-1}\) 使得 \(A \times A^{-1} = A^{-1} \times A = I\) ，则称矩阵 \(A\) 可逆，\(A^{-1}\) 为其逆矩阵。

给出 \(n\) 阶方阵 \(A\) ，求解其逆矩阵的方法如下：

  • 构造一个 \(n \times 2n\) 的矩阵 \((A, I_n)\)
  • 用高斯消元法将其化简为最简形 \((I_n, A^{-1})\) ，即可得到 \(A\) 的逆矩阵 \(A^{-1}\) 。
  • 如果最终最简形的左半部分不是单位矩阵 \(I_n\) ，则矩阵 \(A\) 不可逆。
inline bool Gauss() {
	for (int i = 1; i <= n; ++i)
		g[i][i + n] = 1;

	for (int i = 1; i <= n; ++i) {
		int mxpos = i;

		while (mxpos <= n && !g[mxpos][i])
			++mxpos;

		if (mxpos > n)
			return false;
	
        if (i != mxpos)
			swap(g[i], g[mxpos]);
        
		int inv = mi(g[i][i], Mod - 2);

		for (int j = 1; j <= n; ++j) {
			if (j == i)
				continue;

			int div = 1ll * g[j][i] * inv % Mod;

			for (int k = i; k <= n * 2; ++k)
				g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
		}
	}

	for (int i = 1; i <= n; ++i) {
		int inv = mi(g[i][i], Mod - 2);

		for (int j = 1; j <= 2 * n; ++j)
			g[i][j] = 1ll * g[i][j] * inv % Mod;
	}

	return true;
}

band-matrix

长这个样子：

空白部分都为 \(0\) ，橙色部分可以为任何数，这样中间就形成了一个宽度大约为 \(d\) 的带。

可以发现任意一个 \(i\) 满足从 \((i, i)\) 向右或向下拓展都有不超过 \(d - 1\) 个非零数字，即很多位置根本不需要消。

具体地，假设现在要消第 \(i\) 列，那么从第 \(i\) 行开始往下枚举 \(d - 1\) 行，每行往右消 \(d\) 个数字即可，最后仍能得到一个上三角矩阵。

与普通高斯消元有点不一样的地方在于当主元为 \(0\) 的时候的处理方法。在 band-matrix 中，若直接交换行会破坏 band-matrix 。注意到每次交换完后交换的行右边最多多出 \(d\) 个数，于是每次往右消元 \(2d\) 个数即可。

时间复杂度 \(O(nd^2)\) 。

inline bool Gauss(int n, int band) {
	for (int i = 1; i <= n; ++i) {
        int mxpos = i;
        
        while (mxpos <= min(i + band, n) && fabs(g[mxpos][i]) < eps)
            ++mxpos;
		
		if (mxpos > min(i + band, n))
			return false;
        
        if (i != mxpos)
	        swap(g[i], g[mxpos]);
			
		for (int j = i + 1; j <= min(i + band, n); ++j) {
			double div = g[j][i] / g[i][i];
			
			for (int k = i; k <= min(i + 2 * band, n); ++k)
				g[j][k] -= div * g[i][k];
			
			g[j][n + 1] -= div * g[i][n + 1];
		}
	}
	
	for (int i = n; i; --i) {
		ans[i] = g[i][n + 1];
		
		for (int j = i + 1; j <= min(i + 2 * band, n); ++j)
			ans[i] -= g[i][j] * ans[j];
		
		ans[i] /= g[i][i];
	}
    
    return true;
}

还有一种解决主元为 \(0\) 的方法，普通高消是交换行，这里只要交换列就可以保持 band-matrix 的性质了，时间复杂度也是 \(O(nd^2)\) 。

inline void Gauss(int n, int band) {
    iota(id + 1, id + 1 + n, 1);
    
    for (int i = 1; i <= n; ++i) {
        if (!g[i][i]) {
            for (int j = i + 1; j <= min(n, i + band); ++j)
                if (g[i][j]) {
                    for (int k = 1; k <= n; ++k)
                        swap(g[k][i], g[k][j]);
                    
                    swap(id[i], id[j]);
                    break;
                }
        }
        
        int Inv = mi(g[i][i], Mod - 2);
        
        for (int j = i + 1; j <= n; ++j) {
            int div = 1ll * g[j][i] * Inv % Mod;
            
            for (int k = i; k <= min(n, i + band); ++k)
                g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
            
            g[j][n + 1] = dec(g[j][n + 1], 1ll * g[i][n + 1] * div % Mod);
        }
    }
    
    for (int i = n; i; --i) {
        ans[id[i]] = g[i][n + 1];
        
        for (int j = i + 1; j <= min(n, i + band); ++j)
            ans[id[i]] = dec(ans[i], 1ll * ans[id[j]] * g[i][j] % Mod);
        
        ans[id[i]] = 1ll * ans[id[i]] * mi(g[i][i], Mod - 2) % Mod;
    }
}

应用

Broken robot

\(n\) 行 \(m\) 列的矩阵，现在在 \((x,y)\)，每次等概率向左、右、下走或原地不动，但不能走出去，求走到最后一行期望的步数。

\(n, m \leq 10^3\)

记 \(f_{i, j}\) 表示机器人在 \((i, j)\) 时走到最后一行的期望步数，则：

\(m = 1\) 时有（省略第二维）：

\[f_i = 1 + \dfrac{1}{2} (f_{i + 1} + f_i) \]

即：

\[f_i = f_{i + 1} + 2 \]

\(m > 1\) 时有：

\[f_{i, j} = 1 + \begin{cases} \dfrac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j + 1}), & j = 1 \\ \dfrac{1}{4} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1}), & 1 < j < m \\ \dfrac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1}), & j = m \end{cases} \]

注意到这是一个 \(d = 2\) 的 band-matrix ，直接使用 band-matrix 消元即可，时间复杂度 \(O(nmd^2)\) 。

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int N = 1e3 + 7;

double g[N][N], ans[N][N];

int n, m, x, y;

inline void Gauss(int n, int band, double *ans) {
	for (int i = 1; i <= n; ++i) {
		int mxpos = i;
        
        while (mxpos <= min(i + band, n) && fabs(g[mxpos][i]) < eps)
            ++mxpos;

        if (mxpos > min(i + band, n))
        	continue;
        
        if (i != mxpos)
        	swap(g[i], g[mxpos]);
			
		for (int j = i + 1; j <= min(i + band, n); ++j) {
			double div = g[j][i] / g[i][i];
			
			for (int k = i; k <= min(i + 2 * band, n); ++k)
				g[j][k] -= div * g[i][k];
			
			g[j][n + 1] -= div * g[i][n + 1];
		}
	}
	
	for (int i = n; i; --i) {
		ans[i] = g[i][n + 1];
		
		for (int j = i + 1; j <= min(i + 2 * band, n); ++j)
			ans[i] -= g[i][j] * ans[j];
		
		ans[i] /= g[i][i];
	}
}

signed main() {
	scanf("%d%d%d%d", &n, &m, &x, &y);
	
	if (m == 1)
		return printf("%.10lf\n", 2.0 * (n - x)), 0;
	
	for (int i = n - 1; i; --i) {
		g[1][1] = 2, g[1][2] = -1, g[1][m + 1] = 3 + ans[i + 1][1];
		
		for (int j = 2; j < m; ++j)
			g[j][j] = 3, g[j][j - 1] = g[j][j + 1] = -1, g[j][m + 1] = 4 + ans[i + 1][j];
		
		g[m][m] = 2, g[m][m - 1] = -1, g[m][m + 1] = 3 + ans[i + 1][m];
		Gauss(m, 1, ans[i]);
	}
	
	printf("%.10lf\n", ans[x][y]);
	return 0;
}

Circles of Waiting

平面直角坐标系上有一个点，一开始在 \((0, 0)\) ，每秒钟这个点都会随机移动，如果它在 \((x, y)\) ，下一秒：

  • 在 \((x - 1, y)\) 的概率是 \(p_1\)
  • 在 \((x, y - 1)\) 的概率是 \(p_2\)
  • 在 \((x + 1, y)\) 的概率是 \(p_3\)
  • 在 \((x, y + 1)\) 的概率是 \(p_4\)

保证 \(p_1 + p_2 + p_3 + p_4 = 1\) ，求该点移动至距离原点距离为大于 \(R\) 的点的期望步数

\(0 \leq R \leq 50\)

把所有满足 \(i^2 + j^2 \leq R^2\) 的点依次编号，显然有 \(O(R^2)\) 个点。

设 \(f_{id(i, j)}\) 表示 \((i, j)\) 走出圆的期望步数，\((i, j)\) 只要转移到 \((i, j - 1), (i - 1, j), (i, j + 1), (i + 1, j)\) 。因为是依次编号，所以建出来的矩阵带宽 \(\leq 2R + 1\) 。套用 band-matrix ，可以做到 \(O(R^2 d^2) = O(R^4)\) 。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 1e9 + 7;
const int R = 5e1 + 7, SIZE = 8e3 + 7;

vector<pair<int, int> > Pos;

int g[SIZE][SIZE];
int id[R << 1][R << 1];
int ans[SIZE];

int r, p1, p2, p3, p4;

template <class T = int>
inline T read() {
    char c = getchar();
    bool sign = c == '-';
    
    while (c < '0' || c > '9')
        c = getchar(), sign |= c == '-';
    
    T x = 0;
    
    while ('0' <= c && c <= '9')
        x = (x << 1) + (x << 3) + (c & 15), c = getchar();
    
    return sign ? (~x + 1) : x;
}

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int g, int b) {
    int res = 1;
    
    for (; b; b >>= 1, g = 1ll * g * g % Mod)
        if (b & 1)
            res = 1ll * res * g % Mod;
    
    return res;
}

inline int getid(int x, int y) {
    return id[x + R][y + R];
}

inline void Gauss(int n, int band) {
    for (int i = 1; i <= n; ++i) {
        int mxpos = i;
        
        while (mxpos <= min(i + band, n) && !g[mxpos][i])
            ++mxpos;

        if (mxpos > min(i + band, n))
            continue;
        
        if (i != mxpos)
            swap(g[i], g[mxpos]);
        
        int inv = mi(g[i][i], Mod - 2);
            
        for (int j = i + 1; j <= min(i + band, n); ++j) {
            int div = 1ll * g[j][i] * inv % Mod;
            
            for (int k = i; k <= min(i + 2 * band, n); ++k)
                g[j][k] = dec(g[j][k], 1ll * div * g[i][k] % Mod);
            
            g[j][n + 1] = dec(g[j][n + 1], 1ll * div * g[i][n + 1] % Mod);
        }
    }
    
    for (int i = n; i; --i) {
        ans[i] = g[i][n + 1];
        
        for (int j = i + 1; j <= min(i + 2 * band, n); ++j)
            ans[i] = dec(ans[i], 1ll * g[i][j] * ans[j] % Mod);
        
        ans[i] = 1ll * ans[i] * mi(g[i][i], Mod - 2) % Mod;
    }
}

signed main() {
    r = read(), p1 = read(), p2 = read(), p3 = read(), p4 = read();
    int invsum = mi(add(add(p1, p2), add(p3, p4)), Mod - 2);
    p1 = 1ll * p1 * invsum % Mod, p2 = 1ll * p2 * invsum % Mod;
    p3 = 1ll * p3 * invsum % Mod, p4 = 1ll * p4 * invsum % Mod;
    
    for (int i = -r; i <= r; ++i)
        for (int j = -r; j <= r; ++j)
            if (i * i + j * j <= r * r)
                Pos.emplace_back(make_pair(i, j)), id[i + R][j + R] = Pos.size();
    
    int n = Pos.size(), band = 0;
    
    for (auto it : Pos) {
        int x = it.first, y = it.second;
        g[getid(x, y)][getid(x, y)] = g[getid(x, y)][n + 1] = 1;
        
        if (getid(x - 1, y)) {
            g[getid(x, y)][getid(x - 1, y)] = Mod - p1;
            band = max(band, abs(getid(x, y) - getid(x - 1, y)));
        }
        
        if (getid(x, y - 1)) {
            g[getid(x, y)][getid(x, y - 1)] = Mod - p2;
            band = max(band, abs(getid(x, y) - getid(x, y - 1)));
        }
        
        if (getid(x + 1, y)) {
            g[getid(x, y)][getid(x + 1, y)] = Mod - p3;
            band = max(band, abs(getid(x, y) - getid(x + 1, y)));
        }
        
        if (getid(x, y + 1)) {
            g[getid(x, y)][getid(x, y + 1)] = Mod - p4;
            band = max(band, abs(getid(x, y) - getid(x, y + 1)));
        }
    }
    
    Gauss(n, band);
    printf("%d", ans[getid(0, 0)]);
    return 0;
}

P4457 [BJOI2018] 治疗之雨

你有 \(p\) 滴血量，血量上限为 \(n\) 。每轮操作如下：

  • 先以 \(\dfrac{1}{m + 1}\) 的概率增加 \(1\) 滴血，满血时则概率为 \(0\) 。
  • \(k\) 次判定，每次以 \(\dfrac{1}{m + 1}\) 的概率减少一滴血，死了则概率为 \(0\) 。

求期望几轮死亡。

\(n \leq 1500\) ，\(m, k \leq 10^9\)

设 \(f_i\) 表示血量为 \(i\) 时期望多少轮结束，则：

\[f_i = 1 + \sum_{j = 0}^{\min(i, k)} (\dfrac{m}{m + 1} f_{i - j} + \dfrac{1}{m + 1} f_{i - j + 1}) \times \dbinom{k}{j} (\dfrac{1}{m + 1})^j (\dfrac{m}{m + 1})^{k - j} \]

注意一下 \(i = n\) 时的情况即可。

观察到 \(f_i\) 只与 \(f_{0 \sim i + 1}\) 有关，所以矩阵应该是一个类似下三角矩阵的东西。因为我们高斯消元的时候是拿自己这行去减下面的，所以每一行中只有 \(2\) 个系数要去和下面的相减，时间复杂度就可以做到 \(O(n^2)\) 了。

#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1.5e3 + 7;

int g[N][N];
int inv[N], d[N], id[N], ans[N];

int n, p, m, k;

template <class T = int>
inline T read() {
    char c = getchar();
    bool sign = (c == '-');
    
    while (c < '0' || c > '9')
        c = getchar(), sign |= (c == '-');
    
    T x = 0;
    
    while ('0' <= c && c <= '9')
        x = (x << 1) + (x << 3) + (c & 15), c = getchar();
    
    return sign ? (~x + 1) : x;
}

inline int add(int x, int y) {
    x += y;
    
    if (x >= Mod)
        x -= Mod;
    
    return x;
}

inline int dec(int x, int y) {
    x -= y;
    
    if (x < 0)
        x += Mod;
    
    return x;
}

inline int mi(int a, int b) {
    int res = 1;
    
    for (; b; b >>= 1, a = 1ll * a * a % Mod)
        if (b & 1)
            res = 1ll * res * a % Mod;
    
    return res;
}

inline void prework() {
    inv[0] = inv[1] = 1;
    
    for (int i = 2; i < N; ++i) 
        inv[i] = 1ll * (Mod - Mod / i) * inv[Mod % i] % Mod;
}

inline void Gauss(int n, int band) {
    iota(id + 1, id + 1 + n, 1);
    
    for (int i = 1; i <= n; ++i) {
        if (!g[i][i]) {
            for (int j = i + 1; j <= min(n, i + band); ++j)
                if (g[i][j]) {
                    for (int k = 1; k <= n; ++k)
                        swap(g[k][i], g[k][j]);
                    
                    swap(id[i], id[j]);
                    break;
                }
        }
        
        int Inv = mi(g[i][i], Mod - 2);
        
        for (int j = i + 1; j <= n; ++j) {
            int div = 1ll * g[j][i] * Inv % Mod;
            
            for (int k = i; k <= min(n, i + band); ++k)
                g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
            
            g[j][n + 1] = dec(g[j][n + 1], 1ll * g[i][n + 1] * div % Mod);
        }
    }
    
    for (int i = n; i; --i) {
        ans[id[i]] = g[i][n + 1];
        
        for (int j = i + 1; j <= min(n, i + band); ++j)
            ans[id[i]] = dec(ans[i], 1ll * ans[id[j]] * g[i][j] % Mod);
        
        ans[id[i]] = 1ll * ans[id[i]] * mi(g[i][i], Mod - 2) % Mod;
    }
}

signed main() {
    prework();
    int T = read();
    
    while (T--) {
        n = read(), p = read(), m = read(), k = read();
        
        if (!k) {
            puts("-1");
            continue;
        } else if (!m) {
            if (k == 1)
                puts("-1");
            else
                printf("%d\n", (p + k - 2 - (p == n)) / (k - 1));
            
            continue;
        }
        
        d[0] = 1ll * mi(m, k) * mi(mi(m + 1, k), Mod - 2) % Mod;
        
        for (int i = 1, invm = mi(m, Mod - 2); i <= min(n, k); ++i)
            d[i] = 1ll * d[i - 1] * invm % Mod * (k - i + 1) % Mod * inv[i] % Mod;
        
        memset(g, 0, sizeof(g));
        
        for (int i = 1, invm1 = mi(m + 1, Mod - 2); i < n; ++i) {
            g[i][n + 1] = g[i][i] = 1;
            
            for (int j = 0; j <= min(i, k); ++j) {
                g[i][i - j] = dec(g[i][i - j], 1ll * dec(1, invm1) * d[j] % Mod);
                g[i][i - j + 1] = dec(g[i][i - j + 1], 1ll * invm1 * d[j] % Mod);
            }
        }

        g[n][n + 1] = g[n][n] = 1;

        for (int i = 0; i <= min(n, k); ++i)
        	g[n][n - i] = dec(g[n][n - i], d[i]);
        
        Gauss(n, 1);
        printf("%d\n", ans[p]);
    }
    
    return 0;
}

P6899 [ICPC2014 WF] Pachinko

有一个宽度为 \(w\) 高度为 \(h\) 的方格纸， $ w \times h$ 的格子中，有一些是空的，有一些是洞，有一些是障碍物。从第一行的空的格子中随机选一个放置一个球，向上下左右移动的概率比为 \(p_u : p_d : p_l : p_r\) ，不能移动到有障碍物的格子上。对于每个洞，输出落入该洞的概率。

\(2 \leq 20, h \leq 10^4, p_u + p_d + p_l + p_r = 100\) 。

和上面比较类似。

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int dx[] = {-1, 1, 0, 0};
const int dy[] = {0, 0, -1, 1};
const int N = 1e4 + 7, M = 2e1 + 7;

double g[N * M][M << 1], ans[N * M];
int id[N][M];
char str[N][M];

int n, m, p[4], tot;

inline bool check(int x, int y) {
	return x && x <= n && y && y <= m;
}

inline void Gauss(int n, int band) {
	for (int i = 1; i <= n; ++i) {
		if (fabs(g[i][i]) < eps)
			continue;
		
		for (int j = i + 1; j <= min(i + band, n); ++j)
			g[i][j] /= g[i][i];

		ans[i] /= g[i][i], g[i][i] = 1;

		for (int j = i + 1; j <= min(i + band, n); ++j) {
			double div = g[j][i] / g[i][i];

			for (int k = i; k <= min(i + band, n); ++k)
				g[j][k] -= g[i][k] * div;

			ans[j] -= ans[i] * div;
		}
	}

	for (int i = n; i; --i)
		for (int j = i + 1; j <= min(i + band, n); ++j)
			ans[i] -= g[i][j] * ans[j];
}

signed main() {
	scanf("%d%d", &m, &n);

	for (int i = 0; i < 4; ++i)
		scanf("%d", p + i);

	for (int i = 1; i <= n; ++i) {
		scanf("%s", str[i] + 1);

		for (int j = 1; j <= m; ++j)
			if (str[i][j] != 'X')
				id[i][j] = ++tot;
	}

	int sum = m - count(id[1] + 1, id[1] + 1 + m, 0);

	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= m; ++j) {
			if (!id[i][j])
				continue;

			if (i == 1)
				ans[id[i][j]] = 1.0 / sum;

			if (str[i][j] == 'T')
				continue;

			g[id[i][j]][id[i][j]] = 1;
			int base = 100;

			for (int k = 0; k < 4; ++k) {
				int x = i + dx[k], y = j + dy[k];

				if (!check(x, y) || !id[x][y])
					base -= p[k];
			}

			for (int k = 0; k < 4; ++k) {
				int x = i + dx[k], y = j + dy[k];

				if (check(x, y) && id[x][y])
					g[id[x][y]][id[i][j]] = -1.0 * p[k] / base;
			}
		}

	Gauss(tot, m);

	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= m; ++j)
			if (str[i][j] == 'T')
				printf("%.9lf\n", ans[id[i][j]]);

	return 0;
}

标签：int,矩阵,mxpos,band,inline,高斯消,Mod
From： https://www.cnblogs.com/wshcl/p/18316770/Gauss

相关文章

  • 高斯消元
    #include<iostream>#include<cassert>#include<cstdio>#include<cstring>#include<algorithm>#include<vector>#include<map>#include<cmath>#include<queue>#include<set>#include<cli......
  • 【数学】高斯消元
    1.算法简介高斯消元法（Gauss–Jordanelimination）是求解线性方程组的经典算法。例如求解下列方程组：\(\begin{cases}2x+9y-5z=10\\4x+20y+z=24\\x-2y+3z=8\end{cases}\)形式化的，高斯消元可用于求解类似于\(\begin{cases}a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_n=b......
  • G68 实数线性基+高斯消元法 P3265 [JLOI2015] 装备购买
    视频链接：G68实数线性基+高斯消元法P3265[JLOI2015]装备购买_哔哩哔哩_bilibili  P3265[JLOI2015]装备购买-洛谷|计算机科学教育新生态(luogu.com.cn)//线性基+高斯消元法O(n*m*m)#include<iostream>#include<cstring>#include<algorithm>usingnames......
  • 高斯消元
    高斯-约旦消元解线性方程组例题：线性方程组步骤:选出未被更新的行中第\(k\)列绝对值最大的值，令为主元把主元所在行移到当前行，加减消元消去主元重复12结束后若存在找不到主元（找到是0）的情况，那就遍历没处理的行，如果有常数项非\(0\)则无解，否则无数解点击查看代......
  • 高斯消元
    前言：由于作者未系统过学习线性代数，故下文肯定有不严谨的成分，不过应付算法竞赛是绰绰有余了QAQ。\[\begin{cases}a_{1,1}x_1+a_{1,2}x_2+\cdots+a_{1,n}x_n=b_1\\a_{2,1}x_1+a_{2,2}x_2+\cdots+a_{2,n}x_n=b_2\\\cdots\\a_{n,1}x_1+a_{n,2......
  • 高斯消元和矩阵快速幂
    高斯消元高斯消元是一种能在\(O(N^3)\)的时间内求解\(N\)元一次方程组的算法。其思路大致如下：使第一个未知数只有第一个式子中系数非\(0\)。使第二个未知数只有第二个式子中系数非\(0\)。\(\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\vdots\)使第......
  • AcWing算法基础课笔记——高斯消元
    高斯消元用来求解方程组a11x1+......
  • 高斯消元学习笔记
    引入高斯-约当消元法（Gauss–Jordanelimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。高斯消元法除了用于线性方程组求解外，还可以用于行列式计算、求矩阵的逆，以及其他计算机和工程方面。过程一个经典的问题，给定一......
  • 高斯消元学习笔记
    高斯消元学习笔记其实这个主题能够复活主要还是粘了\(\text{LGV}\)引理的光，不然我还不知道高斯消元其实不光能求解线性方程组。求解线性方程组这个只能说是典中典了，我不相信没有一个人的高斯消元不是从这里开始的。我们考虑求解线性方程组的本质：将每一个式子所有未知数前都......
  • 高斯消元学习笔记——P304题解
    如果你觉得这篇太啰嗦问题[SDOI2006]线性方程组题目描述已知\(n\)元线性一次方程组。\[\begin{cases}a_{1,1}x_1+a_{1,2}x_2+\cdots+a_{1,n}x_n=b_1\\a_{2,1}x_1+a_{2,2}x_2+\cdots+a_{2,n}x_n=b_2\\\cdots\\a_{n,1}x_1+a_{n,2}x......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99