1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
/*++
Copyright (c) 2012 Microsoft Corporation
Module Name:
linear_eq_solver.h
Abstract:
Simple equational solver template for any number field.
No special optimization, just the basics for solving small systems.
It is a solver target to dense system of equations.
Main client: Sparse Modular GCD algorithm.
Author:
Leonardo (leonardo) 2012-01-22
Notes:
--*/
#pragma once
template<typename numeral_manager>
class linear_eq_solver {
typedef typename numeral_manager::numeral numeral;
numeral_manager & m;
unsigned n; // number of variables
vector<svector<numeral> > A;
svector<numeral> b;
public:
linear_eq_solver(numeral_manager & _m):m(_m), n(0) { SASSERT(m.field()); }
~linear_eq_solver() { flush(); }
void flush() {
SASSERT(b.size() == A.size());
auto sz = A.size();
for (unsigned i = 0; i < sz; ++i) {
svector<numeral> & as = A[i];
m.del(b[i]);
SASSERT(as.size() == n);
for (unsigned j = 0; j < n; ++j)
m.del(as[j]);
}
A.reset();
b.reset();
n = 0;
}
void resize(unsigned _n) {
if (n != _n) {
flush();
n = _n;
for (unsigned i = 0; i < n; ++i) {
A.push_back(svector<numeral>());
svector<numeral> & as = A.back();
for (unsigned j = 0; j < n; ++j) {
as.push_back(numeral());
}
b.push_back(numeral());
}
}
}
void reset() {
for (unsigned i = 0; i < n; ++i) {
svector<numeral> & A_i = A[i];
for (unsigned j = 0; j < n; ++j) {
m.set(A_i[j], 0);
}
m.set(b[i], 0);
}
}
// Set row i with _as[0]*x_0 + ... + _as[n-1]*x_{n-1} = b
void add(unsigned i, numeral const * _as, numeral const & _b) {
SASSERT(i < n);
m.set(b[i], _b);
svector<numeral> & A_i = A[i];
for (unsigned j = 0; j < n; ++j) {
m.set(A_i[j], _as[j]);
}
}
// Return true if the system of equations has a solution.
// Return false if the matrix is singular
bool solve(numeral * xs) {
for (unsigned k = 0; k < n; ++k) {
TRACE(linear_eq_solver, tout << "iteration " << k << "\n"; display(tout););
// find pivot
unsigned i = k;
for (; i < n; ++i) {
if (!m.is_zero(A[i][k]))
break;
}
if (i == n)
return false; // matrix is singular
A[k].swap(A[i]); // swap rows
svector<numeral> & A_k = A[k];
numeral & A_k_k = A_k[k];
SASSERT(!m.is_zero(A_k_k));
// normalize row
for (unsigned i = k+1; i < n; ++i)
m.div(A_k[i], A_k_k, A_k[i]);
m.div(b[k], A_k_k, b[k]);
m.set(A_k_k, 1);
// check if first k-1 positions are zero
DEBUG_CODE({ for (unsigned i = 0; i < k; ++i) { SASSERT(m.is_zero(A_k[i])); } });
// for all rows below pivot
for (unsigned i = k+1; i < n; ++i) {
svector<numeral> & A_i = A[i];
numeral & A_i_k = A_i[k];
for (unsigned j = k+1; j < n; ++j) {
m.submul(A_i[j], A_i_k, A_k[j], A_i[j]);
}
m.submul(b[i], A_i_k, b[k], b[i]);
m.set(A_i_k, 0);
}
}
unsigned k = n;
while (k > 0) {
--k;
TRACE(linear_eq_solver, tout << "iteration " << k << "\n"; display(tout););
SASSERT(m.is_one(A[k][k]));
// save result
m.set(xs[k], b[k]);
// back substitute
unsigned i = k;
while (i > 0) {
--i;
m.submul(b[i], A[i][k], b[k], b[i]);
m.set(A[i][k], 0);
}
}
return true;
}
void display(std::ostream & out) const {
for (unsigned i = 0; i < A.size(); ++i) {
SASSERT(A[i].size() == n);
for (unsigned j = 0; j < n; ++j) {
m.display(out, A[i][j]);
out << " ";
}
m.display(out, b[i]); out << "\n";
}
}
};