series.cc

   1 #include <assert.h>
   2 #include <NTL/vec_ZZ.h>
   3 #include <barvinok/polylib.h>
   4 #include <barvinok/barvinok.h>
   5 #include <barvinok/util.h>
   6 #include "genfun_constructor.h"
   7 #include "lattice_width.h"
   8 #include "remove_equalities.h"
   9
  10 using std::cerr;
  11 using std::endl;
  12
  13 static gen_fun *enumerate_series(Polyhedron *P, unsigned nparam,
  14                                     barvinok_options *options)
  15 {
  16     Matrix *CP = NULL;
  17     gen_fun *gf;
  18     Polyhedron *P_orig = P;
  19
  20     if (emptyQ2(P))
  21         return new gen_fun(Empty_Polyhedron(nparam));
  22
  23     if (P->NbEq != 0)
  24         remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
  25     assert(emptyQ2(P) || P->NbEq == 0);
  26     if (CP)
  27         nparam = CP->NbColumns-1;
  28
  29     if (nparam == 0) {
  30         Value c;
  31         value_init(c);
  32         barvinok_count_with_options(P, &c, options);
  33         gf = new gen_fun(c);
  34         value_clear(c);
  35     } else {
  36         POL_ENSURE_VERTICES(P);
  37         if (P->NbEq)
  38             gf = enumerate_series(P, nparam, options);
  39         else {
  40             gf_base *red;
  41             red = gf_base::create(Polyhedron_Project(P, nparam),
  42                                   P->Dimension, nparam, options);
  43             red->start_gf(P, options);
  44             gf = red->gf;
  45             delete red;
  46         }
  47     }
  48     if (CP) {
  49         gf->substitute(CP);
  50         Matrix_Free(CP);
  51     }
  52     if (P != P_orig)
  53         Polyhedron_Free(P);
  54     return gf;
  55 }
  56
  57 gen_fun *barvinok_enumerate_series(Polyhedron *P, unsigned nparam,
  58                                     barvinok_options *options)
  59 {
  60     if (emptyQ2(P))
  61         return new gen_fun(Empty_Polyhedron(nparam));
  62
  63     assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
  64     assert(P->NbBid == 0);
  65     assert(Polyhedron_has_revlex_positive_rays(P, nparam));
  66     return enumerate_series(P, nparam, options);
  67 }
  68
  69 gen_fun * barvinok_series_with_options(Polyhedron *P, Polyhedron* C,
  70                                        barvinok_options *options)
  71 {
  72     Polyhedron *CA;
  73     unsigned nparam = C->Dimension;
  74     gen_fun *gf;
  75
  76     CA = align_context(C, P->Dimension, options->MaxRays);
  77     P = DomainIntersection(P, CA, options->MaxRays);
  78     Polyhedron_Free(CA);
  79
  80     gf = barvinok_enumerate_series(P, nparam, options);
  81     Polyhedron_Free(P);
  82
  83     return gf;
  84 }
  85
  86 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
  87 {
  88     gen_fun *gf;
  89     barvinok_options *options = barvinok_options_new_with_defaults();
  90     options->MaxRays = MaxRays;
  91     gf = barvinok_series_with_options(P, C, options);
  92     barvinok_options_free(options);
  93     return gf;
  94 }
  95
  96 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
  97                                                       barvinok_options *options)
  98 {
  99     Polyhedron *CA;
 100     gen_fun *gf = NULL, *gf2;
 101     unsigned nparam = C->Dimension;
 102
 103     CA = align_context(C, D->Dimension, options->MaxRays);
 104     D = DomainIntersection(D, CA, options->MaxRays);
 105     Polyhedron_Free(CA);
 106
 107     for (Polyhedron *P = D; P; P = P->next) {
 108         assert(P->Dimension == D->Dimension);
 109         gen_fun *P_gf;
 110
 111         P_gf = barvinok_enumerate_series(P, P->Dimension, options);
 112         if (!gf)
 113             gf = P_gf;
 114         else {
 115             gf->add_union(P_gf, options);
 116             delete P_gf;
 117         }
 118     }
 119     /* we actually only need the convex union of the parameter space
 120      * but the reducer classes currently expect a polyhedron in
 121      * the combined space
 122      */
 123     Polyhedron_Free(gf->context);
 124     gf->context = DomainConvex(D, options->MaxRays);
 125
 126     gf2 = gf->summate(D->Dimension - nparam, options);
 127
 128     delete gf;
 129     Domain_Free(D);
 130     return gf2;
 131 }
 132
 133 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
 134                                          unsigned MaxRays)
 135 {
 136     gen_fun *gf;
 137     barvinok_options *options = barvinok_options_new_with_defaults();
 138     options->MaxRays = MaxRays;
 139     gf = barvinok_enumerate_union_series_with_options(D, C, options);
 140     barvinok_options_free(options);
 141     return gf;
 142 }
 143
 144 /* Unimodularly transform the polyhedron P, such that
 145  * the direction specified by dir corresponds to the last
 146  * variable in the transformed polyhedron.
 147  * The number of variables is given by the length of dir.
 148  */
 149 static Polyhedron *put_direction_last(Polyhedron *P, Vector *dir,
 150                                       unsigned MaxRays)
 151 {
 152     Polyhedron *R;
 153     Matrix *T;
 154     int n = dir->Size;
 155
 156     T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 157     T->NbColumns = T->NbRows = n;
 158     Vector_Copy(dir->p, T->p[0], n);
 159     unimodular_complete(T, 1);
 160     Vector_Exchange(T->p[0], T->p[n-1], n);
 161     T->NbColumns = T->NbRows = P->Dimension+1;
 162     for (int j = n; j < P->Dimension+1; ++j)
 163         value_set_si(T->p[j][j], 1);
 164
 165     R = Polyhedron_Image(P, T, MaxRays);
 166     Matrix_Free(T);
 167
 168     return R;
 169 }
 170
 171 /* Do we need to continue shifting and subtracting?
 172  * i is the number of times we shifted so far
 173  * n is the number of coordinates projected out
 174  */
 175 static bool more_shifts_needed(int j, int n,
 176                                 gen_fun *S, gen_fun *S_divide, const vec_ZZ& up,
 177                                 barvinok_options *options)
 178 {
 179     bool empty;
 180     gen_fun *hp;
 181
 182     /* For the 2-dimensional case, we need to subtract at most once */
 183     if (n == 2 && j > 0)
 184         return false;
 185
 186     S_divide->shift(up);
 187
 188     /* Assume that we have to subtract at least once */
 189     if (j == 0)
 190         return true;
 191
 192     hp = S->Hadamard_product(S_divide, options);
 193
 194     empty = hp->is_zero();
 195     delete hp;
 196
 197     return !empty;
 198 }
 199
 200 static gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options,
 201                         int free_P);
 202
 203 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
 204  * project out the first n coordinates.
 205  *
 206  * Assumes P has no equalities.
 207  */
 208 static gen_fun *project_full_dim(Polyhedron *P, unsigned n,
 209         barvinok_options *options)
 210 {
 211     QQ one(1, 1);
 212     QQ mone(-1, 1);
 213     vec_ZZ up;
 214     gen_fun *gf = NULL;
 215     struct width_direction_array *dirs;
 216     Polyhedron *U;
 217
 218     if (n == 0)
 219         return barvinok_enumerate_series(P, P->Dimension, options);
 220
 221     up.SetLength(P->Dimension - (n-1));
 222     up[0] = 1;
 223     for (int i = 1; i < P->Dimension - (n-1); ++i)
 224         up[i] = 0;
 225
 226     if (n == 1) {
 227         gen_fun *S, *S_shift, *hp;
 228
 229         S = barvinok_enumerate_series(P, P->Dimension, options);
 230         S_shift = new gen_fun(S);
 231         S_shift->shift(up);
 232         hp = S->Hadamard_product(S_shift, options);
 233         S->add(mone, hp, options);
 234         delete S_shift;
 235         delete hp;
 236
 237         gf = S->summate(1, options);
 238         delete S;
 239
 240         return gf;
 241     }
 242
 243     U = Universe_Polyhedron(P->Dimension - n);
 244     dirs = Polyhedron_Lattice_Width_Directions(P, U, options);
 245     Polyhedron_Free(U);
 246
 247     for (int i = 0; i < dirs->n; ++i) {
 248         Polyhedron *Pi, *R;
 249         Polyhedron *CA;
 250         gen_fun *S, *S_shift, *S_divide, *sum;
 251
 252         CA = align_context(dirs->wd[i].domain, P->Dimension, options->MaxRays);
 253         R = DomainIntersection(P, CA, options->MaxRays);
 254         Polyhedron_Free(CA);
 255         assert(dirs->wd[i].dir->Size == n);
 256         Pi = put_direction_last(R, dirs->wd[i].dir, options->MaxRays);
 257         Polyhedron_Free(R);
 258
 259         S = project(Pi, n-1, options, 1);
 260
 261         S_shift = new gen_fun(S);
 262         S_divide = new gen_fun(S);
 263         S_divide->divide(up);
 264
 265         for (int j = 0; more_shifts_needed(j, n, S, S_divide, up, options); ++j) {
 266             gen_fun *hp;
 267
 268             S_shift->shift(up);
 269             hp = S->Hadamard_product(S_shift, options);
 270             S->add(mone, hp, options);
 271
 272             delete hp;
 273         }
 274
 275         sum = S->summate(1, options);
 276
 277         delete S_shift;
 278         delete S_divide;
 279         delete S;
 280
 281         if (!gf)
 282             gf = sum;
 283         else {
 284             gf->add(one, sum, options);
 285             delete sum;
 286         }
 287     }
 288     free_width_direction_array(dirs);
 289
 290     return gf;
 291 }
 292
 293 /* Return gf of P projected on last dim(P)-n coordinates, i.e.,
 294  * project out the first n coordinates.
 295  */
 296 static gen_fun *project(Polyhedron *P, unsigned n, barvinok_options *options,
 297                         int free_P)
 298 {
 299     Matrix *CP = NULL;
 300     gen_fun *proj;
 301     unsigned nparam = P->Dimension - n;
 302
 303     if (P->NbEq != 0) {
 304         Polyhedron *Q = P;
 305         remove_all_equalities(&P, NULL, &CP, NULL, nparam, options->MaxRays);
 306         if (CP)
 307             nparam = CP->NbColumns - 1;
 308         n = P->Dimension - nparam;
 309         if (free_P)
 310             Polyhedron_Free(Q);
 311         free_P = 1;
 312     }
 313
 314     if (emptyQ(P))
 315         proj = new gen_fun(Empty_Polyhedron(nparam));
 316     else
 317         proj = project_full_dim(P, n, options);
 318     if (CP) {
 319         proj->substitute(CP);
 320         Matrix_Free(CP);
 321     }
 322
 323     if (free_P)
 324         Polyhedron_Free(P);
 325
 326     return proj;
 327 }
 328
 329 gen_fun *barvinok_enumerate_e_series(Polyhedron *P,
 330                   unsigned exist, unsigned nparam, barvinok_options *options)
 331 {
 332     Polyhedron *P_orig = P;
 333     gen_fun *gf, *proj;
 334     unsigned nvar = P->Dimension - exist - nparam;
 335
 336     if (exist == 0)
 337         return barvinok_enumerate_series(P, nparam, options);
 338
 339     if (emptyQ(P))
 340         return new gen_fun(Empty_Polyhedron(nparam));
 341
 342     assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
 343     assert(P->NbBid == 0);
 344     assert(Polyhedron_has_revlex_positive_rays(P, nparam));
 345
 346     /* Move existentially quantified variables to the front.*/
 347     if (nvar) {
 348         Matrix *T;
 349         T = Matrix_Alloc(exist+nvar+nparam+1, nvar+exist+nparam+1);
 350         for (int i = 0; i < exist; ++i)
 351             value_set_si(T->p[i][nvar+i], 1);
 352         for (int i = 0; i < nvar; ++i)
 353             value_set_si(T->p[exist+i][i], 1);
 354         for (int i = 0; i < nparam+1; ++i)
 355             value_set_si(T->p[exist+nvar+i][nvar+exist+i], 1);
 356         P = Polyhedron_Image(P, T, options->MaxRays);
 357         Matrix_Free(T);
 358     }
 359     proj = project(P, exist, options, P != P_orig);
 360
 361     if (!nvar)
 362         gf = proj;
 363     else {
 364         gf = proj->summate(nvar, options);
 365         delete proj;
 366     }
 367     return gf;
 368 }