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