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