normalization.c

   1 #include <assert.h>
   2 #include "normalization.h"
   3
   4 static int is_unit_row(Value *row, int pos, int len)
   5 {
   6     if (!value_one_p(row[pos]) && !value_mone_p(row[pos]))
   7         return 0;
   8     return First_Non_Zero(row+pos+1, len-(pos+1)) == -1;
   9 }
  10
  11 /* Transform the constraints C to "standard form".
  12  * In particular, if C is described by
  13  *              A x + b(p) >= 0
  14  * then this function returns a matrix H = A U, A = H Q, such
  15  * that D x' = D Q x >= -b(p), with D a diagonal matrix with
  16  * positive entries.  The calling function can then construct
  17  * the standard form H' x' - I s + b(p) = 0, with H' the rows of H
  18  * that are not positive multiples of unit vectors
  19  * (since those correspond to D x' >= -b(p)).
  20  * The number of rows in H' is returned in *rows_p.
  21  * Optionally returns the matrix that maps the new variables
  22  * back to the old variables x = U x'.
  23  * Note that the rows of H (and C) may be reordered by this function.
  24  */
  25 Matrix *standard_constraints(Matrix *C, unsigned nparam, int *rows_p,
  26                              Matrix **T)
  27 {
  28     unsigned dim = C->NbColumns - 2;
  29     unsigned nvar = dim - nparam;
  30     int i, j, d;
  31     int rows;
  32     Matrix *M;
  33     Matrix *H, *U, *Q;
  34
  35     /* move constraints only involving parameters down
  36      * and move unit vectors (if there are any) to the right place.
  37      */
  38     for (d = 0, j = C->NbRows; d < j; ++d) {
  39         int index;
  40         index = First_Non_Zero(C->p[d]+1, nvar);
  41         if (index != -1) {
  42             if (index != d &&
  43                 is_unit_row(C->p[d]+1, index, nvar)) {
  44                 Vector_Exchange(C->p[d], C->p[index], dim+2);
  45                 if (index > d)
  46                     --d;
  47             }
  48             continue;
  49         }
  50         while (d < --j && First_Non_Zero(C->p[j]+1, nvar) == -1)
  51             ;
  52         if (d >= j)
  53             break;
  54         Vector_Exchange(C->p[d], C->p[j], dim+2);
  55     }
  56     M = Matrix_Alloc(d, nvar);
  57     for (j = 0; j < d; ++j)
  58         Vector_Copy(C->p[j]+1, M->p[j], nvar);
  59
  60     neg_left_hermite(M, &H, &Q, &U);
  61     Matrix_Free(M);
  62     Matrix_Free(Q);
  63     if (T)
  64         *T = U;
  65     else
  66         Matrix_Free(U);
  67
  68     /* Rearrange rows such that top of H is lower diagonal and
  69      * compute the number of non (multiple of) unit-vector rows.
  70      */
  71     rows = H->NbRows-nvar;
  72     for (i = 0; i < H->NbColumns; ++i) {
  73         for (j = i; j < H->NbRows; ++j)
  74             if (value_notzero_p(H->p[j][i]))
  75                 break;
  76         if (j != i) {
  77             Vector_Exchange(C->p[i], C->p[j], dim+2);
  78             Vector_Exchange(H->p[i], H->p[j], H->NbColumns);
  79         }
  80         if (First_Non_Zero(H->p[i], i) != -1)
  81             rows++;
  82     }
  83     if (rows_p)
  84         *rows_p = rows;
  85
  86     return H;
  87 }
  88
  89 /*
  90  * Transform the constraints of P, which may have existentially
  91  * quantified variables to "standard form".
  92  *
  93  * We need to make sure that the new variables are affine
  94  * combinations of only the original variables (and not
  95  * of any of the existentially quantified variables).
  96  * We therefore first apply standard_constraints on
  97  * the constraints involving only the variables and
  98  * then apply the same procedure on the remaining
  99  * constraints with the variables as parameters.
 100  *
 101  * Perhaps we are overly careful here.  Moving the constraints
 102  * involving existentially quantified variables down should prevent
 103  * the regular variables from being replaced by linear combinations
 104  * of existential variables.
 105  */
 106 static Matrix *exist_standard_constraints(Polyhedron *P, unsigned nvar)
 107 {
 108     int i, i1, i2;
 109     unsigned exist = P->Dimension - nvar;
 110     Matrix *M = Matrix_Alloc(P->NbConstraints, 2+P->Dimension);
 111     Matrix *M1 = Matrix_Alloc(P->NbConstraints, 2+nvar);
 112     Matrix *M2 = Matrix_Alloc(P->NbConstraints, 2+P->Dimension);
 113     Matrix *S1, *S2;
 114     Value c;
 115
 116     for (i = 0; i < P->NbConstraints; ++i) {
 117         value_assign(M->p[i][0], P->Constraint[i][0]);
 118         Vector_Copy(P->Constraint[i]+1, M->p[i]+1+exist, nvar);
 119         Vector_Copy(P->Constraint[i]+1+nvar, M->p[i]+1, exist);
 120         value_assign(M->p[i][1+P->Dimension], P->Constraint[i][1+P->Dimension]);
 121     }
 122     Gauss(M, P->NbEq, P->Dimension+1);
 123     for (i = i1 = i2 = 0; i < M->NbRows; ++i) {
 124         if (First_Non_Zero(M->p[i]+1, exist) == -1) {
 125             value_assign(M1->p[i1][0], M->p[i][0]);
 126             Vector_Copy(M->p[i]+1+exist, M1->p[i1]+1, nvar+1);
 127             ++i1;
 128         } else {
 129             Vector_Copy(M->p[i], M2->p[i2], M->NbColumns);
 130             ++i2;
 131         }
 132     }
 133     M1->NbRows = i1;
 134     M2->NbRows = i2;
 135     Matrix_Free(M);
 136
 137     S1 = standard_constraints(M1, 0, NULL, NULL);
 138     S2 = standard_constraints(M2, nvar, NULL, NULL);
 139
 140     M = Matrix_Alloc(P->NbConstraints, 2+P->Dimension);
 141     for (i = 0; i < M1->NbRows; ++i) {
 142         value_assign(M->p[i][0], M1->p[i][0]);
 143         Vector_Copy(S1->p[i], M->p[i]+1, nvar);
 144         value_assign(M->p[i][1+P->Dimension], M1->p[i][1+nvar]);
 145     }
 146     for (i = 0; i < M2->NbRows; ++i) {
 147         value_assign(M->p[M1->NbRows+i][0], M2->p[i][0]);
 148         Vector_Copy(M2->p[i]+1+exist, M->p[M1->NbRows+i]+1, nvar);
 149         Vector_Copy(S2->p[i], M->p[M1->NbRows+i]+1+nvar, exist);
 150         value_assign(M->p[M1->NbRows+i][1+P->Dimension],
 151                         M2->p[i][1+P->Dimension]);
 152     }
 153
 154     /* Make sure "standard" constraints for existential variables
 155      * have only non-positive coefficients for the variables.
 156      */
 157     assert(M2->NbRows >= exist);
 158     value_init(c);
 159     for (i = 0; i < exist; ++i) {
 160         assert(value_pos_p(S2->p[i][i]));
 161         for (i1 = 0; i1 < nvar; ++i1) {
 162             if (value_negz_p(M->p[M1->NbRows+i][1+i1]))
 163                 continue;
 164             mpz_cdiv_q(c, M->p[M1->NbRows+i][1+i1], S2->p[i][i]);
 165             value_oppose(c, c);
 166             for (i2 = i; i2 < M2->NbRows; ++i2)
 167                 value_addmul(M->p[M1->NbRows+i2][1+i1], c, S2->p[i2][i]);
 168         }
 169     }
 170     value_clear(c);
 171
 172     Matrix_Free(M1);
 173     Matrix_Free(M2);
 174
 175     return M;
 176 }
 177
 178 /*
 179  * For those variables x that may attain negative values,
 180  * compute a translations t such that x' = x + t is sure
 181  * to attain only non-negative values, i.e., if M is
 182  *
 183  *              A x + b >= 0
 184  *
 185  * then the x' in
 186  *
 187  *              A x' + b - A t >= 0
 188  *
 189  * are sure to attain only non-negative values.
 190  */
 191 static Vector *compute_shifts(Matrix *M)
 192 {
 193     int c, r;
 194     unsigned dim = M->NbColumns-2;
 195     Vector *shifts = Vector_Alloc(M->NbColumns);
 196     for (c = 0, r = 0; c < dim; ++c) {
 197         for ( ; r < M->NbRows; ++r) {
 198             if (value_notzero_p(M->p[r][1+c]))
 199                 break;
 200         }
 201         assert(r < M->NbRows);
 202         assert(value_pos_p(M->p[r][1+c]));
 203         if (c > 0)
 204             Inner_Product(M->p[r]+1, shifts->p+1, c, &shifts->p[1+c]);
 205         value_subtract(shifts->p[1+c], M->p[r][1+dim], shifts->p[1+c]);
 206         if (value_pos_p(shifts->p[1+c]))
 207             mpz_cdiv_q(shifts->p[1+c], shifts->p[1+c], M->p[r][1+c]);
 208         else
 209             value_set_si(shifts->p[1+c], 0);
 210     }
 211     return shifts;
 212 }
 213
 214 /*
 215  * Given a union of polyhedra, with possibly existentially
 216  * quantified variables, but no parameters, apply a unimodular
 217  * transformation with constant translation to the variables,
 218  * and a unimodular transformation with a translation possibly
 219  * depending on the variables to the existentially quantified variables,
 220  * such that all variables and all existentially quantified variables
 221  * attain only non-negative values.
 222  */
 223 Polyhedron *skew_to_positive_orthant(Polyhedron *D, unsigned nvar,
 224                                         unsigned MaxRays)
 225 {
 226     Polyhedron *P;
 227     Polyhedron *R = NULL;
 228     Polyhedron **next = &R;
 229
 230     for (P = D; P; P = P->next) {
 231         Matrix *M, *C;
 232         Vector *shifts, *offsets;
 233         int i;
 234
 235         M = exist_standard_constraints(P, nvar);
 236         shifts = compute_shifts(M);
 237
 238         offsets = Vector_Alloc(P->NbConstraints);
 239         Matrix_Vector_Product(M, shifts->p, offsets->p);
 240         Vector_Free(shifts);
 241
 242         C = Matrix_Alloc(P->NbConstraints, P->Dimension+2);
 243         for (i = 0; i < P->NbConstraints; ++i) {
 244             value_assign(C->p[i][0], M->p[i][0]);
 245             Vector_Copy(M->p[i]+1, C->p[i]+1, P->Dimension);
 246             value_subtract(C->p[i][1+P->Dimension],
 247                             M->p[i][1+P->Dimension], offsets->p[i]);
 248         }
 249         Vector_Free(offsets);
 250         Matrix_Free(M);
 251
 252         *next = Constraints2Polyhedron(C, MaxRays);
 253         Matrix_Free(C);
 254
 255         next = &(*next)->next;
 256     }
 257     return R;
 258 }