remove bernstein
[barvinok.git] / reduce_domain.c
blob5b6fd9bf9cfbe92f2bdbc2637cb404319d62970b
1 #include <assert.h>
2 #include <barvinok/options.h>
3 #include <barvinok/util.h>
4 #include "reduce_domain.h"
6 Polyhedron *true_context(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
8 unsigned nparam = C->Dimension;
9 Polyhedron *tmp = Polyhedron_Project(P, nparam);
10 C = DomainIntersection(tmp, C, MaxRays);
11 Polyhedron_Free(tmp);
12 return C;
15 Vector *inner_point(Polyhedron *P)
17 Vector *average;
18 int nv = 0;
19 Value fc, fv;
20 unsigned dim = P->Dimension;
21 int i;
23 average = Vector_Alloc(1+dim+1);
25 POL_ENSURE_VERTICES(P);
26 for (i = 0; i < P->NbRays; ++i) {
27 if (value_zero_p(P->Ray[i][1+dim]))
28 continue;
29 ++nv;
30 if (nv == 1) {
31 Vector_Copy(P->Ray[i]+1, average->p+1, dim+1);
32 continue;
34 if (nv == 2) {
35 value_init(fc);
36 value_init(fv);
38 value_assign(fc, average->p[1+dim]);
39 value_lcm(average->p[1+dim], fc, P->Ray[i][1+dim]);
40 value_division(fc, average->p[1+dim], fc);
41 value_division(fv, average->p[1+dim], P->Ray[i][1+dim]);
42 Vector_Combine(average->p+1, P->Ray[i]+1, average->p+1, fc, fv, dim);
44 if (nv > 1) {
45 value_set_si(fc, nv);
46 value_multiply(average->p[1+dim], average->p[1+dim], fc);
47 Vector_Normalize(average->p+1, dim+1);
48 value_clear(fc);
49 value_clear(fv);
51 for (i = 0; i < P->NbRays; ++i) {
52 if (value_notzero_p(P->Ray[i][1+dim]))
53 continue;
54 Vector_Add(average->p+1, P->Ray[i]+1, average->p+1, dim);
56 return average;
59 int is_internal(Vector *point, Value *constraint)
61 int p;
62 unsigned dim = point->Size-2;
64 Inner_Product(constraint+1, point->p+1, dim+1, &point->p[0]);
65 if (value_notzero_p(point->p[0]))
66 return value_pos_p(point->p[0]);
68 p = First_Non_Zero(constraint+1, dim);
69 return value_pos_p(constraint[1+p]);
72 Polyhedron *reduce_domain(Polyhedron *D, int nd,
73 Vector *inner, struct barvinok_options *options)
75 Polyhedron *rVD;
76 Value c;
77 int i;
78 Matrix *constraints;
79 int changed;
81 rVD = D->next ? DomainConvex(D, options->MaxRays) : Polyhedron_Copy(D);
83 /* If there is only one chamber, then we don't need to take care
84 * of possible overlaps.
85 * Plus, this decomposition may be the result of a recursive call
86 * and then some of the assumptions used in determining whether
87 * the domain is too small in geometric dimension no longer apply.
89 if (nd == 1) {
90 return rVD;
93 assert(rVD->Dimension == inner->Size-2);
94 constraints = Polyhedron2Constraints(rVD);
95 changed = 0;
96 for (i = 0; i < constraints->NbRows; ++i) {
97 if (!is_internal(inner, constraints->p[i])) {
98 value_decrement(constraints->p[i][1+rVD->Dimension],
99 constraints->p[i][1+rVD->Dimension]);
100 changed = 1;
103 if (changed) {
104 Polyhedron_Free(rVD);
105 rVD = Constraints2Polyhedron(constraints, options->MaxRays);
107 Matrix_Free(constraints);
109 rVD = DomainConstraintSimplify(rVD, options->MaxRays);
110 POL_ENSURE_FACETS(rVD);
111 if (emptyQ(rVD)) {
112 Domain_Free(rVD);
113 return 0;
116 value_init(c);
117 barvinok_count_with_options(rVD, &c, options);
118 if (value_zero_p(c)) {
119 Domain_Free(rVD);
120 rVD = 0;
122 value_clear(c);
124 return rVD;