高斯消元法是用来解线性方程组的一种算法

板子

for(reg i = 1;i < n;i++)
    {
        int r = i;
        while(!a[r][i]) r++;
        swap(a[r],a[i]);
        double pass = a[i][i];
        for(reg k = 1;k <= n;k++) a[i][k] /= pass;
        for(reg j = 1;j < n;j++)
        {
            if(j == i) continue;
            double div = a[j][i];
            for(reg k = 1;k <= n;k++)
                a[j][k] -= a[i][k] * div;
        }
    }

用途

这里主要说下用途

做很多关于期望的题时

\[dp[i]=f(dp[j])\]

\[dp[j]=f(dp[i])\]

出现了以上形式的\(dp\)方程

明显直接\(dp\)是不行的

但如果 我们能求出

\(n\)个关于\(dp[i]=f(dp[j])(i,j\in[1,n],s.t.i,j\in Z)\)

就可以高斯消元来求这些未知量

例题

P3232 [HNOI2013]游走

#include <stack>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define reg register int
#define isdigit(x) ('0' <= x&&x <= '9')
template<typename T>
inline T Read(T Type)
{
    T x = 0,f = 1;
    char a = getchar();
    while(!isdigit(a)) {if(a == '-') f = -1;a = getchar();}
    while(isdigit(a)) {x = (x << 1) + (x << 3) + (a ^ '0');a = getchar();}
    return x * f;
}
const int MAXN = 510,MAXM = MAXN * MAXN;
struct node
{
    int v,next;
}edge[MAXM << 1];
double ans,a[MAXN][MAXN],p[MAXN],edge_p[MAXM];
int cnt,n,m,head[MAXN],deg[MAXN];
inline void addedge(int u,int v)
{
    deg[u]++;
    edge[++cnt].v = v;
    edge[cnt].next = head[u];
    head[u] = cnt;
}
inline void gauss()
{
    a[1][n] = -1;
    for(reg i = 1;i < n;i++)
    {
        for(reg e = head[i];e;e = edge[e].next)
        {
            int v = edge[e].v;
            if(v == n) continue;
            a[i][v] += 1.0 / deg[v];
        }
        a[i][i] -= 1;
    }
    for(reg i = 1;i < n;i++)
    {
        int r = i;
        while(!a[r][i]) r++;
        swap(a[r],a[i]);
        double pass = a[i][i];
        for(reg k = 1;k <= n;k++) a[i][k] /= pass;
        for(reg j = 1;j < n;j++)
        {
            if(j == i) continue;
            double div = a[j][i];
            for(reg k = 1;k <= n;k++)
                a[j][k] -= a[i][k] * div;
        }
    }
    for(reg i = 1;i < n;i++) p[i] = a[i][n] / a[i][i];
    cnt = 0;
    for(reg i = 1;i <= n;i++)
        for(reg e = head[i];e;e = edge[e].next)
        {
            int v = edge[e].v;
            if(v < i) continue;
            edge_p[++cnt] = (i == n?0:p[i] / deg[i]) + (v == n?0:p[v] / deg[v]);
        }
}
int main()
{
    n = Read(1),m = Read(1);
    for(reg i = 1;i <= m;i++)
    {
        int u = Read(1),v = Read(1);
        addedge(u,v),addedge(v,u);
    }
    gauss();
    sort(edge_p + 1,edge_p + 1 + cnt);
    for(reg i = 1;i <= m;i++) ans += edge_p[i] * (m - i + 1);
    printf("%.3f",ans);
    return 0;
}

P3211 [HNOI2011]XOR和路径

