wenbao与高斯消元

消元

高斯消元

 1 typedef double Matrix[maxn][maxn];
 2 void gauss_elimination(Matrix A, int n){
 3     int i, j, k, r;
 4     //消元过程
 5     for(i = 0; i < n; ++i){
 6         //选一行r并与i行交换
 7         r = i;
 8         for(j = i+1; j < n; ++j){
 9             if(fabs(A[j][i]) > fabs(A[r][i])) r = j;
10         }
11         if(r != i){
12             for(j = 0; j <= n; ++j) swap(A[r][j], A[i][j]);
13         }
14         //与i+1~n行进行消元
15         for(k = i+1; k < n; ++k){
16             for(j = n; j >= i; --j){ //必须逆序枚举
17                 A[k][j] -= A[k][i]/A[i][i]*A[i][j];
18             }
19         }
20     }
21     //回代过程
22     for(i = n-1; i >= 0; --i){
23         for(j = i+1; j < n; ++j){
24             A[i][n] -= A[j][n] * A[i][j];
25         }
26         A[i][n] /= A[i][i];
27     }
28 }

高斯——约当消元法

　　运算量比高斯消元略大（将系数矩阵化为对角矩阵），但是代码更简单（少了回调过程）

 1 typedef double Matrix[maxn][maxn];
 2 const double eps = 1e-8;
 3 void gauss_jordan(Matrix A, int n){
 4     int i, j, k, r;
 5     for(i = 0; i < n; i++){
 6         r = i;
 7         for(j = i+1; j < n; j++){
 8             if(fabs(A[j][i]) > fabs(A[r][i])) r = j;
 9         }
10         if(fabs(A[r][i]) < eps) continue;
11         if(r != i){
12             for(j = 0; j <= n; j++){
13                 swap(A[r][j], A[i][j]);
14             }
15         }
16         for(k = 0; k < n; k++) if(k != i){
17             for(j = n; j >= i; j--){
18                 A[k][j] -= A[k][i]/A[i][i] * A[i][j];
19             }
20         }
21     }
22 }

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&category=20&problem=1769&mosmsg=Submission+received+with+ID+18927455

随机程序

求每个节点的期望执行次数

设　i　的出度为　di　期望执行次数为　xi　对于每个有三个前驱点　a　b　c　的节点　i　可以列出方程　xi　=　xa/da　+　xb/db　+　xc/dc　.

所以矩阵就可以构建出来

 1 #include <algorithm>
 2 #include <cmath>
 3 #include <cstdio>
 4 #include <cstring>
 5 #include <vector>
 6 using namespace std;
 7 
 8 const double eps = 1e-8;
 9 const int maxn = 110;
10 typedef double Ma[maxn][maxn];
11 
12 Ma A;
13 int n, d[maxn];
14 vector<int> p[maxn];
15 int inf[maxn];
16 
17 
18 void gass(Ma A, int n){
19     int i, j, k, r;
20     for(i = 0; i < n; ++i){
21         r = i;
22         for(j = i+1; j < n; ++j){
23             if(fabs(A[j][i]) > fabs(A[r][i])) r = j;
24         }        
25         if(fabs(A[r][i]) < eps) continue;
26         if(r != i){
27             for(j = 0; j <= n; ++j) swap(A[r][j], A[i][j]);
28         }
29         for(k = 0; k < n; ++k) if(k != i){
30             for(j = n; j >= i; j--) A[k][j] -= A[k][i] / A[i][i] * A[i][j];
31         }
32     }
33 }
34 
35 
36 int main(){
37     int kase = 0;
38     while(scanf("%d", &n) == 1 && n){
39         memset(d, 0, sizeof(d));
40         for(int i = 0; i < n; ++i) p[i].clear();
41         int a, b;
42         while(scanf("%d%d", &a, &b) == 2 && a){
43             a--, b--;
44             d[a]++;
45             p[b].push_back(a);
46         }
47         memset(A, 0, sizeof(A));
48         for(int i = 0; i < n; ++i){
49             A[i][i] = 1;
50             for(int j = 0; j < p[i].size(); j++){
51                 int xx = p[i][j];
52                 A[i][xx] -= 1.0 / d[xx];
53             }
54             if(i == 0) A[i][n] = 1;
55         }
56         gass(A, n);
57         memset(inf, 0, sizeof(inf));
58         for(int i = n-1; i >= 0; --i){
59             if(fabs(A[i][i]) < eps && fabs(A[i][n]) > eps) inf[i] = 1;
60             for(int j = i+1; j < n; j++){
61                 if(fabs(A[i][j]) > eps && inf[j]) inf[i] = 1;
62             }
63         }
64         
65         int q, u;
66         scanf("%d", &q);
67         printf("Case #%d:\n", ++kase);
68         while(q--){
69             scanf("%d", &u);
70             u--;
71             if(inf[u]) printf("infinity\n");
72             else printf("%.3lf\n", fabs(A[u][u]) < eps ? 0.0 : A[u][n]/A[u][u]);
73         }
74         
75     }
76     return 0;
77 }

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

