source/kernel/Matop.c

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