~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to linalg.h

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#ifndef LINALG_H_INCLUDED
 
2
#define LINALG_H_INCLUDED
 
3
 
 
4
#include <vector>
 
5
#include "field.h"
 
6
#include "field_rationals.h"
 
7
#include "printer.h"
 
8
#include "matrix.h"
 
9
 
 
10
/*
 
11
  Unfortunately it seems that these classes cannot share code with the
 
12
  Vektor and Matrix templates since the classes in this file carry
 
13
  typeinfo around. In particular it is not possible to use the Matrix
 
14
  template on the FieldElement type since a FieldElement has no
 
15
  default initialization from a double -- its initialization really
 
16
  depends on which Field we are using.
 
17
 */
 
18
 
 
19
using namespace std;
 
20
 
 
21
/**
 
22
   @brief A vector containing FieldElement s.
 
23
 
 
24
   When a FieldVector is constructed the Field of which the entries of
 
25
   the vector will be members must be specified either explicitly or
 
26
   implicitly.
 
27
*/
 
28
class FieldVector
 
29
{
 
30
  /**
 
31
     The field of which the entries of the vector are members.
 
32
   */
 
33
  Field theField;
 
34
  /**
 
35
     The length of the vector.
 
36
  */
 
37
  int n;
 
38
  /**
 
39
     The entries of the vector.
 
40
   */
 
41
  vector<FieldElement> v;
 
42
 public:
 
43
   /**
 
44
     Constructs a FieldVector with entries initialized to 0.
 
45
 
 
46
     @param _theField the Field of which the entries of the vector will be members.
 
47
     @param _n the number of entries in the vector.
 
48
   */
 
49
  FieldVector(Field const &_theField, int _n);
 
50
  /**
 
51
   * Produces a vector with all entries equal to 1.
 
52
   */
 
53
  static FieldVector allOnes(Field const &_theField, int _n);
 
54
  /**
 
55
      Returns the Field.
 
56
  */
 
57
  inline Field const&getField()const{return theField;}
 
58
  /**
 
59
     Provides access to i th entry of the vector.
 
60
   */
 
61
  FieldElement& operator[](int i);
 
62
  /**
 
63
     Provides access to i th entry of the vector.
 
64
   */
 
65
  FieldElement const& operator[](int i)const;
 
66
  /**
 
67
     Returns true if all entries of the vector equal the 0-element of the Field, false otherwise.
 
68
   */
 
69
  bool isZero()const;
 
70
  /**
 
71
     Computes the first integer vector on the half ray spanned by the
 
72
     vector. Assumes that the Field is Q, asserts if this is not the
 
73
     case. Asserts if the computed vector has entries that do not fit
 
74
     in the word size.
 
75
   */
 
76
  IntegerVector primitive()const;
 
77
  /**
 
78
     Returns the number of entries in the vector.
 
79
   */
 
80
  int size()const;
 
81
  /**
 
82
     Returns the subvector (v_begin,\dots,v_end-1).
 
83
   */
 
84
  FieldVector subvector(int begin, int end)const;
 
85
  /**
 
86
     Returns the subvector .
 
87
   */
 
88
  FieldVector subvector(list<int> const &l)const;
 
89
  /**
 
90
     Returns the concatenation of the two vectors.
 
91
   */
 
92
  friend FieldVector concatenation(FieldVector const &a, FieldVector const &b);
 
93
  /**
 
94
     Computes the dot product of the vectors a and b.
 
95
   */
 
96
  friend FieldElement dot(FieldVector const &a, FieldVector const &b);
 
97
  /**
 
98
     Computes the vector q scaled by a factor s.
 
99
   */
 
100
  friend FieldVector operator*(FieldElement s, const FieldVector& q);
 
101
  /**
 
102
     Computes the coordinate-wise subtraction a-b.
 
103
   */
 
104
  friend FieldVector operator-(FieldVector const &a, FieldVector const &b);
 
105
  /**
 
106
     Computes the coordinate-wise division of a by b.
 
107
   */
 
108
  friend FieldVector operator/(FieldVector const &a, FieldVector const &b);
 
109
  /**
 
110
     Adds q to the current vector coordinatewise. Asserts if the
 
111
     current vector and q do not have the same number of elements.
 
112
   */
 
113
  FieldVector& operator+=(const FieldVector& q);
 
114
  /**
 
115
     Adds s*q to the current vector coordinatewise. Asserts if the
 
116
     current vector and q do not have the same number of elements.
 
117
   */
 
118
  FieldVector& madd(const FieldElement &s, const FieldVector& q);
 
119
  /**
 
120
     Prints the vector using the Printer P.
 
121
   */
 
122
  void print(Printer &P)const;
 
123
 
 
124
  friend IntegerVector fieldVectorToIntegerVector(FieldVector const &v);//MUST CONTAIN INTEGER ENTRIES ONLY // should this be a member?
 
125
  bool operator==(FieldVector const &b)const;
 
126
};
 