kuangbin大神模板

  1 #include<stdio.h>
  2 #include<algorithm>
  3 #include<iostream>
  4 #include<string.h>
  5 #include<math.h>
  6 using namespace std;
  7 
  8 const int MAXN=50;
  9 
 10 
 11 
 12 int a[MAXN][MAXN];//增广矩阵
 13 int x[MAXN];//解集
 14 bool free_x[MAXN];//标记是否是不确定的变元
 15 
 16 
 17 
 18 /*
 19 void Debug(void)
 20 {
 21     int i, j;
 22     for (i = 0; i < equ; i++)
 23     {
 24         for (j = 0; j < var + 1; j++)
 25         {
 26             cout << a[i][j] << " ";
 27         }
 28         cout << endl;
 29     }
 30     cout << endl;
 31 }
 32 */
 33 
 34 
 35 inline int gcd(int a,int b)
 36 {
 37     int t;
 38     while(b!=0)
 39     {
 40         t=b;
 41         b=a%b;
 42         a=t;
 43     }
 44     return a;
 45 }
 46 inline int lcm(int a,int b)
 47 {
 48     return a/gcd(a,b)*b;//先除后乘防溢出
 49 }
 50 
 51 // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，
 52 //-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
 53 //有equ个方程，var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
 54 int Gauss(int equ,int var)
 55 {
 56     int i,j,k;
 57     int max_r;// 当前这列绝对值最大的行.
 58     int col;//当前处理的列
 59     int ta,tb;
 60     int LCM;
 61     int temp;
 62     int free_x_num;
 63     int free_index;
 64 
 65     for(int i=0;i<=var;i++)
 66     {
 67         x[i]=0;
 68         free_x[i]=true;
 69     }
 70 
 71     //转换为阶梯阵.
 72     col=0; // 当前处理的列
 73     for(k = 0;k < equ && col < var;k++,col++)
 74     {// 枚举当前处理的行.
 75 // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
 76         max_r=k;
 77         for(i=k+1;i<equ;i++)
 78         {
 79             if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
 80         }
 81         if(max_r!=k)
 82         {// 与第k行交换.
 83             for(j=k;j<var+1;j++) swap(a[k][j],a[max_r][j]);
 84         }
 85         if(a[k][col]==0)
 86         {// 说明该col列第k行以下全是0了，则处理当前行的下一列.
 87             k--;
 88             continue;
 89         }
 90         for(i=k+1;i<equ;i++)
 91         {// 枚举要删去的行.
 92             if(a[i][col]!=0)
 93             {
 94                 LCM = lcm(abs(a[i][col]),abs(a[k][col]));
 95                 ta = LCM/abs(a[i][col]);
 96                 tb = LCM/abs(a[k][col]);
 97                 if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
 98                 for(j=col;j<var+1;j++)
 99                 {
100                     a[i][j] = a[i][j]*ta-a[k][j]*tb;
101                 }
102             }
103         }
104     }
105 
106   //  Debug();
107 
108     // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
109     for (i = k; i < equ; i++)
110     { // 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.
111         if (a[i][col] != 0) return -1;
112     }
113     // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
114     // 且出现的行数即为自由变元的个数.
115     if (k < var)
116     {
117         // 首先，自由变元有var - k个，即不确定的变元至少有var - k个.
118         for (i = k - 1; i >= 0; i--)
119         {
120             // 第i行一定不会是(0, 0, ..., 0)的情况，因为这样的行是在第k行到第equ行.
121             // 同样，第i行一定不会是(0, 0, ..., a), a != 0的情况，这样的无解的.
122             free_x_num = 0; // 用于判断该行中的不确定的变元的个数，如果超过1个，则无法求解，它们仍然为不确定的变元.
123             for (j = 0; j < var; j++)
124             {
125                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
126             }
127             if (free_x_num > 1) continue; // 无法求解出确定的变元.
128             // 说明就只有一个不确定的变元free_index，那么可以求解出该变元，且该变元是确定的.
129             temp = a[i][var];
130             for (j = 0; j < var; j++)
131             {
132                 if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
133             }
134             x[free_index] = temp / a[i][free_index]; // 求出该变元.
135             free_x[free_index] = 0; // 该变元是确定的.
136         }
137         return var - k; // 自由变元有var - k个.
138     }
139     // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
140     // 计算出Xn-1, Xn-2 ... X0.
141     for (i = var - 1; i >= 0; i--)
142     {
143         temp = a[i][var];
144         for (j = i + 1; j < var; j++)
145         {
146             if (a[i][j] != 0) temp -= a[i][j] * x[j];
147         }
148         if (temp % a[i][i] != 0) return -2; // 说明有浮点数解，但无整数解.
149         x[i] = temp / a[i][i];
150     }
151     return 0;
152 }
153 int main(void)
154 {
155     freopen("in.txt", "r", stdin);
156     freopen("out.txt","w",stdout);
157     int i, j;
158     int equ,var;
159     while (scanf("%d %d", &equ, &var) != EOF)
160     {
161         memset(a, 0, sizeof(a));
162         for (i = 0; i < equ; i++)
163         {
164             for (j = 0; j < var + 1; j++)
165             {
166                 scanf("%d", &a[i][j]);
167             }
168         }
169 //        Debug();
170         int free_num = Gauss(equ,var);
171         if (free_num == -1) printf("无解!\n");
172    else if (free_num == -2) printf("有浮点数解，无整数解!\n");
173         else if (free_num > 0)
174         {
175             printf("无穷多解! 自由变元个数为%d\n", free_num);
176             for (i = 0; i < var; i++)
177             {
178                 if (free_x[i]) printf("x%d 是不确定的\n", i + 1);
179                 else printf("x%d: %d\n", i + 1, x[i]);
180             }
181         }
182         else
183         {
184             for (i = 0; i < var; i++)
185             {
186                 printf("x%d: %d\n", i + 1, x[i]);
187             }
188         }
189         printf("\n");
190     }
191     return 0;
192 }

