polysign_cdd_template.cc

   1 #define GMPRATIONAL
   2 #include <setoper.h>
   3 #include <cdd.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/options.h>
   6 #include "polysign.h"
   7 #include "initcdd.h"
   8
   9 #define EMPTY_DOMAIN    -2
  10
  11 static DD_LPType *solve_lp(DD_LPObjectiveType obj, Polyhedron *P,
  12                            Value *f)
  13 {
  14     DD_LPType *lp;
  15     DD_rowrange irev = P->NbConstraints;
  16     DD_rowrange rows = irev + P->NbEq + 1;
  17     DD_colrange cols = 1 + P->Dimension;
  18     lp = DD_CreateLPData(obj, DD_Rational, rows, cols);
  19     lp->Homogeneous = DD_FALSE;
  20     lp->objective = obj;
  21
  22     for (DD_rowrange j = 0; j < P->NbConstraints; ++j) {
  23         for (DD_colrange k = 0; k < P->Dimension; ++k) {
  24             DD_set_si(lp->A[j][1+k], 5);
  25             DD_set_z(lp->A[j][1+k], P->Constraint[j][1+k]);
  26         }
  27         DD_set_z(lp->A[j][0], P->Constraint[j][1+P->Dimension]);
  28         if (j < P->NbEq) {
  29             set_addelem(lp->equalityset, j+1);
  30             for (DD_colrange k = 0; k < P->Dimension; ++k)
  31                 DD_neg(lp->A[irev][1+k], lp->A[j][1+k]);
  32             DD_neg(lp->A[irev][0], lp->A[j][0]);
  33             ++irev;
  34         }
  35     }
  36     /* objective function */
  37     for (DD_colrange k = 0; k < P->Dimension; ++k)
  38         DD_set_z(lp->A[rows-1][1+k], f[k]);
  39     DD_set_z(lp->A[rows-1][0], f[P->Dimension]);
  40
  41     DD_ErrorType err = DD_NoError;
  42     DD_LPSolve(lp, DD_DualSimplex, &err);
  43     assert(err == DD_NoError);
  44
  45     return lp;
  46 }
  47
  48 static lp_result polyhedron_affine_minmax(DD_LPObjectiveType obj, Polyhedron *P,
  49                                           Value *f, Value denom, Value *opt)
  50 {
  51     lp_result res = lp_ok;
  52     DD_LPType *lp = solve_lp(obj, P, f);
  53     assert(value_one_p(denom));
  54
  55     switch(lp->LPS) {
  56     case DD_Optimal:
  57         if (obj == DD_LPmin)
  58             DD_ceil(*opt, lp->optvalue);
  59         else
  60             DD_floor(*opt, lp->optvalue);
  61         break;
  62     case DD_DualInconsistent:
  63         res = lp_unbounded;
  64         break;
  65     case DD_Inconsistent:
  66         res = lp_empty;
  67         break;
  68     default:
  69         assert(0);
  70     }
  71     DD_FreeLPData(lp);
  72     return res;
  73 }
  74
  75 static int polyhedron_affine_minmax_sign(DD_LPObjectiveType obj, Polyhedron *P,
  76                                          Matrix *T, bool rational)
  77 {
  78     assert(P->Dimension == T->NbColumns-1);
  79     assert(T->NbRows == 2);
  80     DD_LPType *lp = solve_lp(obj, P, T->p[0]);
  81
  82     int sign;
  83     if (lp->LPS == DD_Optimal) {
  84         if (rational)
  85             DD_rat_sign(sign, obj, lp->optvalue);
  86         else
  87             /* The objective function has integer coefficients,
  88              * so the optimal should be an integer (over the integer points)
  89              */
  90             DD_int_sign(sign, obj, lp->optvalue);
  91     } else if (lp->LPS == DD_DualInconsistent) {
  92         if (obj == DD_LPmin)
  93             sign = -1;
  94         else
  95             sign = 1;
  96     } else if (lp->LPS == DD_Inconsistent) {
  97         sign = EMPTY_DOMAIN;
  98     } else
  99         assert(0);
 100
 101     DD_FreeLPData(lp);
 102     return sign;
 103 }
 104
 105 enum order_sign cdd_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
 106                                             struct barvinok_options *options)
 107 {
 108     if (emptyQ2(D))
 109         return order_undefined;
 110
 111     INIT_CDD;
 112     bool rational = !POL_ISSET(options->MaxRays, POL_INTEGER);
 113     int min = polyhedron_affine_minmax_sign(DD_LPmin, D, T, rational);
 114     if (min == EMPTY_DOMAIN)
 115         return order_undefined;
 116     if (min > 0)
 117         return order_gt;
 118     int max = polyhedron_affine_minmax_sign(DD_LPmax, D, T, rational);
 119     assert(max != EMPTY_DOMAIN);
 120     if (max < 0)
 121         return order_lt;
 122     if (min == max)
 123         return order_eq;
 124     if (max == 0)
 125         return order_le;
 126     if (min == 0)
 127         return order_ge;
 128     return order_unknown;
 129 }
 130
 131 enum lp_result cdd_polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 132                                 Value *min, Value *max,
 133                                 struct barvinok_options *options)
 134 {
 135     lp_result res;
 136
 137     if (emptyQ2(D))
 138         return lp_empty;
 139
 140     INIT_CDD;
 141     res = polyhedron_affine_minmax(DD_LPmin, D, obj, denom, min);
 142     if (res != lp_ok)
 143         return res;
 144     res = polyhedron_affine_minmax(DD_LPmax, D, obj, denom, max);
 145     return res;
 146 }