decomposer.cc

   1 #include <vector>
   2 #include <gmp.h>
   3 #include <NTL/mat_ZZ.h>
   4 #include <NTL/LLL.h>
   5 extern "C" {
   6 #include <polylib/polylibgmp.h>
   7 }
   8 #include <barvinok/barvinok.h>
   9 #include <barvinok/util.h>
  10 #include "conversion.h"
  11 #include "decomposer.h"
  12
  13 #ifdef NTL_STD_CXX
  14 using namespace NTL;
  15 #endif
  16 using std::vector;
  17
  18 /*
  19  * Returns the largest absolute value in the vector
  20  */
  21 static ZZ max(vec_ZZ& v)
  22 {
  23     ZZ max = abs(v[0]);
  24     for (int i = 1; i < v.length(); ++i)
  25         if (abs(v[i]) > max)
  26             max = abs(v[i]);
  27     return max;
  28 }
  29
  30 class cone {
  31 public:
  32     cone(Matrix *M) {
  33         Cone = 0;
  34         Rays = Matrix_Copy(M);
  35         set_det();
  36     }
  37     cone(Polyhedron *C) {
  38         Cone = Polyhedron_Copy(C);
  39         Rays = rays(C);
  40         set_det();
  41     }
  42     void set_det() {
  43         mat_ZZ A;
  44         matrix2zz(Rays, A, Rays->NbRows - 1, Rays->NbColumns - 1);
  45         det = determinant(A);
  46     }
  47
  48     Vector* short_vector(vec_ZZ& lambda) {
  49         Matrix *M = Matrix_Copy(Rays);
  50         Matrix *inv = Matrix_Alloc(M->NbRows, M->NbColumns);
  51         int ok = Matrix_Inverse(M, inv);
  52         assert(ok);
  53         Matrix_Free(M);
  54
  55         ZZ det2;
  56         mat_ZZ B;
  57         mat_ZZ U;
  58         matrix2zz(inv, B, inv->NbRows - 1, inv->NbColumns - 1);
  59         long r = LLL(det2, B, U);
  60
  61         ZZ min = max(B[0]);
  62         int index = 0;
  63         for (int i = 1; i < B.NumRows(); ++i) {
  64             ZZ tmp = max(B[i]);
  65             if (tmp < min) {
  66                 min = tmp;
  67                 index = i;
  68             }
  69         }
  70
  71         Matrix_Free(inv);
  72
  73         lambda = B[index];
  74
  75         Vector *z = Vector_Alloc(U[index].length()+1);
  76         assert(z);
  77         zz2values(U[index], z->p);
  78         value_set_si(z->p[U[index].length()], 0);
  79
  80         Polyhedron *C = poly();
  81         int i;
  82         for (i = 0; i < lambda.length(); ++i)
  83             if (lambda[i] > 0)
  84                 break;
  85         if (i == lambda.length()) {
  86             Value tmp;
  87             value_init(tmp);
  88             value_set_si(tmp, -1);
  89             Vector_Scale(z->p, z->p, tmp, z->Size-1);
  90             value_clear(tmp);
  91         }
  92         return z;
  93     }
  94
  95     ~cone() {
  96         Polyhedron_Free(Cone);
  97         Matrix_Free(Rays);
  98     }
  99
 100     Polyhedron *poly() {
 101         if (!Cone) {
 102             Matrix *M = Matrix_Alloc(Rays->NbRows+1, Rays->NbColumns+1);
 103             for (int i = 0; i < Rays->NbRows; ++i) {
 104                 Vector_Copy(Rays->p[i], M->p[i]+1, Rays->NbColumns);
 105                 value_set_si(M->p[i][0], 1);
 106             }
 107             Vector_Set(M->p[Rays->NbRows]+1, 0, Rays->NbColumns-1);
 108             value_set_si(M->p[Rays->NbRows][0], 1);
 109             value_set_si(M->p[Rays->NbRows][Rays->NbColumns], 1);
 110             Cone = Rays2Polyhedron(M, M->NbRows+1);
 111             assert(Cone->NbConstraints == Cone->NbRays);
 112             Matrix_Free(M);
 113         }
 114         return Cone;
 115     }
 116
 117     ZZ det;
 118     Polyhedron *Cone;
 119     Matrix *Rays;
 120 };
 121
 122 void decomposer::decompose(Polyhedron *C)
 123 {
 124     vector<cone *> nonuni;
 125     cone * c = new cone(C);
 126     ZZ det = c->det;
 127     int s = sign(det);
 128     assert(det != 0);
 129     if (abs(det) > 1) {
 130         nonuni.push_back(c);
 131     } else {
 132         try {
 133             handle(C, 1);
 134             delete c;
 135         } catch (...) {
 136             delete c;
 137             throw;
 138         }
 139     }
 140     vec_ZZ lambda;
 141     while (!nonuni.empty()) {
 142         c = nonuni.back();
 143         nonuni.pop_back();
 144         Vector* v = c->short_vector(lambda);
 145         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 146             if (lambda[i] == 0)
 147                 continue;
 148             Matrix* M = Matrix_Copy(c->Rays);
 149             Vector_Copy(v->p, M->p[i], v->Size);
 150             cone * pc = new cone(M);
 151             assert (pc->det != 0);
 152             if (abs(pc->det) > 1) {
 153                 assert(abs(pc->det) < abs(c->det));
 154                 nonuni.push_back(pc);
 155             } else {
 156                 try {
 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, unsigned MaxRays)
 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, MaxRays);
 186         Polyhedron_Free(tmp);
 187     }
 188     try {
 189         for (Polyhedron *Polar = cone; Polar; Polar = Polar->next)
 190             decomposer::decompose(Polar);
 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                                             unsigned MaxRays)
 206 {
 207     Polyhedron *C = supporting_cone_p(P, V);
 208     vert = _i;
 209     this->V = V;
 210
 211     pd->decompose(C, MaxRays);
 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     Polyhedron *pos = *ppos, *neg = *pneg;
 222     vector<cone *> nonuni;
 223     cone * c = new cone(C);
 224     ZZ det = c->det;
 225     int s = sign(det);
 226     assert(det != 0);
 227     if (abs(det) > 1) {
 228         nonuni.push_back(c);
 229     } else {
 230         Polyhedron *p = Polyhedron_Copy(c->Cone);
 231         p->next = pos;
 232         pos = p;
 233         delete c;
 234     }
 235     vec_ZZ lambda;
 236     while (!nonuni.empty()) {
 237         c = nonuni.back();
 238         nonuni.pop_back();
 239         Vector* v = c->short_vector(lambda);
 240         for (int i = 0; i < c->Rays->NbRows - 1; ++i) {
 241             if (lambda[i] == 0)
 242                 continue;
 243             Matrix* M = Matrix_Copy(c->Rays);
 244             Vector_Copy(v->p, M->p[i], v->Size);
 245             cone * pc = new cone(M);
 246             assert (pc->det != 0);
 247             if (abs(pc->det) > 1) {
 248                 assert(abs(pc->det) < abs(c->det));
 249                 nonuni.push_back(pc);
 250             } else {
 251                 Polyhedron *p = pc->poly();
 252                 pc->Cone = 0;
 253                 if (sign(pc->det) == s) {
 254                     p->next = pos;
 255                     pos = p;
 256                 } else {
 257                     p->next = neg;
 258                     neg = p;
 259                 }
 260                 delete pc;
 261             }
 262             Matrix_Free(M);
 263         }
 264         Vector_Free(v);
 265         delete c;
 266     }
 267     *ppos = pos;
 268     *pneg = neg;
 269 }
 270