----------------------------------------------------------------------------------------------------------------

消元

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cmath>
 4 
 5 using namespace std;
 6 
 7 typedef __int64 lld;
 8 
 9 lld a[205][205];
10 
11 int sign;
12 lld N,MOD;
13 void solved()
14 {
15     lld ans=1;
16     for(int i=0;i<N;i++)//当前行
17     {
18         for(int j=i+1;j<N;j++)//当前之后的每一行，因为每一行的当前第一个数要转化成0（想想线性代数中行列式的计算）
19         {
20             int x=i,y=j;
21             while(a[y][i])//利用gcd的方法，不停地进行辗转相除
22             {
23                 lld t=a[x][i]/a[y][i];
24 
25                 for(int k=i;k<N;k++)
26                     a[x][k]=(a[x][k]-a[y][k]*t)%MOD;
27 
28                 swap(x,y);
29             }
30             if(x!=i)//奇数次交换，则D=-D'整行交换
31             {
32                 for(int k=0;k<N;k++)
33                     swap(a[i][k],a[x][k]);
34                 sign^=1;
35             }
36         }
37         if(a[i][i]==0)//斜对角中有一个0，则结果为0
38         {
39             cout<<0<<endl;
40             return ;
41         }
42 
43         else
44             ans=ans*a[i][i]%MOD;
45 
46     }
47 
48     if(sign!=0)
49         ans*=-1;
50     if(ans<0)
51         ans+=MOD;
52     printf("%I64d\n",ans);
53 }
54 int main()
55 {
56     int t;
57     scanf("%d",&t);
58 
59     while(t--)
60     {
61         sign=0;
62         scanf("%I64d%I64d",&N,&MOD);
63         for(int i=0;i<N;i++)
64             for(int j=0;j<N;j++)
65                 scanf("%I64d",&a[i][j]);
66         solved();
67     }
68     return 0;
69 }

----------------------------------------------------------------------------------------------------------------

湘潭邀请赛

http://202.197.224.59/OnlineJudge2/index.php/Problem/read/id/1260

