source/kernel/compress_parms.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 /**
  19  * $Id: compress_parms.c,v 1.32 2006/11/03 17:34:26 skimo Exp $
  20  *
  21  * The integer points in a parametric linear subspace of Q^n are generally
  22  * lying on a sub-lattice of Z^n.
  23  * Functions here use and compute validity lattices, i.e. lattices induced on a
  24  * set of variables by such equalities involving another set of integer
  25  * variables.
  26  * @author B. Meister 12/2003-2006 meister@icps.u-strasbg.fr
  27  * LSIIT -ICPS
  28  * UMR 7005 CNRS
  29  * Louis Pasteur University (ULP), Strasbourg, France
  30 */
  31
  32 #include <stdlib.h>
  33 #include <polylib/polylib.h>
  34
  35 /**
  36  * debug flags (2 levels)
  37  */
  38 #define dbgCompParm 0
  39 #define dbgCompParmMore 0
  40
  41 #define dbgStart(a) if (dbgCompParmMore) { printf(" -- begin "); \
  42                                            printf(#a);        \
  43                                            printf(" --\n"); }   \
  44                                            while(0)
  45 #define dbgEnd(a) if (dbgCompParmMore) { printf(" -- end "); \
  46                                          printf(#a);      \
  47                                          printf(" --\n"); } \
  48                                          while(0)
  49
  50
  51 /**
  52  * Given a full-row-rank nxm matrix M made of m row-vectors), computes the
  53  * basis K (made of n-m column-vectors) of the integer kernel of the rows of M
  54  * so we have: M.K = 0
  55 */
  56 Matrix * int_ker(Matrix * M) {
  57   Matrix *U, *Q, *H, *H2, *K=NULL;
  58   int i, j, rk;
  59
  60   if (dbgCompParm)
  61     show_matrix(M);
  62   /* eliminate redundant rows : UM = H*/
  63   right_hermite(M, &H, &Q, &U);
  64   for (rk=H->NbRows-1; (rk>=0) && Vector_IsZero(H->p[rk], H->NbColumns); rk--);
  65   rk++;
  66   if (dbgCompParmMore) {
  67     printf("rank = %d\n", rk);
  68   }
  69
  70   /* there is a non-null kernel if and only if the dimension m of
  71      the space spanned by the rows
  72      is inferior to the number n of variables */
  73   if (M->NbColumns <= rk) {
  74     Matrix_Free(H);
  75     Matrix_Free(Q);
  76     Matrix_Free(U);
  77     K = Matrix_Alloc(M->NbColumns, 0);
  78     return K;
  79   }
  80   Matrix_Free(U);
  81   Matrix_Free(Q);
  82   /* fool left_hermite  by giving NbRows =rank of M*/
  83   H->NbRows=rk;
  84   /* computes MU = [H 0] */
  85   left_hermite(H, &H2, &Q, &U);
  86    if (dbgCompParmMore) {
  87     printf("-- Int. Kernel -- \n");
  88     show_matrix(M);
  89     printf(" = \n");
  90     show_matrix(H2);
  91     show_matrix(U);
  92   }
  93   H->NbRows==M->NbRows;
  94   Matrix_Free(H);
  95   /* the Integer Kernel is made of the last n-rk columns of U */
  96   Matrix_subMatrix(U, 0, rk, U->NbRows, U->NbColumns, &K);
  97
  98   /* clean up */
  99   Matrix_Free(H2);
 100   Matrix_Free(U);
 101   Matrix_Free(Q);
 102   return K;
 103 } /* int_ker */
 104
 105
 106 /**
 107  * Computes the intersection of two linear lattices, whose base vectors are
 108  * respectively represented in A and B.
 109  * If I and/or Lb is set to NULL, then the matrix is allocated.
 110  * Else, the matrix is assumed to be allocated already.
 111  * I and Lb are rk x rk, where rk is the rank of A (or B).
 112  * @param A the full-row rank matrix whose column-vectors are the basis for the
 113  * first linear lattice.
 114  * @param B the matrix whose column-vectors are the basis for the second linear
 115  * lattice.
 116  * @param Lb the matrix such that B.Lb = I, where I is the intersection.
 117  * @return their intersection.
 118  */
 119 static void linearInter(Matrix * A, Matrix * B, Matrix ** I, Matrix **Lb) {
 120   Matrix * AB=NULL;
 121   int rk = A->NbRows;
 122   int a = A->NbColumns;
 123   int b = B->NbColumns;
 124   int i,j, z=0;
 125
 126   Matrix * H, *U, *Q;
 127   /* ensure that the spanning vectors are in the same space */
 128   assert(B->NbRows==rk);
 129   /* 1- build the matrix
 130    * (A 0 1)
 131    * (0 B 1)
 132    */
 133   AB = Matrix_Alloc(2*rk, a+b+rk);
 134   Matrix_copySubMatrix(A, 0, 0, rk, a, AB, 0, 0);
 135   Matrix_copySubMatrix(B, 0, 0, rk, b, AB, rk, a);
 136   for (i=0; i< rk; i++) {
 137       value_set_si(AB->p[i][a+b+i], 1);
 138       value_set_si(AB->p[i+rk][a+b+i], 1);
 139   }
 140   if (dbgCompParm) {
 141     show_matrix(AB);
 142   }
 143
 144   /* 2- Compute its left Hermite normal form. AB.U = [H 0] */
 145   left_hermite(AB, &H, &Q, &U);
 146   Matrix_Free(AB);
 147   Matrix_Free(Q);
 148   /* count the number of non-zero colums in H */
 149   for (z=H->NbColumns-1; value_zero_p(H->p[H->NbRows-1][z]); z--);
 150   z++;
 151   if (dbgCompParm) {
 152     show_matrix(H);
 153     printf("z=%d\n", z);
 154   }
 155   Matrix_Free(H);
 156   /* if you split U in 9 submatrices, you have:
 157    * A.U_13 = -U_33
 158    * B.U_23 = -U_33,
 159    * where the nb of cols of U_{*3} equals the nb of zero-cols of H
 160    * U_33 is a (the smallest) combination of col-vectors of A and B at the same
 161    * time: their intersection.
 162   */
 163   Matrix_subMatrix(U, a+b, z, U->NbColumns, U->NbColumns, I);
 164   Matrix_subMatrix(U, a, z, a+b, U->NbColumns, Lb);
 165   if (dbgCompParm) {
 166     show_matrix(U);
 167   }
 168   Matrix_Free(U);
 169 } /* linearInter */
 170
 171
 172 /**
 173  * Given a system of equalities, looks if it has an integer solution in the
 174  * combined space, and if yes, returns one solution.
 175  * <p>pre-condition: the equalities are full-row rank (without the constant
 176  * part)</p>
 177  * @param Eqs the system of equations (as constraints)
 178  * @param I a feasible integer solution if it exists, else NULL. Allocated if
 179  * initially set to NULL, else reused.
 180  */
 181 void Equalities_integerSolution(Matrix * Eqs, Matrix **I) {
 182   Matrix * Hm, *H=NULL, *U, *Q, *M=NULL, *C=NULL, *Hi;
 183   Matrix *Ip;
 184   int i;
 185   Value mod;
 186   unsigned int rk;
 187   if (Eqs==NULL){
 188     if ((*I)!=NULL) Matrix_Free(*I);
 189     I = NULL;
 190     return;
 191   }
 192   /* we use: AI = C = (Ha 0).Q.I = (Ha 0)(I' 0)^T */
 193   /* with I = Qinv.I' = U.I'*/
 194   /* 1- compute I' = Hainv.(-C) */
 195   /* HYP: the equalities are full-row rank */
 196   rk = Eqs->NbRows;
 197   Matrix_subMatrix(Eqs, 0, 1, rk, Eqs->NbColumns-1, &M);
 198   left_hermite(M, &Hm, &Q, &U);
 199   Matrix_Free(M);
 200   Matrix_subMatrix(Hm, 0, 0, rk, rk, &H);
 201   if (dbgCompParmMore) {
 202     show_matrix(Hm);
 203     show_matrix(H);
 204     show_matrix(U);
 205   }
 206   Matrix_Free(Q);
 207   Matrix_Free(Hm);
 208   Matrix_subMatrix(Eqs, 0, Eqs->NbColumns-1, rk, Eqs->NbColumns, &C);
 209   Matrix_oppose(C);
 210   Hi = Matrix_Alloc(rk, rk+1);
 211   MatInverse(H, Hi);
 212   if (dbgCompParmMore) {
 213     show_matrix(C);
 214     show_matrix(Hi);
 215   }
 216   /* put the numerator of Hinv back into H */
 217   Matrix_subMatrix(Hi, 0, 0, rk, rk, &H);
 218   Ip = Matrix_Alloc(Eqs->NbColumns-2, 1);
 219   /* fool Matrix_Product on the size of Ip */
 220   Ip->NbRows = rk;
 221   Matrix_Product(H, C, Ip);
 222   Ip->NbRows = Eqs->NbColumns-2;
 223   Matrix_Free(H);
 224   Matrix_Free(C);
 225   value_init(mod);
 226   for (i=0; i< rk; i++) {
 227     /* if Hinv.C is not integer, return NULL (no solution) */
 228     value_pmodulus(mod, Ip->p[i][0], Hi->p[i][rk]);
 229     if (value_notzero_p(mod)) {
 230       if ((*I)!=NULL) Matrix_Free(*I);
 231       value_clear(mod);
 232       Matrix_Free(U);
 233       Matrix_Free(Ip);
 234       Matrix_Free(Hi);
 235       I = NULL;
 236       return;
 237     }
 238     else {
 239       value_pdivision(Ip->p[i][0], Ip->p[i][0], Hi->p[i][rk]);
 240     }
 241   }
 242   /* fill the rest of I' with zeros */
 243   for (i=rk; i< Eqs->NbColumns-2; i++) {
 244     value_set_si(Ip->p[i][0], 0);
 245   }
 246   value_clear(mod);
 247   Matrix_Free(Hi);
 248   /* 2 - Compute the particular solution I = U.(I' 0) */
 249   ensureMatrix((*I), Eqs->NbColumns-2, 1);
 250   Matrix_Product(U, Ip, (*I));
 251   Matrix_Free(U);
 252   Matrix_Free(Ip);
 253   if (dbgCompParm) {
 254     show_matrix(*I);
 255   }
 256 }
 257
 258
 259 /**
 260  * Computes the validity lattice of a set of equalities. I.e., the lattice
 261  * induced on the last <tt>b</tt> variables by the equalities involving the
 262  * first <tt>a</tt> integer existential variables.  The submatrix of Eqs that
 263  * concerns only the existential variables (so the first a columns) is assumed
 264  * to be full-row rank.
 265  * @param Eqs the equalities
 266  * @param a the number of existential integer variables, placed as first
 267  * variables
 268  * @param vl the (returned) validity lattice, in homogeneous form. It is
 269  * allocated if initially set to null, or reused if already allocated.
 270  */
 271 void Equalities_validityLattice(Matrix * Eqs, int a, Matrix** vl) {
 272   unsigned int b = Eqs->NbColumns-2-a;
 273   unsigned int r = Eqs->NbRows;
 274   Matrix * A=NULL, * B=NULL, *I = NULL, *Lb=NULL, *sol=NULL;
 275   Matrix *H, *U, *Q;
 276   unsigned int i;
 277
 278   if (dbgCompParm) {
 279     printf("Computing validity lattice induced by the %d first variables of:"
 280            ,a);
 281     show_matrix(Eqs);
 282   }
 283   if (b==0) {
 284     ensureMatrix((*vl), 1, 1);
 285     value_set_si((*vl)->p[0][0], 1);
 286     return;
 287   }
 288
 289   /* 1- check that there is an integer solution to the equalities */
 290   /* OPT: could change integerSolution's profile to allocate or not*/
 291   Equalities_integerSolution(Eqs, &sol);
 292   /* if there is no integer solution, there is no validity lattice */
 293   if (sol==NULL) {
 294     if ((*vl)!=NULL) Matrix_Free(*vl);
 295     return;
 296   }
 297   Matrix_subMatrix(Eqs, 0, 1, r, 1+a, &A);
 298   Matrix_subMatrix(Eqs, 0, 1+a, r, 1+a+b, &B);
 299   linearInter(A, B, &I, &Lb);
 300   Matrix_Free(A);
 301   Matrix_Free(B);
 302   Matrix_Free(I);
 303   if (dbgCompParm) {
 304     show_matrix(Lb);
 305   }
 306
 307   /* 2- The linear part of the validity lattice is the left HNF of Lb */
 308   left_hermite(Lb, &H, &Q, &U);
 309   Matrix_Free(Lb);
 310   Matrix_Free(Q);
 311   Matrix_Free(U);
 312
 313   /* 3- build the validity lattice */
 314   ensureMatrix((*vl), b+1, b+1);
 315   Matrix_copySubMatrix(H, 0, 0, b, b, (*vl), 0,0);
 316   Matrix_Free(H);
 317   for (i=0; i< b; i++) {
 318     value_assign((*vl)->p[i][b], sol->p[0][a+i]);
 319   }
 320   Matrix_Free(sol);
 321   Vector_Set((*vl)->p[b],0, b);
 322   value_set_si((*vl)->p[b][b], 1);
 323
 324 } /* validityLattice */
 325
 326
 327 /**
 328  * Eliminate the columns corresponding to a list of eliminated parameters.
 329  * @param M the constraints matrix whose columns are to be removed
 330  * @param nbVars an offset to be added to the ranks of the variables to be
 331  * removed
 332  * @param elimParms the list of ranks of the variables to be removed
 333  * @param newM (output) the matrix without the removed columns
 334  */
 335 void Constraints_removeElimCols(Matrix * M, unsigned int nbVars,
 336                            unsigned int *elimParms, Matrix ** newM) {
 337   unsigned int i, j, k;
 338   if (elimParms[0]==0) {
 339     Matrix_clone(M, newM);
 340     return;
 341   }
 342   if ((*newM)==NULL) {
 343     (*newM) = Matrix_Alloc(M->NbRows, M->NbColumns - elimParms[0]);
 344   }
 345   else {
 346     assert ((*newM)->NbColumns==M->NbColumns - elimParms[0]);
 347   }
 348   for (i=0; i< M->NbRows; i++) {
 349     value_assign((*newM)->p[i][0], M->p[i][0]); /* kind of cstr */
 350     k=0;
 351     Vector_Copy(&(M->p[i][1]), &((*newM)->p[i][1]), nbVars);
 352     for (j=0; j< M->NbColumns-2-nbVars; j++) {
 353       if (j!=elimParms[k+1]) {
 354         value_assign((*newM)->p[i][j-k+nbVars+1], M->p[i][j+nbVars+1]);
 355       }
 356       else {
 357         k++;
 358       }
 359     }
 360     value_assign((*newM)->p[i][(*newM)->NbColumns-1],
 361                  M->p[i][M->NbColumns-1]); /* cst part */
 362   }
 363 } /* Constraints_removeElimCols */
 364
 365
 366 /**
 367  * Eliminates all the equalities in a set of constraints and returns the set of
 368  * constraints defining a full-dimensional polyhedron, such that there is a
 369  * bijection between integer points of the original polyhedron and these of the
 370  * resulting (projected) polyhedron).
 371  * If VL is set to NULL, this funciton allocates it. Else, it assumes that
 372  * (*VL) points to a matrix of the right size.
 373  * <p> The following things are done:
 374  * <ol>
 375  * <li> remove equalities involving only parameters, and remove as many
 376  *      parameters as there are such equalities. From that, the list of
 377  *      eliminated parameters <i>elimParms</i> is built.
 378  * <li> remove equalities that involve variables. This requires a compression
 379  *      of the parameters and of the other variables that are not eliminated.
 380  *      The affine compresson is represented by matrix VL (for <i>validity
 381  *      lattice</i>) and is such that (N I 1)^T = VL.(N' I' 1), where N', I'
 382  *      are integer (they are the parameters and variables after compression).
 383  *</ol>
 384  *</p>
 385  */
 386 void Constraints_fullDimensionize(Matrix ** M, Matrix ** C, Matrix ** VL,
 387                                   Matrix ** Eqs, Matrix ** ParmEqs,
 388                                   unsigned int ** elimVars,
 389                                   unsigned int ** elimParms,
 390                                   int maxRays) {
 391   unsigned int i, j;
 392   Matrix * A=NULL, *B=NULL;
 393   Matrix * Ineqs=NULL;
 394   unsigned int nbVars = (*M)->NbColumns - (*C)->NbColumns;
 395   unsigned int nbParms;
 396   int nbElimVars;
 397   Matrix * fullDim = NULL;
 398
 399   /* variables for permutations */
 400   unsigned int * permutation;
 401   Matrix * permutedEqs=NULL, * permutedIneqs=NULL;
 402
 403   /* 1- Eliminate the equalities involving only parameters. */
 404   (*ParmEqs) = Constraints_removeParmEqs(M, C, 0, elimParms);
 405   /* if the polyehdron is empty, return now. */
 406   if ((*M)->NbColumns==0) return;
 407   /* eliminate the columns corresponding to the eliminated parameters */
 408   if (elimParms[0]!=0) {
 409     Constraints_removeElimCols(*M, nbVars, (*elimParms), &A);
 410     Matrix_Free(*M);
 411     (*M) = A;
 412     Constraints_removeElimCols(*C, 0, (*elimParms), &B);
 413     Matrix_Free(*C);
 414     (*C) = B;
 415     if (dbgCompParm) {
 416       printf("After false parameter elimination: \n");
 417       show_matrix(*M);
 418       show_matrix(*C);
 419     }
 420   }
 421   nbParms = (*C)->NbColumns-2;
 422
 423   /* 2- Eliminate the equalities involving variables */
 424   /*   a- extract the (remaining) equalities from the poyhedron */
 425   split_constraints((*M), Eqs, &Ineqs);
 426   nbElimVars = (*Eqs)->NbRows;
 427   /*    if the polyhedron is already full-dimensional, return */
 428   if ((*Eqs)->NbRows==0) {
 429     Matrix_identity(nbParms+1, VL);
 430     return;
 431   }
 432   /*   b- choose variables to be eliminated */
 433   permutation = find_a_permutation((*Eqs), nbParms);
 434
 435   if (dbgCompParm) {
 436     printf("Permuting the vars/parms this way: [ ");
 437     for (i=0; i< (*Eqs)->NbColumns-2; i++) {
 438       printf("%d ", permutation[i]);
 439     }
 440     printf("]\n");
 441   }
 442
 443   Constraints_permute((*Eqs), permutation, &permutedEqs);
 444   Equalities_validityLattice(permutedEqs, (*Eqs)->NbRows, VL);
 445
 446   if (dbgCompParm) {
 447     printf("Validity lattice: ");
 448     show_matrix(*VL);
 449   }
 450   Constraints_compressLastVars(permutedEqs, (*VL));
 451   Constraints_permute(Ineqs, permutation, &permutedIneqs);
 452   if (dbgCompParmMore) {
 453     show_matrix(permutedIneqs);
 454     show_matrix(permutedEqs);
 455   }
 456   Matrix_Free(*Eqs);
 457   Matrix_Free(Ineqs);
 458   Constraints_compressLastVars(permutedIneqs, (*VL));
 459   if (dbgCompParm) {
 460     printf("After compression: ");
 461     show_matrix(permutedIneqs);
 462   }
 463   /*   c- eliminate the first variables */
 464   assert(Constraints_eliminateFirstVars(permutedEqs, permutedIneqs));
 465   if (dbgCompParmMore) {
 466     printf("After elimination of the variables: ");
 467     show_matrix(permutedIneqs);
 468   }
 469
 470   /*   d- get rid of the first (zero) columns,
 471        which are now useless, and put the parameters back at the end */
 472   fullDim = Matrix_Alloc(permutedIneqs->NbRows,
 473                          permutedIneqs->NbColumns-nbElimVars);
 474   for (i=0; i< permutedIneqs->NbRows; i++) {
 475     value_set_si(fullDim->p[i][0], 1);
 476     for (j=0; j< nbParms; j++) {
 477       value_assign(fullDim->p[i][j+fullDim->NbColumns-nbParms-1],
 478                    permutedIneqs->p[i][j+nbElimVars+1]);
 479     }
 480     for (j=0; j< permutedIneqs->NbColumns-nbParms-2-nbElimVars; j++) {
 481       value_assign(fullDim->p[i][j+1],
 482                    permutedIneqs->p[i][nbElimVars+nbParms+j+1]);
 483     }
 484     value_assign(fullDim->p[i][fullDim->NbColumns-1],
 485                  permutedIneqs->p[i][permutedIneqs->NbColumns-1]);
 486   }
 487   Matrix_Free(permutedIneqs);
 488
 489 } /* Constraints_fullDimensionize */
 490
 491
 492 /**
 493  * Given a matrix that defines a full-dimensional affine lattice, returns the
 494  * affine sub-lattice spanned in the k first dimensions.
 495  * Useful for instance when you only look for the parameters' validity lattice.
 496  * @param lat the original full-dimensional lattice
 497  * @param subLat the sublattice
 498  */
 499 void Lattice_extractSubLattice(Matrix * lat, unsigned int k, Matrix ** subLat) {
 500   Matrix * H, *Q, *U, *linLat = NULL;
 501   unsigned int i;
 502   dbgStart(Lattice_extractSubLattice);
 503   /* if the dimension is already good, just copy the initial lattice */
 504   if (k==lat->NbRows-1) {
 505     if (*subLat==NULL) {
 506       (*subLat) = Matrix_Copy(lat);
 507     }
 508     else {
 509       Matrix_copySubMatrix(lat, 0, 0, lat->NbRows, lat->NbColumns, (*subLat), 0, 0);
 510     }
 511     return;
 512   }
 513   assert(k<lat->NbRows-1);
 514   /* 1- Make the linear part of the lattice triangular to eliminate terms from
 515      other dimensions */
 516   Matrix_subMatrix(lat, 0, 0, lat->NbRows, lat->NbColumns-1, &linLat);
 517   /* OPT: any integer column-vector elimination is ok indeed. */
 518   /* OPT: could test if the lattice is already in triangular form. */
 519   left_hermite(linLat, &H, &Q, &U);
 520   if (dbgCompParmMore) {
 521     show_matrix(H);
 522   }
 523   Matrix_Free(Q);
 524   Matrix_Free(U);
 525   Matrix_Free(linLat);
 526   /* if not allocated yet, allocate it */
 527   if (*subLat==NULL) {
 528     (*subLat) = Matrix_Alloc(k+1, k+1);
 529   }
 530   Matrix_copySubMatrix(H, 0, 0, k, k, (*subLat), 0, 0);
 531   Matrix_Free(H);
 532   Matrix_copySubMatrix(lat, 0, lat->NbColumns-1, k, 1, (*subLat), 0, k);
 533   for (i=0; i<k; i++) {
 534     value_set_si((*subLat)->p[k][i], 0);
 535   }
 536   value_set_si((*subLat)->p[k][k], 1);
 537   dbgEnd(Lattice_extractSubLattice);
 538 } /* Lattice_extractSubLattice */
 539
 540
 541 /**
 542  * Computes the overall period of the variables I for (MI) mod |d|, where M is
 543  * a matrix and |d| a vector. Produce a diagonal matrix S = (s_k) where s_k is
 544  * the overall period of i_k
 545  * @param M the set of affine functions of I (row-vectors)
 546  * @param d the column-vector representing the modulos
 547 */
 548 Matrix * affine_periods(Matrix * M, Matrix * d) {
 549   Matrix * S;
 550   unsigned int i,j;
 551   Value tmp;
 552   Value * periods = (Value *)malloc(sizeof(Value) * M->NbColumns);
 553   value_init(tmp);
 554   for(i=0; i< M->NbColumns; i++) {
 555     value_init(periods[i]);
 556     value_set_si(periods[i], 1);
 557   }
 558   for (i=0; i<M->NbRows; i++) {
 559     for (j=0; j< M->NbColumns; j++) {
 560       value_gcd(tmp, d->p[i][0], M->p[i][j]);
 561       value_divexact(tmp, d->p[i][0], tmp);
 562       value_lcm(periods[j], periods[j], tmp);
 563      }
 564   }
 565   value_clear(tmp);
 566
 567   /* 2- build S */
 568   S = Matrix_Alloc(M->NbColumns, M->NbColumns);
 569   for (i=0; i< M->NbColumns; i++)
 570     for (j=0; j< M->NbColumns; j++)
 571       if (i==j) value_assign(S->p[i][j],periods[j]);
 572       else value_set_si(S->p[i][j], 0);
 573
 574   /* 3- clean up */
 575   for(i=0; i< M->NbColumns; i++) value_clear(periods[i]);
 576   free(periods);
 577   return S;
 578 } /* affine_periods */
 579
 580
 581 /**
 582  * Given an integer matrix B with m rows and integer m-vectors C and d,
 583  * computes the basis of the integer solutions to (BN+C) mod d = 0 (1).
 584  * This is an affine lattice (G): (N 1)^T= G(N' 1)^T, forall N' in Z^b.
 585  * If there is no solution, returns NULL.
 586  * @param B B, a (m x b) matrix
 587  * @param C C, a (m x 1) integer matrix
 588  * @param d d, a (1 x m) integer matrix
 589  * @param imb the affine (b+1)x(b+1) basis of solutions, in the homogeneous
 590  * form. Allocated if initially set to NULL, reused if not.
 591 */
 592 void Equalities_intModBasis(Matrix * B, Matrix * C, Matrix * d, Matrix ** imb) {
 593   int b = B->NbColumns;
 594   /* FIXME: treat the case d=0 as a regular equality B_kN+C_k = 0: */
 595   /* OPT: could keep only equalities for which d>1 */
 596   int nbEqs = B->NbRows;
 597   unsigned int i;
 598
 599   /* 1- buid the problem DI+BN+C = 0 */
 600   Matrix * eqs = Matrix_Alloc(nbEqs, nbEqs+b+1);
 601   for (i=0; i< nbEqs; i++) {
 602     value_assign(eqs->p[i][i], d->p[0][i]);
 603   }
 604   Matrix_copySubMatrix(B, 0, 0, nbEqs, b, eqs, 0, nbEqs);
 605   Matrix_copySubMatrix(C, 0, 0, nbEqs, 1, eqs, 0, nbEqs+b);
 606
 607   /* 2- the solution is the validity lattice of the equalities */
 608   Equalities_validityLattice(eqs, nbEqs, imb);
 609   Matrix_Free(eqs);
 610 } /* Equalities_intModBasis */
 611
 612
 613 /** kept here for backwards compatiblity. Wrapper to Equalities_intModBasis() */
 614 Matrix * int_mod_basis(Matrix * B, Matrix * C, Matrix * d) {
 615   Matrix * imb = NULL;
 616   Equalities_intModBasis(B, C, d, &imb);
 617   return imb;
 618 } /* int_mod_basis */
 619
 620
 621 /**
 622  * Given a parameterized constraints matrix with m equalities, computes the
 623  * compression matrix G such that there is an integer solution in the variables
 624  * space for each value of N', with N = G N' (N are the "nb_parms" parameters)
 625  * @param E a matrix of parametric equalities @param nb_parms the number of
 626  * parameters
 627  * <b>Note: </b>this function is mostly here for backwards
 628  * compatibility. Prefer the use of <tt>Equalities_validityLattice</tt>.
 629 */
 630 Matrix * compress_parms(Matrix * E, int nbParms) {
 631   Matrix * vl=NULL;
 632   Equalities_validityLattice(E, E->NbColumns-2-nbParms, &vl);
 633   return vl;
 634 }/* compress_parms */
 635
 636
 637 /** Removes the equalities that involve only parameters, by eliminating some
 638  * parameters in the polyhedron's constraints and in the context.<p>
 639  * <b>Updates M and Ctxt.</b>
 640  * @param M1 the polyhedron's constraints
 641  * @param Ctxt1 the constraints of the polyhedron's context
 642  * @param renderSpace tells if the returned equalities must be expressed in the
 643  * parameters space (renderSpace=0) or in the combined var/parms space
 644  * (renderSpace = 1)
 645  * @param elimParms the list of parameters that have been removed: an array
 646  * whose 1st element is the number of elements in the list.  (returned)
 647  * @return the system of equalities that involve only parameters.
 648  */
 649 Matrix * Constraints_Remove_parm_eqs(Matrix ** M1, Matrix ** Ctxt1,
 650                                      int renderSpace,
 651                                      unsigned int ** elimParms) {
 652   int i, j, k, nbEqsParms =0;
 653   int nbEqsM, nbEqsCtxt, allZeros, nbTautoM = 0, nbTautoCtxt = 0;
 654   Matrix * M = (*M1);
 655   Matrix * Ctxt = (*Ctxt1);
 656   int nbVars = M->NbColumns-Ctxt->NbColumns;
 657   Matrix * Eqs;
 658   Matrix * EqsMTmp;
 659
 660   /* 1- build the equality matrix(ces) */
 661   nbEqsM = 0;
 662   for (i=0; i< M->NbRows; i++) {
 663     k = First_Non_Zero(M->p[i], M->NbColumns);
 664     /* if it is a tautology, count it as such */
 665     if (k==-1) {
 666       nbTautoM++;
 667     }
 668     else {
 669       /* if it only involves parameters, count it */
 670       if (k>= nbVars+1) nbEqsM++;
 671     }
 672   }
 673
 674   nbEqsCtxt = 0;
 675   for (i=0; i< Ctxt->NbRows; i++) {
 676     if (value_zero_p(Ctxt->p[i][0])) {
 677       if (First_Non_Zero(Ctxt->p[i], Ctxt->NbColumns)==-1) {
 678         nbTautoCtxt++;
 679       }
 680       else {
 681         nbEqsCtxt ++;
 682       }
 683     }
 684   }
 685   nbEqsParms = nbEqsM + nbEqsCtxt;
 686
 687   /* nothing to do in this case */
 688   if (nbEqsParms+nbTautoM+nbTautoCtxt==0) {
 689     (*elimParms) = (unsigned int*) malloc(sizeof(int));
 690     (*elimParms)[0] = 0;
 691     if (renderSpace==0) {
 692       return Matrix_Alloc(0,Ctxt->NbColumns);
 693     }
 694     else {
 695       return Matrix_Alloc(0,M->NbColumns);
 696     }
 697   }
 698
 699   Eqs= Matrix_Alloc(nbEqsParms, Ctxt->NbColumns);
 700   EqsMTmp= Matrix_Alloc(nbEqsParms, M->NbColumns);
 701
 702   /* copy equalities from the context */
 703   k = 0;
 704   for (i=0; i< Ctxt->NbRows; i++) {
 705     if (value_zero_p(Ctxt->p[i][0])
 706                      && First_Non_Zero(Ctxt->p[i], Ctxt->NbColumns)!=-1) {
 707       Vector_Copy(Ctxt->p[i], Eqs->p[k], Ctxt->NbColumns);
 708       Vector_Copy(Ctxt->p[i]+1, EqsMTmp->p[k]+nbVars+1,
 709                   Ctxt->NbColumns-1);
 710       k++;
 711     }
 712   }
 713   for (i=0; i< M->NbRows; i++) {
 714     j=First_Non_Zero(M->p[i], M->NbColumns);
 715     /* copy equalities that involve only parameters from M */
 716     if (j>=nbVars+1) {
 717       Vector_Copy(M->p[i]+nbVars+1, Eqs->p[k]+1, Ctxt->NbColumns-1);
 718       Vector_Copy(M->p[i]+nbVars+1, EqsMTmp->p[k]+nbVars+1,
 719                   Ctxt->NbColumns-1);
 720       /* mark these equalities for removal */
 721       value_set_si(M->p[i][0], 2);
 722       k++;
 723     }
 724     /* mark the all-zero equalities for removal */
 725     if (j==-1) {
 726       value_set_si(M->p[i][0], 2);
 727     }
 728   }
 729
 730   /* 2- eliminate parameters until all equalities are used or until we find a
 731   contradiction (overconstrained system) */
 732   (*elimParms) = (unsigned int *) malloc((Eqs->NbRows+1) * sizeof(int));
 733   (*elimParms)[0] = 0;
 734   allZeros = 0;
 735   for (i=0; i< Eqs->NbRows; i++) {
 736     /* find a variable that can be eliminated */
 737     k = First_Non_Zero(Eqs->p[i], Eqs->NbColumns);
 738     if (k!=-1) { /* nothing special to do for tautologies */
 739
 740       /* if there is a contradiction, return empty matrices */
 741       if (k==Eqs->NbColumns-1) {
 742         printf("Contradiction in %dth row of Eqs: ",k);
 743         show_matrix(Eqs);
 744         Matrix_Free(Eqs);
 745         Matrix_Free(EqsMTmp);
 746         (*M1) = Matrix_Alloc(0, M->NbColumns);
 747         Matrix_Free(M);
 748         (*Ctxt1) = Matrix_Alloc(0,Ctxt->NbColumns);
 749         Matrix_Free(Ctxt);
 750         free(*elimParms);
 751         (*elimParms) = (unsigned int *) malloc(sizeof(int));
 752         (*elimParms)[0] = 0;
 753         if (renderSpace==1) {
 754           return Matrix_Alloc(0,(*M1)->NbColumns);
 755         }
 756         else {
 757           return Matrix_Alloc(0,(*Ctxt1)->NbColumns);
 758         }
 759       }
 760       /* if we have something we can eliminate, do it in 3 places:
 761          Eqs, Ctxt, and M */
 762       else {
 763         k--; /* k is the rank of the variable, now */
 764         (*elimParms)[0]++;
 765         (*elimParms)[(*elimParms[0])]=k;
 766         for (j=0; j< Eqs->NbRows; j++) {
 767           if (i!=j) {
 768             eliminate_var_with_constr(Eqs, i, Eqs, j, k);
 769             eliminate_var_with_constr(EqsMTmp, i, EqsMTmp, j, k+nbVars);
 770           }
 771         }
 772         for (j=0; j< Ctxt->NbRows; j++) {
 773           if (value_notzero_p(Ctxt->p[i][0])) {
 774             eliminate_var_with_constr(Eqs, i, Ctxt, j, k);
 775           }
 776         }
 777         for (j=0; j< M->NbRows; j++) {
 778           if (value_cmp_si(M->p[i][0], 2)) {
 779             eliminate_var_with_constr(EqsMTmp, i, M, j, k+nbVars);
 780           }
 781         }
 782       }
 783     }
 784     /* if (k==-1): count the tautologies in Eqs to remove them later */
 785     else {
 786       allZeros++;
 787     }
 788   }
 789
 790   /* elimParms may have been overallocated. Now we know how many parms have
 791      been eliminated so we can reallocate the right amount of memory. */
 792   if (!realloc((*elimParms), ((*elimParms)[0]+1)*sizeof(int))) {
 793     fprintf(stderr, "Constraints_Remove_parm_eqs > cannot realloc()");
 794   }
 795
 796   Matrix_Free(EqsMTmp);
 797
 798   /* 3- remove the "bad" equalities from the input matrices
 799      and copy the equalities involving only parameters */
 800   EqsMTmp = Matrix_Alloc(M->NbRows-nbEqsM-nbTautoM, M->NbColumns);
 801   k=0;
 802   for (i=0; i< M->NbRows; i++) {
 803     if (value_cmp_si(M->p[i][0], 2)) {
 804       Vector_Copy(M->p[i], EqsMTmp->p[k], M->NbColumns);
 805       k++;
 806     }
 807   }
 808   Matrix_Free(M);
 809   (*M1) = EqsMTmp;
 810
 811   EqsMTmp = Matrix_Alloc(Ctxt->NbRows-nbEqsCtxt-nbTautoCtxt, Ctxt->NbColumns);
 812   k=0;
 813   for (i=0; i< Ctxt->NbRows; i++) {
 814     if (value_notzero_p(Ctxt->p[i][0])) {
 815       Vector_Copy(Ctxt->p[i], EqsMTmp->p[k], Ctxt->NbColumns);
 816       k++;
 817     }
 818   }
 819   Matrix_Free(Ctxt);
 820   (*Ctxt1) = EqsMTmp;
 821
 822   if (renderSpace==0) {/* renderSpace=0: equalities in the parameter space */
 823     EqsMTmp = Matrix_Alloc(Eqs->NbRows-allZeros, Eqs->NbColumns);
 824     k=0;
 825     for (i=0; i<Eqs->NbRows; i++) {
 826       if (First_Non_Zero(Eqs->p[i], Eqs->NbColumns)!=-1) {
 827         Vector_Copy(Eqs->p[i], EqsMTmp->p[k], Eqs->NbColumns);
 828         k++;
 829       }
 830     }
 831   }
 832   else {/* renderSpace=1: equalities rendered in the combined space */
 833     EqsMTmp = Matrix_Alloc(Eqs->NbRows-allZeros, (*M1)->NbColumns);
 834     k=0;
 835     for (i=0; i<Eqs->NbRows; i++) {
 836       if (First_Non_Zero(Eqs->p[i], Eqs->NbColumns)!=-1) {
 837         Vector_Copy(Eqs->p[i], &(EqsMTmp->p[k][nbVars]), Eqs->NbColumns);
 838         k++;
 839       }
 840     }
 841   }
 842   Matrix_Free(Eqs);
 843   Eqs = EqsMTmp;
 844
 845   return Eqs;
 846 } /* Constraints_Remove_parm_eqs */
 847
 848
 849 /** Removes equalities involving only parameters, but starting from a
 850  * Polyhedron and its context.
 851  * @param P the polyhedron
 852  * @param C P's context
 853  * @param renderSpace: 0 for the parameter space, =1 for the combined space.
 854  * @maxRays Polylib's usual <i>workspace</i>.
 855  */
 856 Polyhedron * Polyhedron_Remove_parm_eqs(Polyhedron ** P, Polyhedron ** C,
 857                                         int renderSpace,
 858                                         unsigned int ** elimParms,
 859                                         int maxRays) {
 860   Matrix * Eqs;
 861   Polyhedron * Peqs;
 862   Matrix * M = Polyhedron2Constraints((*P));
 863   Matrix * Ct = Polyhedron2Constraints((*C));
 864
 865   /* if the Minkowski representation is not computed yet, do not compute it in
 866      Constraints2Polyhedron */
 867   if (F_ISSET((*P), POL_VALID | POL_INEQUALITIES) &&
 868       (F_ISSET((*C), POL_VALID | POL_INEQUALITIES))) {
 869     FL_INIT(maxRays, POL_NO_DUAL);
 870   }
 871
 872   Eqs = Constraints_Remove_parm_eqs(&M, &Ct, renderSpace, elimParms);
 873   Peqs = Constraints2Polyhedron(Eqs, maxRays);
 874   Matrix_Free(Eqs);
 875
 876   /* particular case: no equality involving only parms is found */
 877   if (Eqs->NbRows==0) {
 878     Matrix_Free(M);
 879     Matrix_Free(Ct);
 880     return Peqs;
 881   }
 882   Polyhedron_Free(*P);
 883   Polyhedron_Free(*C);
 884   (*P) = Constraints2Polyhedron(M, maxRays);
 885   (*C) = Constraints2Polyhedron(Ct, maxRays);
 886   Matrix_Free(M);
 887   Matrix_Free(Ct);
 888   return Peqs;
 889 } /* Polyhedron_Remove_parm_eqs */
 890
 891
 892 /**
 893  * Given a matrix with m parameterized equations, compress the nb_parms
 894  * parameters and n-m variables so that m variables are integer, and transform
 895  * the variable space into a n-m space by eliminating the m variables (using
 896  * the equalities) the variables to be eliminated are chosen automatically by
 897  * the function.
 898  * <b>Deprecated.</b> Try to use Constraints_fullDimensionize instead.
 899  * @param M the constraints
 900  * @param the number of parameters
 901  * @param validityLattice the the integer lattice underlying the integer
 902  * solutions.
 903 */
 904 Matrix * full_dimensionize(Matrix const * M, int nbParms,
 905                            Matrix ** validityLattice) {
 906   Matrix * Eqs, * Ineqs;
 907   Matrix * permutedEqs, * permutedIneqs;
 908   Matrix * Full_Dim;
 909   Matrix * WVL; /* The Whole Validity Lattice (vars+parms) */
 910   unsigned int i,j;
 911   int nbElimVars;
 912   unsigned int * permutation, * permutationInv;
 913   /* 0- Split the equalities and inequalities from each other */
 914   split_constraints(M, &Eqs, &Ineqs);
 915
 916   /* 1- if the polyhedron is already full-dimensional, return it */
 917   if (Eqs->NbRows==0) {
 918     Matrix_Free(Eqs);
 919     (*validityLattice) = Identity_Matrix(nbParms+1);
 920     return Ineqs;
 921   }
 922   nbElimVars = Eqs->NbRows;
 923
 924   /* 2- put the vars to be eliminated at the first positions,
 925      and compress the other vars/parms
 926      -> [ variables to eliminate / parameters / variables to keep ] */
 927   permutation = find_a_permutation(Eqs, nbParms);
 928   if (dbgCompParm) {
 929     printf("Permuting the vars/parms this way: [ ");
 930     for (i=0; i< Eqs->NbColumns; i++) {
 931       printf("%d ", permutation[i]);
 932     }
 933     printf("]\n");
 934   }
 935   permutedEqs = mpolyhedron_permute(Eqs, permutation);
 936   WVL = compress_parms(permutedEqs, Eqs->NbColumns-2-Eqs->NbRows);
 937   /* Check for empty WVL (no solution) */
 938   if( !WVL )
 939   {
 940           fprintf(stderr,"full_dimensionize > parameters compression failed.\n");
 941           Matrix_Free(Eqs);
 942           Matrix_Free(Ineqs);
 943           (*validityLattice) = Identity_Matrix(nbParms+1);
 944           return NULL;
 945   }
 946   if (dbgCompParm) {
 947     printf("Whole validity lattice: ");
 948     show_matrix(WVL);
 949   }
 950   mpolyhedron_compress_last_vars(permutedEqs, WVL);
 951   permutedIneqs = mpolyhedron_permute(Ineqs, permutation);
 952   if (dbgCompParm) {
 953     show_matrix(permutedEqs);
 954   }
 955   Matrix_Free(Eqs);
 956   Matrix_Free(Ineqs);
 957   mpolyhedron_compress_last_vars(permutedIneqs, WVL);
 958   if (dbgCompParm) {
 959     printf("After compression: ");
 960     show_matrix(permutedIneqs);
 961   }
 962   /* 3- eliminate the first variables */
 963   if (!mpolyhedron_eliminate_first_variables(permutedEqs, permutedIneqs)) {
 964     fprintf(stderr,"full_dimensionize > variable elimination failed.\n");
 965          Matrix_Free(permutedIneqs);
 966          (*validityLattice) = Identity_Matrix(nbParms+1);
 967     return NULL;
 968   }
 969   if (dbgCompParm) {
 970     printf("After elimination of the variables: ");
 971     show_matrix(permutedIneqs);
 972   }
 973
 974   /* 4- get rid of the first (zero) columns,
 975      which are now useless, and put the parameters back at the end */
 976   Full_Dim = Matrix_Alloc(permutedIneqs->NbRows,
 977                           permutedIneqs->NbColumns-nbElimVars);
 978   for (i=0; i< permutedIneqs->NbRows; i++) {
 979     value_set_si(Full_Dim->p[i][0], 1);
 980     for (j=0; j< nbParms; j++)
 981       value_assign(Full_Dim->p[i][j+Full_Dim->NbColumns-nbParms-1],
 982                    permutedIneqs->p[i][j+nbElimVars+1]);
 983     for (j=0; j< permutedIneqs->NbColumns-nbParms-2-nbElimVars; j++)
 984       value_assign(Full_Dim->p[i][j+1],
 985                    permutedIneqs->p[i][nbElimVars+nbParms+j+1]);
 986     value_assign(Full_Dim->p[i][Full_Dim->NbColumns-1],
 987                  permutedIneqs->p[i][permutedIneqs->NbColumns-1]);
 988   }
 989   Matrix_Free(permutedIneqs);
 990
 991   /* 5- Keep only the the validity lattice restricted to the parameters */
 992   *validityLattice = Matrix_Alloc(nbParms+1, nbParms+1);
 993   for (i=0; i< nbParms; i++) {
 994     for (j=0; j< nbParms; j++)
 995       value_assign((*validityLattice)->p[i][j],
 996                    WVL->p[i][j]);
 997     value_assign((*validityLattice)->p[i][nbParms],
 998                  WVL->p[i][WVL->NbColumns-1]);
 999   }
1000   for (j=0; j< nbParms; j++)
1001     value_set_si((*validityLattice)->p[nbParms][j], 0);
1002   value_assign((*validityLattice)->p[nbParms][nbParms],
1003                WVL->p[WVL->NbColumns-1][WVL->NbColumns-1]);
1004
1005   /* 6- Clean up */
1006   Matrix_Free(WVL);
1007   return Full_Dim;
1008 } /* full_dimensionize */
1009
1010 #undef dbgCompParm
1011 #undef dbgCompParmMore