README: update note on shared libraries and NTL
[barvinok.git] / polysign_cdd_template.cc
blob1d325ac502490e7c9dd04a9ee1bcf312e5314322
1 #include <assert.h>
2 #define GMPRATIONAL
3 #include <setoper.h>
4 #include <cdd.h>
5 #include <barvinok/util.h>
6 #include <barvinok/options.h>
7 #include "polysign.h"
8 #include "initcdd.h"
10 #define EMPTY_DOMAIN -2
12 static DD_LPType *solve_lp(DD_LPObjectiveType obj, Matrix *C,
13 Value *f)
15 unsigned dim = C->NbColumns - 2;
16 DD_LPType *lp;
17 DD_rowrange irev = C->NbRows;
18 DD_rowrange rows, cols;
19 unsigned NbEq = 0;
21 for (int i = 0; i < C->NbRows; ++i)
22 if (value_zero_p(C->p[i][0]))
23 NbEq++;
25 rows = irev + NbEq + 1;
26 cols = 1 + dim;
27 lp = DD_CreateLPData(obj, DD_Rational, rows, cols);
28 lp->Homogeneous = DD_FALSE;
29 lp->objective = obj;
31 for (DD_rowrange j = 0; j < C->NbRows; ++j) {
32 for (DD_colrange k = 0; k < dim; ++k)
33 DD_set_z(lp->A[j][1+k], C->p[j][1+k]);
34 DD_set_z(lp->A[j][0], C->p[j][1+dim]);
35 if (value_zero_p(C->p[j][0])) {
36 set_addelem(lp->equalityset, j+1);
37 for (DD_colrange k = 0; k < dim; ++k)
38 DD_neg(lp->A[irev][1+k], lp->A[j][1+k]);
39 DD_neg(lp->A[irev][0], lp->A[j][0]);
40 ++irev;
43 /* objective function */
44 for (DD_colrange k = 0; k < dim; ++k)
45 DD_set_z(lp->A[rows-1][1+k], f[k]);
46 DD_set_z(lp->A[rows-1][0], f[dim]);
48 DD_ErrorType err = DD_NoError;
49 DD_LPSolve(lp, DD_DualSimplex, &err);
50 assert(err == DD_NoError);
52 return lp;
55 static lp_result constraints_affine_minmax(DD_LPObjectiveType obj, Matrix *C,
56 Value *f, Value denom, Value *opt)
58 lp_result res = lp_ok;
59 DD_LPType *lp = solve_lp(obj, C, f);
60 assert(value_one_p(denom));
62 switch(lp->LPS) {
63 case DD_Optimal:
64 if (obj == DD_LPmin)
65 DD_ceil(*opt, lp->optvalue);
66 else
67 DD_floor(*opt, lp->optvalue);
68 break;
69 case DD_DualInconsistent:
70 res = lp_unbounded;
71 break;
72 case DD_Inconsistent:
73 res = lp_empty;
74 break;
75 default:
76 assert(0);
78 DD_FreeLPData(lp);
79 return res;
82 static int constraints_affine_minmax_sign(DD_LPObjectiveType obj, Matrix *C,
83 Matrix *T, bool rational)
85 unsigned dim = C->NbColumns-2;
86 assert(dim == T->NbColumns-1);
87 assert(T->NbRows == 2);
88 DD_LPType *lp = solve_lp(obj, C, T->p[0]);
90 int sign;
91 if (lp->LPS == DD_Optimal) {
92 if (rational)
93 DD_rat_sign(sign, obj, lp->optvalue);
94 else
95 /* The objective function has integer coefficients,
96 * so the optimal should be an integer (over the integer points)
98 DD_int_sign(sign, obj, lp->optvalue);
99 } else if (lp->LPS == DD_DualInconsistent) {
100 if (obj == DD_LPmin)
101 sign = -1;
102 else
103 sign = 1;
104 } else if (lp->LPS == DD_Inconsistent) {
105 sign = EMPTY_DOMAIN;
106 } else
107 assert(0);
109 DD_FreeLPData(lp);
110 return sign;
113 enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
114 struct barvinok_options *options)
116 Matrix M;
118 if (emptyQ2(D))
119 return order_undefined;
121 INIT_CDD;
122 bool rational = !POL_ISSET(options->MaxRays, POL_INTEGER);
123 Polyhedron_Matrix_View(D, &M, D->NbConstraints);
124 int min = constraints_affine_minmax_sign(DD_LPmin, &M, T, rational);
125 if (min == EMPTY_DOMAIN)
126 return order_undefined;
127 if (min > 0)
128 return order_gt;
129 int max = constraints_affine_minmax_sign(DD_LPmax, &M, T, rational);
130 assert(max != EMPTY_DOMAIN);
131 if (max < 0)
132 return order_lt;
133 if (min == max)
134 return order_eq;
135 if (max == 0)
136 return order_le;
137 if (min == 0)
138 return order_ge;
139 return order_unknown;
142 enum lp_result cdd_constraints_opt(Matrix *C, Value *obj, Value denom,
143 enum lp_dir dir, Value *opt)
145 DD_LPObjectiveType cdd_dir = dir == lp_min ? DD_LPmin : DD_LPmax;
147 INIT_CDD;
148 return constraints_affine_minmax(cdd_dir, C, obj, denom, opt);