消元
高斯消元
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 }
----------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------
只有不断学习才能进步！
原文地址：https://www.cnblogs.com/wenbao/p/6498612.html
时间： 2024-10-13 02:54:50