CDOJ 1330 柱爷与远古法阵（高斯消元）

CDOJ 1330 柱爷与远古法阵（高斯消元）

柱爷与远古法阵

Input

Output

Sample input and output

100 0

33.0476190476

100 2
2 3
99 100

29.8571428571

Hint

Source

1 #include <bits/stdc++.h>
2 using namespace std;
3 const int maxn=305;
4 const long double eps=1e-14;
5 long double a[maxn][maxn];//构造的高斯消元的矩阵，代表第i个方程式的第j个系数是多少 ，精度要求很高
6 int n,m,f[maxn],x,y;
7 inline int read()//读入优化
8 {
9     int x=0,f=1;
10     char ch=getchar();
11     while(ch<'0'||ch>'9')
12     {
13         if(ch=='-')
14             f=-1;
15         ch=getchar();
16     }
17     while(ch>='0'&&ch<='9')
18     {
19         x=x*10+ch-'0';
20         ch=getchar();
21     }
22     return x*f;
23 }
24 inline void write(int x)//输出优化
25 {
26     if(x<0)
27     {
28         putchar('-');
29         x=-x;
30     }
31     if(x>9)
32         write(x/10);
33     putchar(x%10+'0');
34 }
35 int main()
36 {
37     n=read();
38     m=read();
39     for(int i=1;i<=n;i++)
40         f[i]=i;
41     for(int i=1;i<=m;i++)//如果有传送的话，到哪里
42         f[read()]=read();
43      //建立增广矩阵的过程
44     for(int i=1;i<n;i++)
45     {
46         a[i][i]=6;//第一个方程
47         if(f[i]!=i)
48             a[i][f[i]]=-6;//如果有传送门 系数直接抵消 x-y=0 相当于 x=y
49         else
50         {
51             a[i][n+1]=6;//方程右边的常数
52             for(int j=1;j<=6;j++)
53             {
54                 if(i+j<=n)
55                     a[i][i+j]-=1.0;
56                 else
57                     a[i][i]-=1.0;//另外一个方程
58             }
59         }
60     }
61     a

=1.0;//最后的方程
62     a
[n+1]=0;
63     //高斯消元的过程
64     for(int i=1;i<=n;i++)
65     {
66         int p=i;
67         for(int j=i+1;j<=n;j++)
68         {
69             if(fabs(a[j][i])>eps)//向下查找第j个系数不为0的方程
70                 p=j;
71         }
72         if(fabs(a[p][i])>eps)
73         {
74             for(int j=i;j<=n+1;j++)
75                 swap(a[i][j],a[p][j]);//把方程移上来
76             for(int j=i+1;j<=n;j++)//向下消元 同时除去其他的系数
77             {
78                 if(fabs(a[j][i])>eps)
79                 {
80                     long double k=a[j][i]/a[i][i];//消元
81                     for(int t=i;t<=n+1;t++)
82                         a[j][t]-=a[i][t]*k;//系数相减
83                 }
84             }
85         }
86     }
87     //回代过程
88     for(int i=n;i>=1;i--)
89     {
90         for(int j=i+1;j<=n;j++)
91         {
92             if(fabs(a[i][j])>eps)
93                 a[i][n+1]-=a[i][j]*a[j][n+1];//用已知的解求未知解
94         }
95         if(abs(a[i][i])<=eps&&abs(a[i][n+1])>eps)//如果出现矛盾
96         {
97             printf("-1\n");
98             return 0;
99         }
100         a[i][n+1]/=a[i][i];//求出当前的解
101     }
102     printf("%.12lf\n",(double)a[1][n+1]);//a[i][n+1]就是第i个未知数的解
103     return 0;
104 }