DomainConstraintSimplify.c

   1 #include <polylib/polylibgmp.h>
   2
   3 /*
   4  * Replaces constraint a x >= c by x >= ceil(c/a)
   5  * where "a" is a common factor in the coefficients
   6  * old is the constraint; v points to an initialized
   7  * value that this procedure can use.
   8  * Return non-zero if something changed.
   9  * Result is placed in new.
  10  */
  11 int ConstraintSimplify(Value *old, Value *new, int len, Value* v)
  12 {
  13     Vector_Gcd(old+1, len - 2, v);
  14
  15     if (value_one_p(*v))
  16         return 0;
  17
  18     Vector_AntiScale(old+1, new+1, *v, len-2);
  19     mpz_fdiv_q(new[len-1], old[len-1], *v);
  20     return 1;
  21 }
  22
  23 static Polyhedron *p_simplify_constraints(Polyhedron *P, Vector *row,
  24                                           Value *g, unsigned MaxRays)
  25 {
  26     Polyhedron *T, *R = P;
  27     int len = P->Dimension+2;
  28     int r;
  29
  30     /* Also look at equalities.
  31      * If an equality can be "simplified" then there
  32      * are no integer solutions anyway and the following loop
  33      * will add a conflicting constraint
  34      */
  35     for (r = 0; r < R->NbConstraints; ++r) {
  36         if (ConstraintSimplify(R->Constraint[r], row->p, len, g)) {
  37             T = R;
  38             R = AddConstraints(row->p, 1, R, MaxRays);
  39             if (T != P)
  40                 Polyhedron_Free(T);
  41             r = -1;
  42         }
  43     }
  44     if (R != P)
  45         Polyhedron_Free(P);
  46     return R;
  47 }
  48
  49 /*
  50  * Replaces constraint a x >= c by x >= ceil(c/a)
  51  * where "a" is a common factor in the coefficients
  52  * Destroys P and returns a newly allocated Polyhedron
  53  * or just returns P in case no changes were made
  54  */
  55 Polyhedron *DomainConstraintSimplify(Polyhedron *P, unsigned MaxRays)
  56 {
  57     Polyhedron **prev;
  58     int len = P->Dimension+2;
  59     Vector *row = Vector_Alloc(len);
  60     Value g;
  61     Polyhedron *R = P, *N;
  62     value_set_si(row->p[0], 1);
  63     value_init(g);
  64
  65     for (prev = &R; P; P = N) {
  66         Polyhedron *T;
  67         N = P->next;
  68         T = p_simplify_constraints(P, row, &g, MaxRays);
  69
  70         if (emptyQ(T) && prev != &R) {
  71             Polyhedron_Free(T);
  72             *prev = NULL;
  73             continue;
  74         }
  75
  76         if (T != P)
  77             T->next = N;
  78         *prev = T;
  79         prev = &T->next;
  80     }
  81
  82     if (R->next && emptyQ(R)) {
  83         N = R->next;
  84         Polyhedron_Free(R);
  85         R = N;
  86     }
  87
  88     value_clear(g);
  89     Vector_Free(row);
  90     return R;
  91 }