127
 
 
128
 
 
129
/**
 
130
   @brief A matrix containing FieldElement s.
 
131
 
 
132
   When a FieldMatrix is constructed the Field of which the entries of
 
133
   the matrix will be members must be specified either explicitly or
 
134
   implicitly. The FieldMatrix class contains routines for performing
 
135
   Gauss reductions.
 
136
*/
 
137
class FieldMatrix
 
138
{
 
139
  Field theField;
 
140
  int width,height;
 
141
  vector<FieldVector> rows;
 
142
 public:
 
143
   /**
 
144
     Constructs a FieldMatrix with entries initialized to 0.
 
145
 
 
146
     @param _theField the Field of which the entries of the matrix will be members.
 
147
     @param _height the number of rows of the matrix.
 
148
     @param _width the number of columns of the matrix.
 
149
   */
 
150
  FieldMatrix(Field const &_theField, int _height, int _width);
 
151
  /**
 
152
     Returns the nxn identity matrix.
 
153
  */
 
154
  static FieldMatrix identity(Field const &_theField, int _n);
 
155
  /**
 
156
     Returns the Field.
 
157
  */
 
158
  inline Field const&getField()const{return theField;}
 
159
  /**
 
160
     Provides access to i th row of the matrix which is a FieldVector.
 
161
  */
 
162
  FieldVector& operator[](int i);
 
163
  /**
 
164
     Provides access to i th row of the matrix which is a FieldVector.
 
165
  */
 
166
  FieldVector const& operator[](int i)const;
 
167
  /**
 
168
     Returns the number of rows of the matrix.
 
169
  */
 
170
  int getHeight()const;
 
171
  /**
 
172
     Returns the number of columns of the matrix.
 
173
  */
 
174
  int getWidth()const;
 
175
  /**
 
176
     Swaps the i th and the j th row.
 
177
   */
 
178
  void swapRows(int i, int j);
 
179
  /**
 
180
     Returns the submatrix consisting of rows index by rowIndices.
 
181
   */
 
182
  FieldMatrix submatrixRows(list<int> const &rowIndices)const;
 
183
  /**
 
184
     Adds a times the i th row to the j th row.
 
185
  */
 
186
  void madd(int i, FieldElement a, int j);
 
187
  /**
 
188
     Returns the index of a row whose index is at least currentRow and
 
189
     which has a non-zero element on the column th entry. If no such
 
190
     row exists then -1 is returned. This routine is used in the Gauss
 
191
     reduction. To make the reduction more efficient the routine
 
192
     chooses its row with as few non-zero entries as possible.
 
193
   */
 
194
  int findRowIndex(int column, int currentRow)const;
 
195
  /**
 
196
     Prints the matrix using the Printer P.
 
197
   */
 
198
  void printMatrix(Printer &P)const;
 
199
  /**
 
200
     This method is used for iterating through the pivots in a matrix
 
201
     in row echelon form. To find the first pivot put i=-1 and
 
202
     j=-1 and call this routine. The (i,j) th entry of the matrix is a
 
203
     pivot. Call the routine again to get the next pivot. When no more
 
204
     pivots are found the routine returns false.
 
205
  */
 
206
  bool nextPivot(int &i, int &j)const;
 
207
  /**
 
208
     Performs a Gauss reduction and returns the number of row swaps
 
209
     done.  The result is a matrix in row echelon form. The pivots may
 
210
     not be all 1.  In terms of Groebner bases, what is computed is a
 
211
     minimal (not necessarily reduced) Groebner basis of the linear
 
212
     ideal generated by the rows.  The number of swaps is need if one
 
213
     wants to compute the determinant afterwards. In this case it is
 
214
     also a good idea to set the flag returnIfZeroDeterminant which
 
215
     make the routine terminate before completion if it discovers that
 
216
     the determinant is zero.
 
217
  */
 
218
  int reduce(bool returnIfZeroDeterminant=false, bool hermite=false); //bool reducedRowEcholonForm,
 
219
  /**
 
220
     Calls reduce() and returns the determinant of the matrix. If
 
221
     the matrix is not square the routine asserts.
 
222
   */
 
223
  FieldElement reduceAndComputeDeterminant();
 
224
  /**
 
225
     Calls reduce() and returns the rank of the matrix.
 
226
   */
 
227
  int reduceAndComputeRank();
 
228
  /**
 
229
     Scales every row of the matrix so that the first non-zero entry of
 
230
     each row is 1.
 
231
   */
 
232
  void scaleFirstNonZeroEntryToOne();
 
233
  /**
 
234
     Assumes that the matrix has a kernel of dimension 1.
 
235
     Calls reduce() and returns a non-zero vector in the kernel.
 
236
     If the matrix is an (n-1)x(n) matrix then the returned vector has
 
237
     the property that if appended as a row to the original matrix
 
238
     then the determinant of that matrix would be positive. Of course
 
239
     this property, as described here, only makes sense for ordered fields.
 
240
   */
 
241
  FieldVector reduceAndComputeVectorInKernel();
 
242
  /**
 
243
     Calls reduce() and constructs matrix whose rows forms a basis of
 
244
     the kernel of the linear map defined by the original matrix. The
 
245
     return value is the new matrix.
 
246
   */
 
247
  FieldMatrix reduceAndComputeKernel();
 
248
  /**
 
249
     Returns the transposed of the the matrix.
 
250
   */
 
251
  FieldMatrix transposed()const;
 
252
  /**
 
253
     Returns the matrix with rows in reversed order.
 
254
   */
 
255
  FieldMatrix flipped()const;
 
256
  /**
 
257
     Takes two matrices with the same number of columns and construct
 
258
     a new matrix which has the rows of the matrix top on the top and
 
259
     the rows of the matrix bottom at the bottom. The return value is
 
260
     the constructed matrix.
 
261
   */
 
262
  friend FieldMatrix combineOnTop(FieldMatrix const &top, FieldMatrix const &bottom);
 
263
  /**
 
264
     Takes two matrices with the same number of rows and construct
 
265
     a new matrix which has the columns of the matrix left on the left and
 
266
     the columns of the matrix right on the right. The return value is
 
267
     the constructed matrix.
 
268
   */
 
269
  friend FieldMatrix combineLeftRight(FieldMatrix const &left, FieldMatrix const &right);
 
270
  /**
 
271
     Computes a reduced row echelon form from a row echelon form. In
 
272
     Groebner basis terms this is the same as tranforming a minimal
 
273
     Groebner basis to a reduced one except that we do not force
 
274
     pivots to be 1 unless the scalePivotsToOne parameter is set.
 
275
   */
 
276
  int REformToRREform(bool scalePivotsToOne=false);
 
277
  /**
 
278
     This function may be called if the FieldMatrix is in Row Echelon
 
279
     Form. The input is a FieldVector which is rewritten modulo the
 
280
     rows of the matrix. The result is unique and is the same as the
 
281
     normal form of the input vector modulo the Groebner basis of the
 
282
     linear ideal generated by the rows of the matrix.
 
283
  */
 
284
  FieldVector canonicalize(FieldVector v)const;
 
285
  /**
 
286
     This routine removes the zero rows of the matrix.
 
287
   */
 
288
  void removeZeroRows();
 
289
  /**
 
290
     This routine produces a matrix in row Echelon form solving a
 
291
     system of the kind Ax=b, where A equals *this. The way this is
 
292
     done is by considering [A^t -I] and computing its RE form.  To
 
293
     solve the system afterwards simply call canonicalize of the new
 
294
     matrix on (b_1,\dots,b_m,0,\dots,0) to get
 
295
     (0,\dots,0,x_1,\dots,x_n).  If no solution exists then the first
 
296
     m entries of the returned vector are not all 0.
 
297
   */
 
298
  FieldMatrix solver()const;
 
299
  /**
 
300
     Computes the vector m scaled by a factor s.
 
301
   */
 
302
  friend FieldMatrix operator*(FieldElement s, const FieldMatrix& m);
 
303
  /**
 
304
     Returns the matrix product a*b.
 
305
   */
 
306
  friend FieldMatrix operator*(FieldMatrix const &a, const FieldMatrix &b);
 
307
  /**
 
308
     Returns the vector a*b, where a is considered as a row vector.
 
309
   */
 
310
  friend FieldVector operator*(FieldVector const &a, const FieldMatrix &b);
 
311
  /**
 
312
     Compares two matrices.
 
313
  */
 
314
  bool operator==(FieldMatrix const &b)const;
 
315
  /**
 
316
     Returns the rank of the matrix. This function is const and slow -
 
317
     will copy the matrix and perform Gauss reduction.
 
318
   */
 
319
  int rank()const;
 
320
  /**
 
321
     Returns the indices of the columns containing a pivot.
 
322
     The returned list is sorted.
 
323
     The matrix must be in row echelon form.
 
324
   */
 
325
  list<int> pivotColumns()const;
 
326
  /**
 
327
     Returns the indices of the columns not containing a pivot.
 
328
     The returned list is sorted.
 
329
     The matrix must be in row echelon form.
 
330
   */
 
331
  list<int> nonPivotColumns()const;
 
332
  /**
 
333
   * Returns the specified submatrix with size (endRow-startRow)x(endColumn-startColumn).
 
334
   */
 
335
  FieldMatrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const;
 
336
  /**
 
337
   * Thinking of an upper triangular matrix as a Groebner basis this routine returns the normal form
 
338
   * of the linear polynomial represented by v. If coeff is non-zero, then that vector is assigned
 
339
   * the coefficients used in the reduction. That is, if the normal form is zero, then multiplying coeff
 
340
   * and *this gives v. The matrix must be in row-echelon form. coeff does not have to have the right length before the call.
 
341
   */
 
342
  FieldVector normalForm(FieldVector const &v, FieldVector *coeff)const;
 
343
};
 
344
 
 
345
 
 
346
FieldMatrix integerMatrixToFieldMatrix(IntegerMatrix const &m, Field f);
 
347
IntegerMatrix fieldMatrixToIntegerMatrix(FieldMatrix const &m);//MUST CONTAIN INTEGER ENTRIES ONLY
 
348
IntegerMatrix fieldMatrixToIntegerMatrixPrimitive(FieldMatrix const &m);//Scales to primitive vector
 
349
FieldVector integerVectorToFieldVector(IntegerVector const &v, Field f);
 
350
 
 
351
/**
 
352
   Computes a non-zero vector in the kernel of the given matrix. The
 
353
   dimension of the kernel is assumed to be 1. The returned vector is primitive.
 
354
*/
 
355
IntegerVector vectorInKernel(IntegerMatrix const &m);
 
356
 
 
357
/**
 
358
 * Using exact arithmetics this routine chooses a subset of the vectors in m, which generate the span of m.
 
359
 */
 
360
IntegerVectorList subsetBasis(IntegerVectorList const &m);
 
361
 
 
362
#endif