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 */ getField()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 */ getField()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 * Computes a matrix with rows forming a lattice basis of the lattice 250 * kernel of the current matrix assuming that its entries are integral. 251 */ 252 FieldMatrix basisOfLatticeKernel()const; 253 /** 254 Computes and returns the inverse of the matrix. The matrix must be invertible. 255 */ 256 FieldMatrix inverse()const; 257 /** 258 Returns the transposed of the the matrix. 259 */ 260 FieldMatrix transposed()const; 261 /** 262 Returns the matrix with rows in reversed order. 263 */ 264 FieldMatrix flipped()const; 265 /** 266 Takes two matrices with the same number of columns and construct 267 a new matrix which has the rows of the matrix top on the top and 268 the rows of the matrix bottom at the bottom. The return value is 269 the constructed matrix. 270 */ 271 friend FieldMatrix combineOnTop(FieldMatrix const &top, FieldMatrix const &bottom); 272 /** 273 Takes two matrices with the same number of rows and construct 274 a new matrix which has the columns of the matrix left on the left and 275 the columns of the matrix right on the right. The return value is 276 the constructed matrix. 277 */ 278 friend FieldMatrix combineLeftRight(FieldMatrix const &left, FieldMatrix const &right); 279 /** 280 Computes a reduced row echelon form from a row echelon form. In 281 Groebner basis terms this is the same as tranforming a minimal 282 Groebner basis to a reduced one except that we do not force 283 pivots to be 1 unless the scalePivotsToOne parameter is set. 284 */ 285 int REformToRREform(bool scalePivotsToOne=false); 286 /** 287 This function may be called if the FieldMatrix is in Row Echelon 288 Form. The input is a FieldVector which is rewritten modulo the 289 rows of the matrix. The result is unique and is the same as the 290 normal form of the input vector modulo the Groebner basis of the 291 linear ideal generated by the rows of the matrix. 292 */ 293 FieldVector canonicalize(FieldVector v)const; 294 /** 295 This routine removes the zero rows of the matrix. 296 */ 297 void removeZeroRows(); 298 /** 299 This routine produces a matrix in row Echelon form solving a 300 system of the kind Ax=b, where A equals *this. The way this is 301 done is by considering [A^t -I] and computing its RE form. To 302 solve the system afterwards simply call canonicalize of the new 303 matrix on (b_1,\dots,b_m,0,\dots,0) to get 304 (0,\dots,0,x_1,\dots,x_n). If no solution exists then the first 305 m entries of the returned vector are not all 0. 306 */ 307 FieldMatrix solver()const; 308 /** 309 Returns the matrix sum a*b. 310 */ 311 friend FieldMatrix operator+(FieldMatrix const &a, const FieldMatrix &b); 312 /** 313 Computes the vector m scaled by a factor s. 314 */ 315 friend FieldMatrix operator*(FieldElement s, const FieldMatrix& m); 316 /** 317 Returns the matrix product a*b. 318 */ 319 friend FieldMatrix operator*(FieldMatrix const &a, const FieldMatrix &b); 320 /** 321 Returns the vector a*b, where a is considered as a row vector. 322 */ 323 friend FieldVector operator*(FieldVector const &a, const FieldMatrix &b); 324 /** 325 Compares two matrices. 326 */ 327 bool operator==(FieldMatrix const &b)const; 328 /** 329 Returns the rank of the matrix. This function is const and slow - 330 will copy the matrix and perform Gauss reduction. 331 */ 332 int rank()const; 333 /** 334 Returns the indices of the columns containing a pivot. 335 The returned list is sorted. 336 The matrix must be in row echelon form. 337 */ 338 list<int> pivotColumns()const; 339 /** 340 Returns the indices of the columns not containing a pivot. 341 The returned list is sorted. 342 The matrix must be in row echelon form. 343 */ 344 list<int> nonPivotColumns()const; 345 /** 346 * Returns the specified submatrix with size (endRow-startRow)x(endColumn-startColumn). 347 */ 348 FieldMatrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const; 349 /** 350 * Thinking of an upper triangular matrix as a Groebner basis this routine returns the normal form 351 * of the linear polynomial represented by v. If coeff is non-zero, then that vector is assigned 352 * the coefficients used in the reduction. That is, if the normal form is zero, then multiplying coeff 353 * and *this gives v. The matrix must be in row-echelon form. coeff does not have to have the right length before the call. 354 */ 355 FieldVector normalForm(FieldVector const &v, FieldVector *coeff)const; 356 appendRow(FieldVector const & r)357 void appendRow(FieldVector const &r){assert(width==r.size());rows.push_back(r);height++;} 358 }; 359 360 361 FieldMatrix integerMatrixToFieldMatrix(IntegerMatrix const &m, Field f); 362 IntegerMatrix fieldMatrixToIntegerMatrix(FieldMatrix const &m);//MUST CONTAIN INTEGER ENTRIES ONLY 363 IntegerMatrix fieldMatrixToIntegerMatrixPrimitive(FieldMatrix const &m);//Scales to primitive vector 364 FieldVector integerVectorToFieldVector(IntegerVector const &v, Field f); 365 366 /** 367 Computes a non-zero vector in the kernel of the given matrix. The 368 dimension of the kernel is assumed to be 1. The returned vector is primitive. 369 */ 370 IntegerVector vectorInKernel(IntegerMatrix const &m); 371 372 /** 373 * Using exact arithmetics this routine chooses a subset of the vectors in m, which generate the span of m. 374 */ 375 IntegerVectorList subsetBasis(IntegerVectorList const &m); 376 377 /** 378 * Computes a basis of the quotient lattice of two saturated lattices. 379 * The lattice kernel of "equations" is a lattice A which contains the sublattice B which is the 380 * lattice kernel of "euqations union additionalEquations". 381 * The routine returns a minimal set of vectors which will extend a lattice basis of B to a lattice basis of A. 382 */ 383 IntegerVectorList basisOfQuotientLattice(IntegerVectorList const &equations, IntegerVectorList const &additionalEquations, int n); 384 385 386 387 #endif 388