source/kernel/compress_parms.c

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