source/kernel/Matop.c

   1 /*
   2     This file is part of PolyLib.
   3
   4     PolyLib is free software: you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation, either version 3 of the License, or
   7     (at your option) any later version.
   8
   9     PolyLib is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  16 */
  17
  18 #include <stdlib.h>
  19 #include <polylib/polylib.h>
  20
  21 /* computes c = lcm(a,b) using Gcd(a,b,&c) */
  22 void Lcm3(Value a, Value b, Value *c)
  23 {
  24   Value tmp;
  25
  26   if (value_zero_p(a)) {
  27     value_assign(*c, b);
  28     return;
  29   }
  30   if (value_zero_p(b)) {
  31     value_assign(*c, a);
  32     return;
  33   }
  34   value_init(tmp);
  35   value_multiply(tmp, a, b);
  36   value_absolute(tmp, tmp);
  37   Gcd(a,b,c);
  38   value_division(*c, tmp, *c);
  39   value_clear(tmp);
  40 }
  41
  42 /*
  43  * Return the lcm of 'i' and 'j'
  44  */
  45 Value *Lcm (Value i, Value j)
  46 {
  47   Value *tmp;
  48
  49   tmp = (Value *) malloc (sizeof(Value));
  50   value_init(*tmp);
  51   Lcm3(i, j, tmp);
  52   return tmp;
  53 } /* Lcm */
  54
  55 /*
  56  * Return an identity matrix of size 'size'.
  57  */
  58 Matrix *Identity(unsigned size)
  59 {
  60   unsigned i;
  61   Matrix *A;
  62
  63   A = Matrix_Alloc(size, size);
  64   for (i = 0; i < size; i++)
  65     value_set_si(A->p[i][i], 1);
  66   return A;
  67 } /* Identity */
  68
  69 /*
  70  * Exchange the rows 'Row1' and 'Row2' of the matrix 'M'
  71  */
  72 void ExchangeRows(Matrix *M, int Row1, int Row2) {
  73
  74   int i;
  75   Value temp;
  76
  77   value_init(temp);
  78   for (i = 0; i < (int)M->NbColumns;  i++) {
  79     value_assign(temp,M->p[Row1][i]);
  80     value_assign(M->p[Row1][i],M->p[Row2][i]);
  81     value_assign(M->p[Row2][i],temp);
  82   }
  83   value_clear(temp);
  84   return ;
  85 } /* ExchangeRows */
  86
  87 /*
  88  * Exchange the columns 'Column1' and 'Column2' of the matrix 'M'
  89  */
  90 void ExchangeColumns(Matrix *M, int Column1, int Column2) {
  91
  92   int i;
  93
  94   for (i = 0; i < M->NbRows; i++)
  95     value_swap(M->p[i][Column1],M->p[i][Column2]);
  96
  97   return;
  98 } /* ExchangeColumns */
  99
 100 /*
 101  * Return the Transpose of a matrix 'A'
 102  */
 103 Matrix *Transpose (Matrix *A) {
 104
 105   int i,j;
 106   Matrix *transpose;
 107
 108   transpose = Matrix_Alloc (A->NbColumns, A->NbRows);
 109   for (i = 0; i < (int)A->NbRows; i++)
 110     for (j = 0; j < (int)A->NbColumns; j++)
 111       value_assign(transpose->p[j][i],A->p[i][j]);
 112   return transpose;
 113 } /* Transpose */
 114
 115 /*
 116  * Return a copy of the contents of a matrix 'Src'.
 117  */
 118 Matrix *Matrix_Copy(Matrix const *Src )
 119 {
 120   Matrix *Dst;
 121   unsigned i, j;
 122
 123   Dst = Matrix_Alloc(Src->NbRows, Src->NbColumns);
 124
 125   for (i = 0; i < Src->NbRows; i++)
 126     for (j = 0; j < Src->NbColumns; j++)
 127       value_assign(Dst->p[i][j],Src->p[i][j]);
 128   return Dst;
 129 } /* Matrix_copy */
 130
 131 /*
 132  * Test if the input matrix is integral or not.
 133  */
 134 Bool isIntegral (Matrix *A) {
 135
 136   unsigned i, j;
 137   Value divisor, tmp;
 138
 139   value_init(divisor); value_init(tmp);
 140   value_assign(divisor,A->p[A->NbRows-1][A->NbColumns-1]);
 141
 142   for (i = 0; i < A->NbRows ; i++)
 143     for (j = 0; j < A->NbColumns ; j++) {
 144       value_modulus(tmp,A->p[i][j],divisor);
 145       if (value_notzero_p(tmp)) {
 146         value_clear(divisor); value_clear(tmp);
 147           return False;
 148       }
 149     }
 150   value_clear(divisor); value_clear(tmp);
 151   return True ;
 152 } /* isIntegral */
 153
 154 /*
 155  * Check if the matrix 'A' is in Hermite normal form or not.
 156  */
 157 Bool isinHnf(Matrix *A) {
 158
 159   Matrix *temp ;
 160   unsigned i, j ;
 161   Value rem;
 162
 163   value_init(rem);
 164   temp = Homogenise(A,True) ;
 165   for (i = 0; i < temp->NbRows; i++) {
 166     value_assign(rem,temp->p[i][i]);
 167     for (j = 0; j < i ; j++)
 168       if (value_ge(temp->p[i][j],rem)) {
 169         Matrix_Free(temp);
 170         value_clear(rem);
 171         return False ;
 172       }
 173     for (j = i+1; j < temp->NbColumns ; j++)
 174       if (value_notzero_p(temp->p[i][j])) {
 175         Matrix_Free(temp);
 176         value_clear(rem);
 177         return False ;
 178       }
 179   }
 180   Matrix_Free(temp);
 181   value_clear(rem);
 182   return True ;
 183 } /* isinHnf */
 184
 185 /*
 186  * Remove the row 'Rownumber' and place it at the end of the matrix 'X'
 187  */
 188 void PutRowLast (Matrix *X, int Rownumber) {
 189
 190   int i, j ;
 191   Value temp;
 192
 193   if (Rownumber == X->NbRows-1)
 194     return;
 195
 196   value_init(temp);
 197   for (j = 0; j < X->NbColumns; j++) {
 198     value_assign(temp,X->p[Rownumber][j]);
 199     for (i = Rownumber; i < X->NbRows-1; i++)
 200       value_assign(X->p[i][j],X->p[i+1][j]);
 201     value_assign(X->p[i][j],temp);
 202   }
 203   value_clear(temp);
 204   return;
 205 } /* PutRowLast */
 206
 207 /*
 208  * Remove the row 'Rownumber' and place it at the begining of the matrix 'X'
 209  */
 210 void PutRowFirst(Matrix *X, int Rownumber) {
 211
 212   int i, j ;
 213   Value temp;
 214
 215   value_init(temp);
 216   for (j = 0; j < X->NbColumns; j++) {
 217     value_assign(temp,X->p[Rownumber][j]);
 218     for (i = Rownumber; i > 0; i--)
 219       value_assign(X->p[i][j],X->p[i-1][j]);
 220     value_assign(X->p[i][j],temp);
 221   }
 222   value_clear(temp);
 223   return;
 224 } /* PutRowFirst */
 225
 226 /*
 227  * Remove the column 'Columnnumber' and place it at the begining of the matrix
 228  * 'X'
 229  */
 230 void PutColumnFirst (Matrix *X, int Columnnumber) {
 231
 232   int i, j ;
 233   Value temp;
 234
 235   value_init(temp);
 236   for (i = 0; i < X->NbRows; i ++) {
 237     value_assign(temp,X->p[i][Columnnumber]);
 238     for (j = Columnnumber; j > 0; j --)
 239       value_assign(X->p[i][j],X->p[i][j-1]);
 240     value_assign(X->p[i][0],temp);
 241   }
 242   value_clear(temp);
 243   return ;
 244 } /* PutColumnFirst */
 245
 246 /*
 247  * Remove the column 'Columnnumber' and place it at the end of the matrix 'X'
 248  */
 249 void PutColumnLast (Matrix *X, int Columnnumber) {
 250
 251   int i, j ;
 252   Value temp;
 253
 254   value_init(temp);
 255   for (i = 0; i < X->NbRows; i++) {
 256     value_assign(temp,X->p[i][Columnnumber]);
 257     for (j = Columnnumber; j < X->NbColumns-1; j++)
 258       value_assign(X->p[i][j],X->p[i][j+1]);
 259     value_assign(X->p[i][X->NbColumns-1],temp);
 260   }
 261   value_clear(temp);
 262   return ;
 263 } /* PutColumnLast */
 264
 265 /*
 266  * Add a row of zeros at the end of the matrix 'M' and return the new matrix
 267  */
 268 Matrix *AddANullRow (Matrix *M) {
 269
 270   int i,j;
 271   Matrix *Result;
 272
 273   Result = Matrix_Alloc(M->NbRows+1,M->NbColumns);
 274   for (i = 0;i < M->NbRows; i++)
 275     for (j = 0; j < M->NbColumns ; j++)
 276       value_assign(Result->p[i][j],M->p[i][j]);
 277   for (j = 0; j < M->NbColumns; j++)
 278     value_set_si(Result->p[i][j],0);
 279   return Result;
 280 } /* AddANullRow */
 281
 282 /*
 283  * Add a column of zeros at the end of the matrix 'M' and return the new
 284  * matrix
 285  */
 286 Matrix *AddANullColumn(Matrix *M) {
 287
 288   int i,j;
 289   Matrix *Result;
 290
 291   Result = Matrix_Alloc(M->NbRows, M->NbColumns+1);
 292   for (i = 0;i < M->NbRows; i++)
 293     for (j = 0; j < M->NbColumns ; j++)
 294       value_assign(Result->p[i][j],M->p[i][j]);
 295   for (i = 0; i < M->NbRows; i++)
 296     value_set_si(Result->p[i][M->NbColumns],0);
 297   return Result;
 298 } /* AddANullColumn */
 299
 300 /*
 301  * Remove a row 'Rownumber' from matrix 'M' and return the new matrix.
 302  */
 303 Matrix *RemoveRow(Matrix *M, int Rownumber) {
 304
 305   Matrix *Result;
 306   int i;
 307
 308   Result = Matrix_Alloc(M->NbRows-1, M->NbColumns);
 309
 310   for (i = 0; i < Rownumber; i++)
 311     Vector_Copy(M->p[i], Result->p[i], M->NbColumns);
 312   for ( ; i < Result->NbRows; i++)
 313     Vector_Copy(M->p[i+1], Result->p[i], M->NbColumns);
 314
 315   return Result;
 316 } /* RemoveRow */
 317
 318 /*
 319  * Remove NumColumns columns starting at column number 'FirstColumnnumber' from matrix 'M'
 320  * and return the new matrix.
 321  */
 322 Matrix *RemoveNColumns (Matrix *M, int FirstColumnnumber, int NumColumns) {
 323
 324   Matrix *Result;
 325   int i;
 326
 327   Result = Matrix_Alloc (M->NbRows, M->NbColumns-NumColumns);
 328
 329   for (i = 0; i < Result->NbRows; i++) {
 330     Vector_Copy(M->p[i], Result->p[i], FirstColumnnumber);
 331     Vector_Copy(M->p[i]+FirstColumnnumber+NumColumns, Result->p[i]+FirstColumnnumber,
 332                 M->NbColumns-NumColumns-FirstColumnnumber);
 333   }
 334   return Result;
 335 } /* RemoveColumn */
 336
 337 /*
 338  * Remove a column 'Columnnumber' from matrix 'M' and return the new matrix.
 339  */
 340 Matrix *RemoveColumn (Matrix *M, int Columnnumber) {
 341
 342   Matrix *Result;
 343   int i;
 344
 345   Result = Matrix_Alloc (M->NbRows, M->NbColumns-1);
 346
 347   for (i = 0; i < Result->NbRows; i++) {
 348     Vector_Copy(M->p[i], Result->p[i], Columnnumber);
 349     Vector_Copy(M->p[i]+Columnnumber+1, Result->p[i]+Columnnumber,
 350                 M->NbColumns-1-Columnnumber);
 351   }
 352   return Result;
 353 } /* RemoveColumn */
 354
 355 /*
 356  * Given a Matrix M of dimesnion n * l and rank l1, find a unimodular matrix
 357  * 'Result' such that the Vector Space spanned by M is the subset of the vector
 358  * Space spanned by the first l1 Rows of Result. The function returns the rank
 359  * l1 and the Matrix Result.
 360  */
 361 int findHermiteBasis(Matrix *M, Matrix **Result) {
 362
 363   int i, j;
 364   Matrix *C, *curMat, *temp1, *temp2;
 365   Matrix *H, *U;
 366   Vector *V;
 367   int dim, curDim, curVect,rank;
 368
 369   if (M->NbRows == 0) {
 370     Result[0] = Identity (M->NbColumns);
 371     return 0;
 372   }
 373
 374   if (M->NbRows <= M->NbColumns) {
 375     Hermite(M, &H, &U);
 376
 377     for (i = 0; i < H->NbRows; i++)
 378       if (value_zero_p(H->p[i][i]))
 379         break;
 380
 381     if (i == H->NbRows) {
 382       Result[0] = Transpose(U);
 383       Matrix_Free(H);
 384       Matrix_Free(U);
 385       return(i);
 386     }
 387     Matrix_Free (H);
 388     Matrix_Free (U);
 389   }
 390
 391   /* Eliminating the Zero Rows  */
 392
 393   C = Matrix_Copy (M);
 394   for (i = 0; i < C->NbRows; i++) {
 395     for (j = 0; j < C->NbColumns; j++)
 396       if (value_notzero_p(C->p[i][j]))
 397         break;
 398     if (j == C->NbColumns) {
 399       Matrix *temp;
 400       temp = RemoveRow(C, i);
 401       Matrix_Free(C);
 402       C = Matrix_Copy(temp);
 403       Matrix_Free(temp);
 404       i --;
 405     }
 406   }
 407
 408   /* Eliminating the Redundant Rows */
 409
 410   curDim = 1;
 411   curVect = 1;
 412   dim = C->NbColumns;
 413
 414   curMat = Matrix_Alloc(1,C->NbColumns);
 415   for (i = 0; i < C->NbColumns; i ++)
 416     value_assign(curMat->p[0][i],C->p[0][i]);
 417
 418   while((curVect < C->NbRows) && (curDim < dim)) {
 419     Matrix *temp;
 420     temp = AddANullRow(curMat);
 421     for (i = 0; i < C->NbColumns; i++)
 422       value_assign(temp->p[temp->NbRows-1][i],C->p[curVect][i]);
 423
 424     temp1 = AddANullRow(temp);
 425     temp2 = AddANullColumn(temp1);
 426     rank = SolveDiophantine(temp2, &U, &V);
 427     if (rank == temp->NbRows) {
 428       Matrix_Free(curMat);
 429       curMat = Matrix_Copy(temp);
 430       curDim ++;
 431     }
 432     curVect ++;
 433     Matrix_Free (U);
 434     Vector_Free (V);
 435     Matrix_Free (temp1);
 436     Matrix_Free (temp);
 437     Matrix_Free (temp2);
 438   }
 439   Matrix_Free(C);
 440
 441   Hermite(curMat, &H, &U);
 442   rank = curMat->NbRows;
 443   Matrix_Free(curMat);
 444
 445   Result[0] = Transpose (U);
 446   Matrix_Free (H);
 447   Matrix_Free (U);
 448   return (rank);
 449 } /* findHermiteBasis */
 450
 451
 452
 453