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