decomposer.cc

   1 #include <vector>
   2 #include <assert.h>
   3 #include <gmp.h>
   4 #include <NTL/mat_ZZ.h>
   5 #include <NTL/LLL.h>
   6 #include <barvinok/barvinok.h>
   7 #include <barvinok/util.h>
   8 #include "conversion.h"
   9 #include "decomposer.h"
  10
  11 #ifdef NTL_STD_CXX
  12 using namespace NTL;
  13 #endif
  14 using std::vector;
  15
  16 /*
  17  * Returns the largest absolute value in the vector
  18  */
  19 static ZZ max(vec_ZZ& v)
  20 {
  21     ZZ max = abs(v[0]);
  22     for (int i = 1; i < v.length(); ++i)
  23         if (abs(v[i]) > max)
  24             max = abs(v[i]);
  25     return max;
  26 }
  27
  28 class cone {
  29 public:
  30     cone(Matrix *M) {
  31         Cone = 0;
  32         Rays = Matrix_Copy(M);
  33         set_det();
  34     }
  35     cone(Polyhedron *C) {
  36         Cone = Polyhedron_Copy(C);
  37         Rays = rays(C);
  38         set_det();
  39     }
  40     void set_det() {
  41         mat_ZZ A;
  42         matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
  43         det = determinant(A);
  44     }
  45
  46     Vector* short_vector(vec_ZZ& lambda, barvinok_options *options) {
  47         Matrix *M = Matrix_Copy(Rays);
  48         Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
  49         int ok = Matrix_Inverse(M, inv);
  50         assert(ok);
  51         Matrix_Free(M);
  52
  53         ZZ det2;
  54         mat_ZZ B;
  55         mat_ZZ U;
  56         matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
  57         long r = LLL(det2, B, U, options->LLL_a, options->LLL_b);
  58
  59         ZZ min = max(B[0]);
  60         int index = 0;
  61         for (int i = 1; i < B.NumRows(); ++i) {
  62             ZZ tmp = max(B[i]);
  63             if (tmp < min) {
  64                 min = tmp;
  65                 index = i;
  66             }
  67         }
  68
  69         Matrix_Free(inv);
  70
  71         lambda = B[index];
  72
  73         Vector *z = Vector_Alloc(U[index].length()+1);
  74         assert(z);
  75         zz2values(U[index], z->p);
  76         value_set_si(z->p[U[index].length()], 0);
  77
  78         Polyhedron *C = poly();
  79         int i;
  80         for (i = 0; i < lambda.length(); ++i)
  81             if (lambda[i] > 0)
  82                 break;
  83         if (i == lambda.length()) {
  84             Value tmp;
  85             value_init(tmp);
  86             value_set_si(tmp, -1);
  87             Vector_Scale(z->p, z->p, tmp, z->Size-1);
  88             value_clear(tmp);
  89         }
  90         return z;
  91     }
  92
  93     ~cone() {
  94         Polyhedron_Free(Cone);
  95         Matrix_Free(Rays);
  96     }
  97
  98     Polyhedron *poly() {
  99         if (!Cone) {
 100             Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
 101             for (int i = 0; i < Rays->NbRows; ++i) {
 102                 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
 103                 value_set_si(M->p[i][0], 1);
 104             }
 105             Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
 106             value_set_si(M->p[Rays->NbRows][0], 1);
 107             value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
 108             Cone = Rays2Polyhedron(M, M->NbRows+1);
 109             assert(Cone->NbConstraints == Cone->NbRays);
 110             Matrix_Free(M);
 111         }
 112         return Cone;
 113     }
 114
 115     ZZ det;
 116     Polyhedron *Cone;
 117     Matrix *Rays;
 118 };
 119
 120 void decomposer::decompose(Polyhedron *C, barvinok_options *options)
 121 {
 122     vector<cone *> nonuni;
 123     cone * c = new cone(C);
 124     ZZ det = c->det;
 125     int s = sign(det);
 126     assert(det != 0);
 127     if (abs(det) > 1) {
 128         nonuni.push_back(c);
 129     } else {
 130         try {
 131             options->stats.unimodular_cones++;
 132             handle(C, 1);
 133             delete c;
 134         } catch (...) {
 135             delete c;
 136             throw;
 137         }
 138     }
 139     vec_ZZ lambda;
 140     while (!nonuni.empty()) {
 141         c = nonuni.back();
 142         nonuni.pop_back();
 143         Vector* v = c->short_vector(lambda, options);
 144         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 145             if (lambda[i] == 0)
 146                 continue;
 147             Matrix* M = Matrix_Copy(c->Rays);
 148             Vector_Copy(v->p, M->p[i], v->Size);
 149             cone * pc = new cone(M);
 150             assert (pc->det != 0);
 151             if (abs(pc->det) > 1) {
 152                 assert(abs(pc->det) < abs(c->det));
 153                 nonuni.push_back(pc);
 154             } else {
 155                 try {
 156                     options->stats.unimodular_cones++;
 157                     handle(pc->poly(), sign(pc->det) * s);
 158                     delete pc;
 159                 } catch (...) {
 160                     delete c;
 161                     delete pc;
 162                     while (!nonuni.empty()) {
 163                         c = nonuni.back();
 164                         nonuni.pop_back();
 165                         delete c;
 166                     }
 167                     Matrix_Free(M);
 168                     Vector_Free(v);
 169                     throw;
 170                 }
 171             }
 172             Matrix_Free(M);
 173         }
 174         Vector_Free(v);
 175         delete c;
 176     }
 177 }
 178
 179 void polar_decomposer::decompose(Polyhedron *cone, barvinok_options *options)
 180 {
 181     POL_ENSURE_VERTICES(cone);
 182     Polyhedron_Polarize(cone);
 183     if (cone->NbRays - 1 != cone->Dimension) {
 184         Polyhedron *tmp = cone;
 185         cone = triangulate_cone(cone, options->MaxRays);
 186         Polyhedron_Free(tmp);
 187     }
 188     try {
 189         for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
 190             decomposer::decompose(Polar, options);
 191         Domain_Free(cone);
 192     } catch (...) {
 193         Domain_Free(cone);
 194         throw;
 195     }
 196 }
 197
 198 void polar_decomposer::handle(Polyhedron *P, int sign)
 199 {
 200     Polyhedron_Polarize(P);
 201     handle_polar(P, sign);
 202 }
 203
 204 void vertex_decomposer::decompose_at_vertex(Param_Vertices *V, int _i,
 205                                             barvinok_options *options)
 206 {
 207     Polyhedron *C = supporting_cone_p(P, V);
 208     vert = _i;
 209     this->V = V;
 210
 211     pd->decompose(C, options);
 212 }
 213
 214 /*
 215  * Barvinok's Decomposition of a simplicial cone
 216  *
 217  * Returns two lists of polyhedra
 218  */
 219 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
 220 {
 221     barvinok_options *options = barvinok_options_new_with_defaults();
 222     Polyhedron *pos = *ppos, *neg = *pneg;
 223     vector<cone *> nonuni;
 224     cone * c = new cone(C);
 225     ZZ det = c->det;
 226     int s = sign(det);
 227     assert(det != 0);
 228     if (abs(det) > 1) {
 229         nonuni.push_back(c);
 230     } else {
 231         Polyhedron *p = Polyhedron_Copy(c->Cone);
 232         p->next = pos;
 233         pos = p;
 234         delete c;
 235     }
 236     vec_ZZ lambda;
 237     while (!nonuni.empty()) {
 238         c = nonuni.back();
 239         nonuni.pop_back();
 240         Vector* v = c->short_vector(lambda, options);
 241         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 242             if (lambda[i] == 0)
 243                 continue;
 244             Matrix* M = Matrix_Copy(c->Rays);
 245             Vector_Copy(v->p, M->p[i], v->Size);
 246             cone * pc = new cone(M);
 247             assert (pc->det != 0);
 248             if (abs(pc->det) > 1) {
 249                 assert(abs(pc->det) < abs(c->det));
 250                 nonuni.push_back(pc);
 251             } else {
 252                 Polyhedron *p = pc->poly();
 253                 pc->Cone = 0;
 254                 if (sign(pc->det) == s) {
 255                     p->next = pos;
 256                     pos = p;
 257                 } else {
 258                     p->next = neg;
 259                     neg = p;
 260                 }
 261                 delete pc;
 262             }
 263             Matrix_Free(M);
 264         }
 265         Vector_Free(v);
 266         delete c;
 267     }
 268     *ppos = pos;
 269     *pneg = neg;
 270     free(options);
 271 }
 272