A-1 = A*/|A|

 1 #include <string.h>
 2 #include <iostream>
 3 #include <stdio.h>
 4 using namespace std;
 5 #define ll long long
 6 const int Mod = 1e9+7;
 7 int n;
 8 int a[210][210], b[210][210];
 9 int mul(int x){
10     int xx = Mod - 2, sum = 1;
11     while(xx){
12         if(xx&1) sum = 1LL * sum * x % Mod;
13         x = 1LL * x * x % Mod;
14         xx >>= 1;
15     }
16     return sum;
17 }
18 void d(){
19     int cnt = 1;
20     for(int i = 0; i < n; ++i){
21         for(int j = 0; j < n; ++j){
22             b[i][j] = (i == j);
23         }
24     }
25     for(int i = 0; i < n; ++i){
26         int t = i;
27         for(int j = i; j < n; ++j){
28             if(!a[j][i]){
29                 t = j;
30             }
31         }
32         if(t != i) cnt *= -1;
33         for(int j = 0; j < n; ++j){
34             swap(a[i][j], a[t][j]);
35             swap(b[i][j], b[t][j]);
36         }
37         cnt = 1LL * cnt * a[i][i] % Mod;
38         int xx = mul(a[i][i]);  //求逆
39         for(int j = 0; j < n; ++j){
40             a[i][j] = 1LL * a[i][j]*xx%Mod;
41             b[i][j] = 1LL * b[i][j]*xx%Mod;
42         }
43         for(int k = 0; k < n; ++k){
44             if(k == i) continue;
45             ll tm = a[k][i];
46             for(int j = 0; j < n; ++j){
47                 a[k][j] = (a[k][j] - 1LL * tm*a[i][j]%Mod + Mod)%Mod;
48                 b[k][j] = (b[k][j] - 1LL * tm*b[i][j]%Mod + Mod)%Mod;
49             }
50         }
51     }
52     cnt = (Mod+cnt)%Mod;
53     for(int i = 0; i < n; ++i){
54         for(int j = 0; j < n; ++j){
55             b[i][j] = 1LL * b[i][j]*cnt%Mod;
56         }
57     }
58 }
59 int main(){
60     while(~scanf("%d", &n)){
61         for(int i = 0; i < n; ++i){
62             a[0][i] = 1;
63         }
64         for(int i = 1; i < n; ++i){
65             for(int j = 0; j < n; ++j){
66                 scanf("%d", &a[i][j]);
67             }
68         }
69         d();
70         for(int i = 0; i < n; ++i){
71             printf("%d%c", (i&1 ? (Mod-b[i][0])%Mod : b[i][0]), (i == n-1 ? '\n' : ' '));
72         }
73     }
74     return 0;
75 }

----------------------------------------------------------------------------------------------------------------

https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=27&page=show_problem&problem=2537

乘积是平方数

 1 #include <iostream>
 2 #include <string.h>
 3 using namespace std;
 4 #define ll long long
 5 const int maxn = 555;
 6 int num = 0, n, ma;
 7 bool vis[maxn];
 8 int pr[maxn], A[185][105];
 9 void p(){
10     for(int i = 2; i <= 500; ++i){
11         if(!vis[i]) pr[num++] = i;
12         for(int j = 0; j < num && i*pr[j] <= 500; ++j){
13             if(i%pr[j]) vis[i*pr[j]] = true;
14             else{
15                 vis[i*pr[j]] = true;
16                 break;
17             }
18         }
19     }
20 }
21 void solve(){
22     int i = 0, j = 0, k, r, u;
23     while(i < ma && j < n){
24         r = i;
25         for(k = i; k < ma; ++k){
26             if(A[k][j]){
27                 r = k;
28                 break;
29             }
30         }
31         if(A[r][j]){
32             if(r != i){
33                 for(k = 0; k <= n; ++k) swap(A[r][k], A[i][k]);
34             }
35             for(k = i+1; k < ma; ++k) if(A[k][j]){
36                 for(u = i; u <= n; ++u){
37                     A[k][u] ^= A[i][u];
38                 }
39             }
40             ++i;
41         }
42         ++j;
43     }
44     //cout<<"****************"<<i<<endl;
45     printf("%lld\n", (1LL << (n-i))-1LL);
46 }
47 int main(){
48     int t;
49     p();
50     //cout<<num<<endl;
51     //cout<<pr[0]<<endl;
52     scanf("%d", &t);
53     while(t--){
54         memset(A, 0, sizeof(A));
55         scanf("%d", &n);
56         ma = -1;
57         for(int i = 0; i < n; ++i){
58             ll x;
59             scanf("%lld", &x);
60             for(int j = 0; j < num; ++j){
61                 while(x%pr[j] == 0){
62                     A[j][i] ^= 1, ma = max(ma, j+1), x/=pr[j];
63                 }
64             }
65         }
66         /*
67         for(int i = 0; i < ma; ++i){
68             for(int j = 0; j < n; ++j){
69                 printf("%d%c", A[i][j], (j == n-1 ? '\n' : ' '));
70             }
71         }
72         */
73         solve();
74     }
75     return 0;
76 }

----------------------------------------------------------------------------------------------------------------

只有不断学习才能进步！

posted @ 2018-04-14 13:53 wenbao 阅读(190) 评论(0) 编辑收藏举报

刷新页面返回顶部

wenbao

只有不断学习才能进步！新博客地址(http://oabnew.com/myblog)

wenbao与高斯消元

公告