evalue.c: evalue_frac2polynomial: improve accuracy
[barvinok.git] / DomainConstraintSimplify.c
blob64e00246cdcc33c0edfd86e709ed05112b9e6fdb
1 #include <barvinok/polylib.h>
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.
11 int ConstraintSimplify(Value *old, Value *new, int len, Value* v)
13 Vector_Gcd(old+1, len - 2, v);
15 if (value_one_p(*v))
16 return 0;
18 Vector_AntiScale(old+1, new+1, *v, len-2);
19 mpz_fdiv_q(new[len-1], old[len-1], *v);
20 return 1;
23 static Polyhedron *p_simplify_constraints(Polyhedron *P, Vector *row,
24 Value *g, unsigned MaxRays)
26 Polyhedron *T, *R = P;
27 int len = P->Dimension+2;
28 int r;
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
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;
44 if (R != P)
45 Polyhedron_Free(P);
46 return R;
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
55 Polyhedron *DomainConstraintSimplify(Polyhedron *P, unsigned MaxRays)
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);
65 for (prev = &R; P; P = N) {
66 Polyhedron *T;
67 N = P->next;
68 T = p_simplify_constraints(P, row, &g, MaxRays);
70 if (emptyQ(T) && prev != &R) {
71 Polyhedron_Free(T);
72 *prev = NULL;
73 continue;
76 if (T != P)
77 T->next = N;
78 *prev = T;
79 prev = &T->next;
82 if (R->next && emptyQ(R)) {
83 N = R->next;
84 Polyhedron_Free(R);
85 R = N;
88 value_clear(g);
89 Vector_Free(row);
90 return R;