高斯消元 模板 hdu2262 判断无穷解 接受无穷解 两种

有几个人找我要高斯消元模板，还提出不同要求，一劳永逸，发一个。

最初自己只用 不接受有某个x无穷多解的模板，后来有人来问我，于是自己弄了一个。

//模板1：每个x都有解则解出所有x，并返回1，如果有x无穷解，返回0（不解方程）。
bool gauss1(int n)
{
    int i,j,k,r;
    for(i=0;i<n;++i)
    {
        //数据稳定性优化
        r = i;
        for(j=i+1;j<n;++j) if(fabs(A[j][i]>fabs(A[r][i]))) r=j;
        //if(fabs(a[r][i]) < eps) continue;
        if(r!=i)
        {
            for(j=0;j<n;++j)swap(A[r][j],A[i][j]);
            swap(B[r],B[i]);
        }
        //消元
        for(j=n;j>=i;j--)
        for(k=i+1;k<n;++k)
            if(j==n) B[k]-=A[k][i]/A[i][i]*B[i];
            else A[k][j]-=A[k][i]/A[i][i]*A[i][j];
    }
    //判断是唯一解,
    for(i=0;i<n;++i)if(fabs(A[i][i])<eps) return 0;
    //回代
    for(i=n-1;i>=0;--i)
{
//if(fabs(a[i][i]) < eps) continue;
        for(j=i+1;j<n;++j)B[i]-=A[i][j]*B[j];
        B[i]/=A[i][i];
    }
    return 1;
}

//模板2：方程中有某列全是0，即该x无穷多解，而剩余x仍有唯一解，则解出能解的x
bool gauss2(int n)
{
    int i,j,k,r;
    for(i=0;i<n;++i)
    {
        //数据稳定性优化
        r = i;
        for(j=i+1;j<n;++j) if(fabs(A[j][i]>fabs(A[r][i]))) r=j;
if(fabs(A[r][i]) < eps) continue;
        if(r!=i)
        {
            for(j=0;j<n;++j)swap(A[r][j],A[i][j]);
            swap(B[r],B[i]);
        }
        //消元
        for(j=n;j>=i;j--)
        for(k=i+1;k<n;++k)
            if(j==n) B[k]-=A[k][i]/A[i][i]*B[i];
            else A[k][j]-=A[k][i]/A[i][i]*A[i][j];
    }
  //  for(i=0;i<n;++i)if(fabs(A[i][i])<eps) return 0;
    //回代
    for(i=n-1;i>=0;--i)
{
        if(fabs(A[i][i]) < eps) continue;
        for(j=i+1;j<n;++j)B[i]-=A[i][j]*B[j];
        B[i]/=A[i][i];
    }
    return 1;
}

建议是用模板1

