series.cc

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