【高斯消元】兼 【期望dp】例题

【总览】

高斯消元基本思想是将方程式的系数和常数化为矩阵，通过将矩阵通过行变换成为阶梯状（三角形），然后从小往上逐一求解。

如：$3X_1 + 2X_2 + 1X_3 = 3$

　　$              X_2 + 2X_3 = 1$

　　$2X_1 + X_3 = 0$

化为矩阵为：--->----->----->

然后就可以通过最后一行直接求出$X_3 = ...$，将其带回第二行，算出$X_2$，同理算出$X_1$。

代码很好理解：

inline void gauss(){
    int i, j, k, l;
    for(i = 1; i <= n; i++){
        l = i;
        for(j = i + 1; j <= n; j++)
            if(fabs(matrix[j][i]) > fabs(matrix[l][i])) l = j;
        if(l != i) for(j = i; j <= n + 1; j++)
            swap(matrix[i][j], matrix[l][j]);
        for(j = i + 1; j <= n; j++){
            double tmp = matrix[j][i] / matrix[i][i];
            for(k = i; k <= n + 1; k++)
                matrix[j][k] -= matrix[i][k] * tmp;
        }
    }
    for(i = n; i >= 1; i--){
        double t = matrix[i][n + 1];
        for(j = n; j > i; j--)
            t -= ans[j] * matrix[i][j];
        ans[i] = t / matrix[i][i];
    }
}

 高斯消元最常应用在  期望DP  中。下面是几道例题。

“期望dp讲解及例题”

【BZOJ1013】球形空间产生器sphere

由给出的$n + 1$个坐标，可以列出 $n$个方程，剩下的模板。

【CODE】

#include<iostream>
#include<cstring>
#include<string>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;

const int N = 20;
double matrix[N][N], last[N], t, ans[N];
int n;

inline int read(){
    int i = 0, f = 1; char ch = getchar();
    for(; (ch < '0' || ch > '9') && ch != '-'; ch = getchar());
    if(ch == '-') f = -1, ch = getchar();
    for(; ch >= '0' && ch <= '9'; ch = getchar())
        i = (i << 3) + (i << 1) + (ch -'0');
    return i * f;
}

inline void wr(int x){
    if(x < 0) putchar('-'), x = -x;
    if(x > 9) wr(x / 10);
    putchar(x % 10 + '0');
}

inline void gauss(){
    int i, j, l, k;
    for(i = 1; i <= n; i++){
        l = i;
        for(j = i + 1; j <= n; j++) 
            if(fabs(matrix[j][i]) > fabs(matrix[l][i])) l = j;
        if(l != i) for(j = i; j <= n + 1; j++)
            swap(matrix[i][j], matrix[l][j]);
        for(j = i + 1; j <= n; j++){
            double tmp = matrix[j][i] / matrix[i][i];
            for(k = i; k <= n + 1; k++) 
                matrix[j][k] -= matrix[i][k] * tmp;
        }
    }
    for(i = n; i >= 1; i--){
        double tmp = matrix[i][n + 1];
        for(j = n; j > i; j--)
            tmp -= ans[j] * matrix[i][j];
        ans[i] = tmp / matrix[i][i];
    }
}