hdu2262，根据两种模板，分别对应两种代码，一种需要构造方程式去掉无穷解的方程（通过bfs实现），另一种则把所有点都拿来构造方程

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#define MS(x,y) memset(x,y,sizeof(x))
using namespace std;
const int MAXN = 0x3f3f3f3f;
const int N = 1010;
const double eps = 1e-10;
void fre(){freopen("t.txt","r",stdin);}
typedef long long LL;
int dir[4][2] = {1,0,-1,0,0,-1,0,1};
char mp[20][20];
bool e[20][20];
int sx,sy,tot,cnt[20][20],vis[20][20],n,m;
double A[405][405];
bool bfs()
{
    queue<int>q;
    int u,i,j,x,y,xx,yy;
    tot = 0;
    vis[sx][sy] = tot++;
    q.push(sx*m+sy);
    while(!q.empty())
    {
        u = q.front(); q.pop();
        x = u/m; y = u%m;
        for(i = 0; i < 4; ++i)
        {
            xx = x + dir[i][0];
            yy = y + dir[i][1];
            if(xx < 0 || yy < 0 || xx >= n || yy >= m || mp[xx][yy] == '#') continue;
            cnt[x][y]++;
            if(vis[xx][yy]!=-1) continue;
            vis[xx][yy] = tot++;
            q.push(xx*m+yy);
        }
    }
    //judge
    for(i = 0; i < n; ++i)
        for(j = 0; j < m; ++j)
        if(vis[i][j]!=-1 && e[i][j]) return 1;

    return 0;
}
bool gauss(int n)
{
    int i,j,k,r;
    for(i=0;i<n;++i)
    {
        //数据稳定性优化
        r = i;
        for(j=i+1;j<n;++j) if(fabs(A[j][i]>fabs(A[r][i]))) r=j;
        if(r!=i) for(j=0;j<=n;++j)swap(A[r][j],A[i][j]);
        //消元
        for(j=n;j>=i;j--)
        for(k=i+1;k<n;++k) A[k][j]-=A[k][i]/A[i][i]*A[i][j];
    }
    //判断是唯一解
    for(i=0;i<n;++i)if(fabs(A[i][i])<eps) return 0;
    //回代
    for(i=n-1;i>=0;--i)
    {
        for(j=i+1;j<n;++j)A[i][n]-=A[i][j]*A[j][n];
        A[i][n]/=A[i][i];
    }
    return 1;
}
int main()
{
    //fre();
    double p;
    int i,j,u,xx,yy,k;
    while(~scanf("%d%d",&n,&m))
    {
        //init
        MS(e,0); MS(cnt,0); MS(vis,-1);
        //read
        for(i = 0; i < n; ++i)
        {
            scanf("%s",mp[i]);
            for(j = 0; j < m; ++j)
            {
                if(mp[i][j]=='@') {sx=i;sy=j;}
                if(mp[i][j]=='$') e[i][j] = 1;
            }
        }
        //bfs
        if(!bfs()) {printf("-1\n");continue;}
        //make matrix
        MS(A,0);
        for(i = 0; i < n; ++i)
        {
            for(j = 0; j < m; ++j)
            {
                if(vis[i][j]==-1)continue;
                u = vis[i][j];
                A[u][u] = 1;
                if(e[i][j]) A[u][tot] = 0;
                else
                {
                    A[u][tot] = 1;
                    p = 1.0/cnt[i][j];//wrong!
                    for(k = 0; k < 4; ++k)
                    {
                        xx = i + dir[k][0]; yy = j + dir[k][1];
                        if(xx < 0 || yy < 0 || xx >= n || yy >= m || mp[xx][yy] == '#') continue;
                        A[u][vis[xx][yy]] = -p;
                    }
                }
            }
        }
        //gauss
        gauss(tot);
        printf("%.6lf\n",A[vis[sx][sy]][tot]);
    }
}

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<queue>
#include<algorithm>
#define eps 1e-18
using namespace std;
struct Node{
int x;
int y;
};
const int N = 810;
const int dirx[4] = {0, 0, 1, -1};
const int diry[4] = {-1, 1, 0, 0};
int equ, var, col, max_r, x, y, cnt, sx, sy, ex, ey, n, m, u;
double tmp;
double a[N][N], ans[N], b[N];
int flag[N][N], t[N];
string ch[N];
queue<Node> q;
struct Node v, p;
void Gauss(int n)
{
    int i,j,k,r;
    for(i=0;i<n;++i)
    {
        //数据稳定性优化
        r = i;
        for(j=i+1;j<n;++j) if(fabs(a[j][i]>fabs(a[r][i]))) r=j;
        if(fabs(a[r][i]) < eps) continue;
        if(r!=i)
        {
            for(j=0;j<n;++j)swap(a[r][j],a[i][j]);
            swap(b[r],b[i]);
        }
        //消元
        for(j=n;j>=i;j--)
        for(k=i+1;k<n;++k)
            if(j==n) b[k]-=a[k][i]/a[i][i]*b[i];
            else a[k][j]-=a[k][i]/a[i][i]*a[i][j];
    }
    //判断是唯一解
    //for(i=0;i<n;++i)if(fabs(a[i][i])<eps) return 0;
    //回代
    for(i=n-1;i>=0;--i)
    {
        if(fabs(a[i][i]) < eps) continue;
        for(j=i+1;j<n;++j)b[i]-=a[i][j]*b[j];
        b[i]/=a[i][i];
    }

   // return 1;
}


int bfs()
{
    memset(flag, 0, sizeof(flag));
    while (!q.empty()) q.pop();
    p.x = sx;
    p.y = sy;
    int t = 0;
    q.push(p);
    flag[sx][sy] = 1;
    while (!q.empty())
    {
        v = q.front();
        q.pop();
        for (int i = 0; i < 4; i++)
        {
            x = v.x + dirx[i];
            y = v.y + diry[i];
            if (x >= 0 && x < n && y >= 0 && y < m && !flag[x][y] && ch[x][y] != '#')
            {
                p.x = x;
                p.y = y;
                if (ch[x][y] == '$') t = 1;
                flag[x][y] = 1;
                q.push(p);
            }
        }
    }
    return t;
}
void work()
{
    memset(a, 0, sizeof(a));memset(b, 0, sizeof(a));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
        if (flag[i][j] == 1)
        {
            cnt = 0;

            if (ch[i][j] == '$')
            {
                a[i*m+j][i*m+j] = 1;
                continue;
            }
            for (int k = 0; k < 4; k++)
            {
                x = i + dirx[k];
                y = j + diry[k];
                if (x >= 0 && x < n && y >= 0 && y < m && flag[x][y])
                {
                    cnt++;
                    a[i*m+j][x*m+y] = 1;
                }
            }
            b[i*m+j] = -1*cnt;
            a[i*m+j][i*m+j] = -1*cnt;
        }
}
int main()
{
    //freopen("t.txt","r",stdin);
    while (~scanf("%d%d", &n, &m))
    {
        memset(t, 0, sizeof(t));
        for (int i = 0; i < n; i++)
        {
            cin >> ch[i];
            for (int j = 0; j < m; j++)
            {
                if (ch[i][j] == '@')
                {
                    sx = i;
                    sy = j;
                }

            }
        }
        if (!bfs()) printf("-1\n");
        else
        {
            work();
            Gauss(n*m);
            printf("%.6f\n", b[sx*m+sy]);
        }
    }
    return 0;
}