#include <stack>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define reg register int
#define isdigit(x) ('0' <= x&&x <= '9')
template<typename T>
inline T Read(T Type)
{
    T x = 0,f = 1;
    char a = getchar();
    while(!isdigit(a)) {if(a == '-') f = -1;a = getchar();}
    while(isdigit(a)) {x = (x << 1) + (x << 3) + (a ^ '0');a = getchar();}
    return x * f;
}
const int MAXN = 110,MAXM = 10010;
struct node
{
    int v,w,next;
}edge[MAXM << 1];
double ans,a[MAXN][MAXN];
int n,m,maxt,cnt,head[MAXN],deg[MAXN];
inline void addedge(int u,int v,int w)
{
    edge[++cnt].v = v;
    edge[cnt].w = w;
    edge[cnt].next = head[u];
    head[u] = cnt;
}
inline int getlen(int x) {int res = 0;while(x) x >>= 1,res++;return res;}
inline int getx(int x,int t) {return (x >> (t - 1) & 1);}
inline double fabs(double x) {if(x < 0) return -x;return x;}
inline void gauss(int t)
{
    for(reg i = 1;i <= n;i++)
        for(reg j = 1;j <= n + 1;j++)
            a[i][j] = 0;
    for(reg i = 1;i < n;i++)
    {
        for(reg j = head[i];j;j = edge[j].next)
        {
            int v = edge[j].v,w = getx(edge[j].w,t);
            if(w) a[i][v] -= 1.0,a[i][n + 1] -= 1.0;
            else a[i][v] += 1;
        }
        a[i][i] -= deg[i];
    }
    a[n][n] = 1;
    for(reg i = 1;i <= n;i++)
    {
        int r = i;
        while(!a[r][i]) r++;
        swap(a[r],a[i]);
        double pass = a[i][i];
        for(reg j = 1;j <= n + 1;j++) a[i][j] /= pass;
        for(reg j = 1;j <= n;j++)
        {
            if(j == i) continue;
            double pass = a[j][i];
            for(reg k = 1;k <= n + 1;k++)
                a[j][k] -= a[i][k] * pass;
        }
    }
    ans += a[1][n + 1] / a[1][1] * (1 << (t - 1));
}
int main()
{
    n = Read(1),m = Read(1);
    for(reg i = 1;i <= m;i++)
    {
        int u = Read(1),v = Read(1),w = Read(1);
        addedge(u,v,w),deg[u]++;
        maxt = max(maxt,getlen(w));
        if(u != v) addedge(v,u,w),deg[v]++;
    }
    for(reg i = 1;i <= maxt;i++) gauss(i);
    printf("%.3f",ans);
    return 0;
}

P4035 [JSOI2008]球形空间产生器

#include <stack>
#include <cstdio>
#include <vector>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define reg register int
#define isdigit(x) ('0' <= x&&x <= '9')
template<typename T>
inline T Read(T Type)
{
    T x = 0,f = 1;
    char a = getchar();
    while(!isdigit(a)) {if(a == '-') f = -1;a = getchar();}
    while(isdigit(a)) {x = (x << 1) + (x << 3) + (a ^ '0');a = getchar();}
    return x * f;
}
const int MAXN = 15;
double pow[MAXN],a[MAXN][MAXN];
inline double fabs(double x) {if(x < 0) return -x;return x;}
struct Matrix
{
    int n,m;
    double c[MAXN][MAXN];
}T;
int main()
{
    int n = Read(1);
    for(reg i = 1;i <= n + 1;i++)
    {
        for(reg j = 1;j <= n;j++)
        {
            scanf("%lf",&a[i][j]);
            pow[i] += a[i][j] * a[i][j];
            if(i != 1) T.c[i - 1][j] = 2 * (a[i - 1][j] - a[i][j]);
        }
        if(i != 1) T.c[i - 1][n + 1] = pow[i - 1] - pow[i];
    }
    for(reg i = 1;i <= n;i++)
    {
        int r = i;
        for(reg j = i + 1;j <= n;j++)
            if(fabs(T.c[r][i]) > fabs(T.c[j][i])) r = j;
        swap(T.c[r],T.c[i]);
        double pass = T.c[i][i];
        for(reg j = 1;j <= n + 1;j++) T.c[i][j] /= pass;
        for(reg j = 1;j <= n;j++)
        {
            if(i == j) continue;
            double pass = T.c[j][i];
            for(reg k = 1;k <= n + 1;k++) T.c[j][k] -= T.c[i][k] * pass;
        }
    }
    for(reg i = 1;i <= n;i++)
        printf("%.3f ",T.c[i][n + 1]);
    return 0;
}
01-01 06:44