柱爷与远古法阵

 #include <bits/stdc++.h>

 using namespace std;

 const int maxn=;

 const long double eps=1e-;

 long double a[maxn][maxn];//构造的高斯消元的矩阵，代表第i个方程式的第j个系数是多少 ，精度要求很高

 int n,m,f[maxn],x,y;

 inline int read()//读入优化

 {

     int x=,f=;

     char ch=getchar();

     while(ch<''||ch>'')

     {

         if(ch=='-')

             f=-;

         ch=getchar();

     }

     while(ch>=''&&ch<='')

     {

         x=x*+ch-'';

         ch=getchar();

     }

     return x*f;

 }

 inline void write(int x)//输出优化

 {

     if(x<)

     {

         putchar('-');

         x=-x;

     }

     if(x>)

         write(x/);

     putchar(x%+'');

 }

 int main()

 {

     n=read();

     m=read();

     for(int i=;i<=n;i++)

         f[i]=i;

     for(int i=;i<=m;i++)//如果有传送的话，到哪里

         f[read()]=read();

      //建立增广矩阵的过程

     for(int i=;i<n;i++)

     {

         a[i][i]=;//第一个方程

         if(f[i]!=i)

             a[i][f[i]]=-;//如果有传送门 系数直接抵消 x-y=0 相当于 x=y

         else

         {

             a[i][n+]=;//方程右边的常数

             for(int j=;j<=;j++)

             {

                 if(i+j<=n)

                     a[i][i+j]-=1.0;

                 else

                     a[i][i]-=1.0;//另外一个方程

             }

         }

     }

     a[n][n]=1.0;//最后的方程

     a[n][n+]=;

     //高斯消元的过程

     for(int i=;i<=n;i++)

     {

         int p=i;

         for(int j=i+;j<=n;j++)

         {

             if(fabs(a[j][i])>eps)//向下查找第j个系数不为0的方程

                 p=j;

         }

         if(fabs(a[p][i])>eps)

         {

             for(int j=i;j<=n+;j++)

                 swap(a[i][j],a[p][j]);//把方程移上来

             for(int j=i+;j<=n;j++)//向下消元 同时除去其他的系数

             {

                 if(fabs(a[j][i])>eps)

                 {

                     long double k=a[j][i]/a[i][i];//消元

                     for(int t=i;t<=n+;t++)

                         a[j][t]-=a[i][t]*k;//系数相减

                 }

             }

         }

     }

     //回代过程

     for(int i=n;i>=;i--)

     {

         for(int j=i+;j<=n;j++)

         {

             if(fabs(a[i][j])>eps)

                 a[i][n+]-=a[i][j]*a[j][n+];//用已知的解求未知解

         }

         if(abs(a[i][i])<=eps&&abs(a[i][n+])>eps)//如果出现矛盾

         {

             printf("-1\n");

             return ;

         }

         a[i][n+]/=a[i][i];//求出当前的解

     }

     printf("%.12lf\n",(double)a[][n+]);//a[i][n+1]就是第i个未知数的解

     return ;

 }

这个算法和之前谈及的有点儿不同，它由绝对值最大的部分开始做起，这样可以改善算法上的稳定性。将经过调换后的第一列作为起点，这算法由左至右地计算。每作出以下两个步骤，才跳到下一列：

1.定出每列的最后一个非0的数，将每行的数字除以该数，使到每行的第一个数成为1；

2.将每行的数字减去第一行的第一个数的某个倍数。

所有步骤完成后，这个矩阵会变成一个行梯阵式，再用代入法就可解决这个方程组。

线性方程组其实是相当容易解决的，基本的思想就是消元。但是当未知数较多时，解起来也蛮头疼的。在这里向大家介绍高斯消元法。例如解如下四元一次方程组：

除去各未知数，将各数排在一起，成为矩阵：

为方便起见，用r2+r1表示把第一行各数加到第二行对应数上。

r2-2*r1, r3-3*r1, r4-4*r1,r3+4*r2, r4+3*r2,r4-2*r3, 可得：

化为方程组形式则为：

从而，可得：

X4=4

X3=3

X2=2

X1=1.

说明：对于矩阵采取的变换的合理性，可对照相应方程组的变换加以理解。