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 #include "reduce_domain.h"
  11
  12 #ifdef NTL_STD_CXX
  13 using namespace NTL;
  14 #endif
  15 using std::vector;
  16 using std::cerr;
  17 using std::endl;
  18
  19 /*
  20  * Returns the largest absolute value in the vector
  21  */
  22 static ZZ max(vec_ZZ& v)
  23 {
  24     ZZ max = abs(v[0]);
  25     for (int i = 1; i < v.length(); ++i)
  26         if (abs(v[i]) > max)
  27             max = abs(v[i]);
  28     return max;
  29 }
  30
  31 /* Remove common divisor of elements of cols of B */
  32 static void normalize_cols(mat_ZZ& B)
  33 {
  34     ZZ gcd;
  35     for (int i = 0; i < B.NumCols(); ++i) {
  36         gcd = B[0][i];
  37         for (int j = 1 ; gcd != 1 && j < B.NumRows(); ++j)
  38             GCD(gcd, gcd, B[j][i]);
  39         if (gcd != 1)
  40             for (int j = 0; j < B.NumRows(); ++j)
  41                 B[j][i] /= gcd;
  42     }
  43 }
  44
  45 /* Remove common divisor of elements of B */
  46 static void normalize_matrix(mat_ZZ& B)
  47 {
  48     ZZ gcd;
  49     for (int i = 0; i < B.NumCols(); ++i)
  50         for (int j = 0 ; j < B.NumRows(); ++j) {
  51             GCD(gcd, gcd, B[i][j]);
  52             if (IsOne(gcd))
  53                 return;
  54         }
  55     for (int i = 0; i < B.NumCols(); ++i)
  56         for (int j = 0; j < B.NumRows(); ++j)
  57             B[i][j] /= gcd;
  58 }
  59
  60 class cone {
  61 public:
  62     cone(const mat_ZZ& r, int row, const vec_ZZ& w, int s) {
  63         sgn = s;
  64         rays = r;
  65         rays[row] = w;
  66         set_det();
  67     }
  68     cone(const signed_cone& sc) {
  69         rays = sc.rays;
  70         sgn = sc.sign;
  71         set_det();
  72     }
  73     void set_det() {
  74         det = determinant(rays);
  75         assert(!IsZero(det));
  76     }
  77     bool needs_split(barvinok_options *options) {
  78         index = abs(det);
  79         if (IsOne(index))
  80             return false;
  81         if (options->primal && index <= options->max_index)
  82             return false;
  83
  84         inv(det, B, rays);
  85         normalize_matrix(B);
  86         if (sign(det) < 0)
  87             negate(B, B);
  88
  89         if (!options->primal && options->max_index > 1) {
  90             mat_ZZ B2 = B;
  91             normalize_cols(B2);
  92             index = abs(determinant(B2));
  93             if (index <= options->max_index)
  94                 return false;
  95         }
  96
  97         return true;
  98     }
  99
 100     void short_vector(vec_ZZ& v, vec_ZZ& lambda, barvinok_options *options) {
 101         ZZ det2;
 102         mat_ZZ U;
 103         long r = LLL(det2, B, U, options->LLL_a, options->LLL_b);
 104
 105         ZZ min = max(B[0]);
 106         int index = 0;
 107         for (int i = 1; i < B.NumRows(); ++i) {
 108             ZZ tmp = max(B[i]);
 109             if (tmp < min) {
 110                 min = tmp;
 111                 index = i;
 112             }
 113         }
 114
 115         lambda = B[index];
 116
 117         v = U[index];
 118
 119         int i;
 120         for (i = 0; i < lambda.length(); ++i)
 121             if (lambda[i] > 0)
 122                 break;
 123         if (i == lambda.length()) {
 124             v = -v;
 125             lambda = -lambda;
 126         }
 127     }
 128
 129     ZZ det;
 130     ZZ index;
 131     mat_ZZ rays;
 132     mat_ZZ B;
 133     int sgn;
 134 };
 135
 136 std::ostream & operator<<(std::ostream & os, const cone& c)
 137 {
 138     os << c.rays << endl;
 139     os << "det: " << c.det << endl;
 140     os << "sign: " << c.sgn << endl;
 141     return os;
 142 }
 143
 144 static void decompose(const signed_cone& sc, signed_cone_consumer& scc,
 145                       bool primal, barvinok_options *options)
 146 {
 147     vector<cone *> nonuni;
 148     cone *c = new cone(sc);
 149     if (c->needs_split(options)) {
 150         nonuni.push_back(c);
 151     } else {
 152         try {
 153             options->stats->base_cones++;
 154             scc.handle(signed_cone(sc.C, sc.rays, sc.sign, to_ulong(c->index)),
 155                                    options);
 156             delete c;
 157         } catch (...) {
 158             delete c;
 159             throw;
 160         }
 161         return;
 162     }
 163     vec_ZZ lambda;
 164     vec_ZZ v;;
 165     while (!nonuni.empty()) {
 166         c = nonuni.back();
 167         nonuni.pop_back();
 168         c->short_vector(v, lambda, options);
 169         for (int i = 0; i < c->rays.NumRows(); ++i) {
 170             if (lambda[i] == 0)
 171                 continue;
 172             cone *pc = new cone(c->rays, i, v, sign(lambda[i]) * c->sgn);
 173             if (primal) {
 174                 for (int j = 0; j <= i; ++j) {
 175                     if ((j == i && sign(lambda[i]) < 0) ||
 176                         (j < i && sign(lambda[i]) == sign(lambda[j]))) {
 177                         pc->rays[j] = -pc->rays[j];
 178                         pc->sgn = -pc->sgn;
 179                     }
 180                 }
 181             }
 182             if (pc->needs_split(options)) {
 183                 assert(abs(pc->det) < abs(c->det));
 184                 nonuni.push_back(pc);
 185             } else {
 186                 try {
 187                     options->stats->base_cones++;
 188                     scc.handle(signed_cone(pc->rays, pc->sgn,
 189                                            to_ulong(pc->index)),
 190                                            options);
 191                     delete pc;
 192                 } catch (...) {
 193                     delete c;
 194                     delete pc;
 195                     while (!nonuni.empty()) {
 196                         c = nonuni.back();
 197                         nonuni.pop_back();
 198                         delete c;
 199                     }
 200                     throw;
 201                 }
 202             }
 203         }
 204         delete c;
 205     }
 206 }
 207
 208 struct polar_signed_cone_consumer : public signed_cone_consumer {
 209     signed_cone_consumer& scc;
 210     mat_ZZ r;
 211     polar_signed_cone_consumer(signed_cone_consumer& scc) : scc(scc) {}
 212     virtual void handle(const signed_cone& sc, barvinok_options *options) {
 213         Polyhedron *C = sc.C;
 214         if (!sc.C) {
 215             Matrix *M = Matrix_Alloc(sc.rays.NumRows()+1, sc.rays.NumCols()+2);
 216             for (int i = 0; i < sc.rays.NumRows(); ++i) {
 217                 zz2values(sc.rays[i], M->p[i]+1);
 218                 value_set_si(M->p[i][0], 1);
 219             }
 220             value_set_si(M->p[sc.rays.NumRows()][0], 1);
 221             value_set_si(M->p[sc.rays.NumRows()][1+sc.rays.NumCols()], 1);
 222             C = Rays2Polyhedron(M, M->NbRows+1);
 223             assert(C->NbConstraints == C->NbRays);
 224             Matrix_Free(M);
 225         }
 226         Polyhedron_Polarize(C);
 227         rays(C, r);
 228         scc.handle(signed_cone(C, r, sc.sign, sc.det), options);
 229         if (!sc.C)
 230             Polyhedron_Free(C);
 231     }
 232 };
 233
 234 /* Remove common divisor of elements of rows of B */
 235 static void normalize_rows(mat_ZZ& B)
 236 {
 237     ZZ gcd;
 238     for (int i = 0; i < B.NumRows(); ++i) {
 239         gcd = B[i][0];
 240         for (int j = 1 ; gcd != 1 && j < B.NumCols(); ++j)
 241             GCD(gcd, gcd, B[i][j]);
 242         if (gcd != 1)
 243             for (int j = 0; j < B.NumCols(); ++j)
 244                 B[i][j] /= gcd;
 245     }
 246 }
 247
 248 static void polar_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 249                             barvinok_options *options)
 250 {
 251     POL_ENSURE_VERTICES(cone);
 252     Polyhedron_Polarize(cone);
 253     if (cone->NbRays - 1 != cone->Dimension) {
 254         Polyhedron *tmp = cone;
 255         cone = triangulate_cone_with_options(cone, options);
 256         Polyhedron_Free(tmp);
 257     }
 258     polar_signed_cone_consumer pssc(scc);
 259     mat_ZZ r;
 260     try {
 261         for (Polyhedron *Polar = cone; Polar; Polar = Polar->next) {
 262             rays(Polar, r);
 263             normalize_rows(r);
 264             decompose(signed_cone(Polar, r, 1), pssc, false, options);
 265         }
 266         Domain_Free(cone);
 267     } catch (...) {
 268         Domain_Free(cone);
 269         throw;
 270     }
 271 }
 272
 273 static void primal_decompose(Polyhedron *cone, signed_cone_consumer& scc,
 274                              barvinok_options *options)
 275 {
 276     POL_ENSURE_VERTICES(cone);
 277     Polyhedron *parts;
 278     if (cone->NbRays - 1 == cone->Dimension)
 279         parts = cone;
 280     else
 281         parts = triangulate_cone_with_options(cone, options);
 282     Vector *average = NULL;
 283     Value tmp;
 284     if (parts != cone) {
 285         value_init(tmp);
 286         average = inner_point(cone);
 287     }
 288     mat_ZZ ray;
 289     try {
 290         for (Polyhedron *simple = parts; simple; simple = simple->next) {
 291             int sign = 1;
 292             Matrix *Rays = rays2(simple);
 293             for (int i = 0; i < Rays->NbRows; ++i) {
 294                 if (simple == cone) {
 295                     continue;
 296                 } else {
 297                     int f;
 298                     for (f = 0; f < simple->NbConstraints; ++f) {
 299                         Inner_Product(Rays->p[i], simple->Constraint[f]+1,
 300                                       simple->Dimension, &tmp);
 301                         if (value_notzero_p(tmp))
 302                             break;
 303                     }
 304                     assert(f < simple->NbConstraints);
 305                     if (!is_internal(average, simple->Constraint[f])) {
 306                         Vector_Oppose(Rays->p[i], Rays->p[i], Rays->NbColumns);
 307                         sign = -sign;
 308                     }
 309                 }
 310             }
 311             matrix2zz(Rays, ray, Rays->NbRows, Rays->NbColumns);
 312             Matrix_Free(Rays);
 313             decompose(signed_cone(simple, ray, sign), scc, true, options);
 314         }
 315         Domain_Free(parts);
 316         if (parts != cone) {
 317             Domain_Free(cone);
 318             value_clear(tmp);
 319             Vector_Free(average);
 320         }
 321     } catch (...) {
 322         Domain_Free(parts);
 323         if (parts != cone) {
 324             Domain_Free(cone);
 325             value_clear(tmp);
 326             Vector_Free(average);
 327         }
 328         throw;
 329     }
 330 }
 331
 332 void barvinok_decompose(Polyhedron *C, signed_cone_consumer& scc,
 333                         barvinok_options *options)
 334 {
 335     if (options->primal)
 336         primal_decompose(C, scc, options);
 337     else
 338         polar_decompose(C, scc, options);
 339 }
 340
 341 void vertex_decomposer::decompose_at_vertex(Param_Vertices *V, int _i,
 342                                             barvinok_options *options)
 343 {
 344     Polyhedron *C = supporting_cone_p(P, V);
 345     vert = _i;
 346     this->V = V;
 347
 348     barvinok_decompose(C, scc, options);
 349 }
 350
 351 struct posneg_collector : public signed_cone_consumer {
 352     posneg_collector(Polyhedron *pos, Polyhedron *neg) : pos(pos), neg(neg) {}
 353     virtual void handle(const signed_cone& sc, barvinok_options *options) {
 354         Polyhedron *p = Polyhedron_Copy(sc.C);
 355         if (sc.sign > 0) {
 356             p->next = pos;
 357             pos = p;
 358         } else {
 359             p->next = neg;
 360             neg = p;
 361         }
 362     }
 363     Polyhedron *pos;
 364     Polyhedron *neg;
 365 };
 366
 367 /*
 368  * Barvinok's Decomposition of a simplicial cone
 369  *
 370  * Returns two lists of polyhedra
 371  */
 372 void barvinok_decompose(Polyhedron *C, Polyhedron **ppos, Polyhedron **pneg)
 373 {
 374     barvinok_options *options = barvinok_options_new_with_defaults();
 375     posneg_collector pc(*ppos, *pneg);
 376     POL_ENSURE_VERTICES(C);
 377     mat_ZZ r;
 378     rays(C, r);
 379     decompose(signed_cone(C, r, 1), pc, false, options);
 380     *ppos = pc.pos;
 381     *pneg = pc.neg;
 382     barvinok_options_free(options);
 383 }
 384