int main(){
    n = read();
    for(int i = 1; i <= n; i++) scanf("%lf", &last[i]);
    for(int i = 1; i <= n; i++){
        t = 0;
        for(int j = 1; j <= n; j++){
            double tmp; scanf("%lf", &tmp);
            matrix[i][j] = 2 * (tmp - last[j]);
            t += tmp * tmp - last[j] * last[j];
            last[j] = tmp;
        }
        matrix[i][n + 1] = t;
    }
    gauss();
    for(int i = 1; i <= n; i++){
        if(i < n) printf("%.3lf ", ans[i]);
        else printf("%.3lf
", ans[i]);
    }
    return 0;
}

【BZOJ3143】游走

　　因为要求期望的最小值，那么走的次数多的边肯定要让花费（编号）尽可能小，所以可以先求出从每个点出发次数的期望值$E_i$，那么对于一条边而言，走这条边的期望次数就是$E_i / degree[i] + E_j / degree[j]$，只要排一遍序就好。

　　求点的期望：$E_i = sum (E_{son[i]} / degree[i])$

【code】

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#include<vector>
#define eps 1e-10
using namespace std;

const int N = 600;
double matrix[N][N], ans[N], ret;
int n, m, num, degree[N];
int st[500000], ed[500000];
double gg[500000];

inline void addEdge(const int &u, const int &v){
    degree[u]++;
    degree[v]++;
    st[++num] = u, ed[num] = v;
}

inline int read(){
    int i = 0, f = 1; char ch = getchar();
    for(; (ch < '0' || ch > '9') && ch != '-'; ch = getchar());
    if(ch == '-') f = -1, ch = getchar();
    for(; ch >= '0' && ch <= '9'; ch = getchar())
        i = (i << 3) + (i << 1) + (ch -'0');
    return i * f;
}

inline void gauss(){
    int i, j, l, k;
    for(i = 1; i <= n; i++){
        l = i;
        for(j = i + 1; j <= n; j++) 
            if(fabs(matrix[j][i]) > fabs(matrix[l][i])) l = j;
        if(l != i) for(j = i; j <= n + 1; j++)
            swap(matrix[i][j], matrix[l][j]);
        for(j = i + 1; j <= n; j++){
            double tmp = matrix[j][i] / matrix[i][i];
            for(k = i; k <= n + 1; k++) 
                matrix[j][k] -= matrix[i][k] * tmp;
        }
    }
    for(i = n; i >= 1; i--){
        double tmp = matrix[i][n + 1];
        for(j = n; j > i; j--)
            tmp -= ans[j] * matrix[i][j];
        ans[i] = tmp / matrix[i][i];
    }
}

inline bool cmp (double a, double b){
    return a > b;
}

int main(){
    n = read(), m = read();
    for(int i = 1; i <= m; i++){
        int u = read(), v = read();
        addEdge(u, v);
    }
    int i, j;
    for(i = 1; i <= m; i++){
        matrix[st[i]][ed[i]] += 1.0 /degree[ed[i]];
        matrix[ed[i]][st[i]] += 1.0 /degree[st[i]];
    }
    for(i = 1; i <= n; i++) matrix[n][i] = 0;  
    for(i = 1; i <= n; i++) matrix[i][i] = -1.0;
    matrix[1][n + 1] = -1.0;
    gauss();
    for(i = 1; i <= m; i++)
        gg[i] = ans[st[i]] / degree[st[i]] + ans[ed[i]] / degree[ed[i]];
    sort(gg + 1, gg + m + 1, cmp);
    for(i = 1; i <= m; i++) ret += gg[i] * i;
    printf("%.3f
", ret);
    return 0;
}

 【bzoj2337】XOR和路径

　　学到了！看见求异或和$----->$按位计算：即一位一位的计算答案每一位上为1的期望值，这样就可以轻松统计出答案。

　　每一位都要重新构造矩阵求期望，设$a[i]$表示从$i$到$n$的路径异或和（这一位）为$1$的期望概率（总是≤$1$）

　　对于当前第$i + 1$位，若$(dis >> i) & 1$（这一位为1），那么要异或和为$1$，要求他从关联点异或和为$0$转移来，

　　同理，若这一位为$0$，要求从$1$转移来。即：$$a[i] = sum a[son[i]](dis这一位为0) / degree[i]  + sum a[son[i]](dis这一位为1) / degree[i]$$。

【code】

#include<iostream>
#include<cstring>
#include<string>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cmath>
using namespace std;

const int N = 105, M = 10005;
int n, m;
int ecnt, st[M << 1], ed[M << 1], len[M << 1], degree[N];
double matrix[N][N], ans[N], ret;

inline void addEdge(const int &u, const int &v, const int &l){
    st[++ecnt] = u, ed[ecnt] = v, len[ecnt] = l; degree[u]++;
    if(u != v) st[++ecnt] = v, ed[ecnt] = u, len[ecnt] = l, degree[v]++;
}

inline void gauss(){
    int i, j, k, l;
    for(i = 1; i <= n; i++){
        l = i;
        for(j = i + 1; j <= n; j++)
            if(fabs(matrix[j][i]) > fabs(matrix[l][i])) l = j;
        if(l != i) for(j = i; j <= n + 1; j++)
            swap(matrix[i][j], matrix[l][j]);
        for(j = i + 1; j <= n; j++){
            double tmp = matrix[j][i] / matrix[i][i];
            for(k = i; k <= n + 1; k++)
                matrix[j][k] -= matrix[i][k] * tmp;
        }
    }
    for(i = n; i >= 1; i--){
        double t = matrix[i][n + 1];
        for(j = n; j > i; j--)
            t -= ans[j] * matrix[i][j];
        ans[i] = t / matrix[i][i];
    }
}

int main(){
    scanf("%d%d", &n, &m);
    int i, j, k;
    for(i = 1; i <= m; i++){
        int u, v, w; scanf("%d%d%d", &u, &v, &w);
        addEdge(u, v, w);
    }
    for(i = 0; i <= 30; i++){
        memset(matrix, 0, sizeof matrix);
        memset(ans, 0, sizeof ans);
        for(j = 1; j <= n; j++) matrix[j][j] = 1;
        for(j = 1; j <= ecnt; j++){
            int l = len[j], u = st[j], v = ed[j];
            if(u == n) continue;
            if((l >> i) & 1){
                matrix[u][v] += 1.0 / degree[u];
                matrix[u][n + 1] += 1.0 / degree[u];
            }
            else matrix[u][v] -= 1.0 / degree[u];
        }
        gauss();
        ret += ans[1] * (1 << i);
    }
    printf("%.3f
", ret);
    return 0;
}