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