add tests/ehrhart.README explaining origin of some of the ehrhart inputs
[barvinok.git] / polysign_polylib.c
blob9287d00195b6b26c0050fbb7e9c16b8c35fa3c76
1 #include <assert.h>
2 #include <barvinok/util.h>
3 #include <barvinok/options.h>
4 #include "polysign.h"
6 static enum order_sign interval_minmax(Polyhedron *I)
8 int i;
9 int min = 1;
10 int max = -1;
11 assert(I->Dimension == 1);
12 POL_ENSURE_VERTICES(I);
13 for (i = 0; i < I->NbRays; ++i) {
14 if (value_pos_p(I->Ray[i][1]))
15 max = 1;
16 else if (value_neg_p(I->Ray[i][1]))
17 min = -1;
18 else {
19 if (max < 0)
20 max = 0;
21 if (min > 0)
22 min = 0;
25 if (min > 0)
26 return order_gt;
27 if (max < 0)
28 return order_lt;
29 if (min == max)
30 return order_eq;
31 if (max == 0)
32 return order_le;
33 if (min == 0)
34 return order_ge;
35 return order_unknown;
38 /* Returns the sign of the affine function specified by T on the polyhedron D */
39 enum order_sign PL_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
40 struct barvinok_options *options)
42 enum order_sign sign;
43 Polyhedron *I = Polyhedron_Image(D, T, options->MaxRays);
44 if (POL_ISSET(options->MaxRays, POL_INTEGER))
45 I = DomainConstraintSimplify(I, options->MaxRays);
46 if (emptyQ2(I)) {
47 Polyhedron_Free(I);
48 I = Polyhedron_Image(D, T, options->MaxRays);
50 sign = interval_minmax(I);
51 Polyhedron_Free(I);
52 return sign;
55 enum lp_result PL_constraints_opt(Matrix *C, Value *obj, Value denom,
56 enum lp_dir dir, Value *opt,
57 unsigned MaxRays)
59 enum lp_result res;
60 Polyhedron *P = Constraints2Polyhedron(C, MaxRays);
61 res = PL_polyhedron_opt(P, obj, denom, dir, opt);
62 Polyhedron_Free(P);
63 return res;
66 enum lp_result PL_polyhedron_opt(Polyhedron *P, Value *obj, Value denom,
67 enum lp_dir dir, Value *opt)
69 int i;
70 int first = 1;
71 Value val, d;
72 enum lp_result res = lp_empty;
74 POL_ENSURE_VERTICES(P);
75 if (emptyQ(P))
76 return res;
78 value_init(val);
79 value_init(d);
80 for (i = 0; i < P->NbRays; ++ i) {
81 Inner_Product(P->Ray[i]+1, obj, P->Dimension+1, &val);
82 if (value_zero_p(P->Ray[i][0]) && value_notzero_p(val)) {
83 res = lp_unbounded;
84 break;
86 if (value_zero_p(P->Ray[i][1+P->Dimension])) {
87 if ((dir == lp_min && value_neg_p(val)) ||
88 (dir == lp_max && value_pos_p(val))) {
89 res = lp_unbounded;
90 break;
92 } else {
93 res = lp_ok;
94 value_multiply(d, denom, P->Ray[i][1+P->Dimension]);
95 if (dir == lp_min)
96 mpz_cdiv_q(val, val, d);
97 else
98 mpz_fdiv_q(val, val, d);
99 if (first || (dir == lp_min ? value_lt(val, *opt) :
100 value_gt(val, *opt)))
101 value_assign(*opt, val);
102 first = 0;
105 value_clear(d);
106 value_clear(val);
108 return res;