series.cc

   1 #include <assert.h>
   2 #include <barvinok/barvinok.h>
   3 #include <barvinok/util.h>
   4 #include "genfun_constructor.h"
   5 #include "remove_equalities.h"
   6
   7 /* Check whether all rays point in the positive directions
   8  * for the parameters
   9  */
  10 static bool Polyhedron_has_positive_rays(Polyhedron *P, unsigned nparam)
  11 {
  12     int r;
  13     for (r = 0; r < P->NbRays; ++r)
  14         if (value_zero_p(P->Ray[r][P->Dimension+1])) {
  15             int i;
  16             for (i = P->Dimension - nparam; i < P->Dimension; ++i)
  17                 if (value_neg_p(P->Ray[r][i+1]))
  18                     return false;
  19         }
  20     return true;
  21 }
  22
  23 /*
  24  * remove equalities that require a "compression" of the parameters
  25  */
  26 static Polyhedron *remove_more_equalities(Polyhedron *P, unsigned nparam,
  27                                           Matrix **CP, unsigned MaxRays)
  28 {
  29     Polyhedron *Q = P;
  30     remove_all_equalities(&P, NULL, CP, NULL, nparam, MaxRays);
  31     if (P != Q)
  32         Polyhedron_Free(Q);
  33     return P;
  34 }
  35
  36 /* frees P */
  37 static gen_fun *series(Polyhedron *P, unsigned nparam, barvinok_options *options)
  38 {
  39     Matrix *CP = NULL;
  40     gen_fun *gf;
  41
  42     if (emptyQ2(P)) {
  43         Polyhedron_Free(P);
  44         return new gen_fun;
  45     }
  46
  47     assert(!Polyhedron_is_unbounded(P, nparam, options->MaxRays));
  48     assert(P->NbBid == 0);
  49     assert(Polyhedron_has_revlex_positive_rays(P, nparam));
  50     if (P->NbEq != 0)
  51         P = remove_more_equalities(P, nparam, &CP, options->MaxRays);
  52     assert(P->NbEq == 0);
  53     if (CP)
  54         nparam = CP->NbColumns-1;
  55
  56     if (nparam == 0) {
  57         Value c;
  58         value_init(c);
  59         barvinok_count_with_options(P, &c, options);
  60         gf = new gen_fun(c);
  61         value_clear(c);
  62     } else {
  63         gf_base *red;
  64         red = gf_base::create(Polyhedron_Project(P, nparam),
  65                               P->Dimension, nparam, options);
  66         POL_ENSURE_VERTICES(P);
  67         red->start_gf(P, options);
  68         gf = red->gf;
  69         delete red;
  70     }
  71     if (CP) {
  72         gf->substitute(CP);
  73         Matrix_Free(CP);
  74     }
  75     Polyhedron_Free(P);
  76     return gf;
  77 }
  78
  79 gen_fun * barvinok_series_with_options(Polyhedron *P, Polyhedron* C,
  80                                        barvinok_options *options)
  81 {
  82     Polyhedron *CA;
  83     unsigned nparam = C->Dimension;
  84     gen_fun *gf;
  85
  86     CA = align_context(C, P->Dimension, options->MaxRays);
  87     P = DomainIntersection(P, CA, options->MaxRays);
  88     Polyhedron_Free(CA);
  89
  90     gf = series(P, nparam, options);
  91
  92     return gf;
  93 }
  94
  95 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
  96 {
  97     gen_fun *gf;
  98     barvinok_options *options = barvinok_options_new_with_defaults();
  99     options->MaxRays = MaxRays;
 100     gf = barvinok_series_with_options(P, C, options);
 101     barvinok_options_free(options);
 102     return gf;
 103 }
 104
 105 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
 106                                               unsigned MaxRays)
 107 {
 108     Matrix *M = NULL;
 109     Value tmp;
 110     value_init(tmp);
 111     for (Polyhedron *P = D; P; P = P->next) {
 112         POL_ENSURE_VERTICES(P);
 113         assert(!Polyhedron_is_unbounded(P, nparam, MaxRays));
 114         assert(P->NbBid == 0);
 115         assert(Polyhedron_has_positive_rays(P, nparam));
 116
 117         for (int r = 0; r < P->NbRays; ++r) {
 118             if (value_notzero_p(P->Ray[r][P->Dimension+1]))
 119                 continue;
 120             for (int i = 0; i < nparam; ++i) {
 121                 int j;
 122                 if (value_posz_p(P->Ray[r][i+1]))
 123                     continue;
 124                 if (!M) {
 125                     M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 126                     for (int i = 0; i < D->Dimension+1; ++i)
 127                         value_set_si(M->p[i][i], 1);
 128                 } else {
 129                     Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
 130                     if (value_posz_p(tmp))
 131                         continue;
 132                 }
 133                 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
 134                     if (value_pos_p(P->Ray[r][j+1]))
 135                         break;
 136                 assert(j < P->Dimension);
 137                 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
 138                 value_subtract(M->p[i][j], M->p[i][j], tmp);
 139             }
 140         }
 141     }
 142     value_clear(tmp);
 143     if (M) {
 144         D = DomainImage(D, M, MaxRays);
 145         Matrix_Free(M);
 146     }
 147     return D;
 148 }
 149
 150 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
 151                                                       barvinok_options *options)
 152 {
 153     Polyhedron *conv, *D2;
 154     Polyhedron *CA;
 155     gen_fun *gf = NULL, *gf2;
 156     unsigned nparam = C->Dimension;
 157     ZZ one, mone;
 158     one = 1;
 159     mone = -1;
 160
 161     CA = align_context(C, D->Dimension, options->MaxRays);
 162     D = DomainIntersection(D, CA, options->MaxRays);
 163     Polyhedron_Free(CA);
 164
 165     D2 = skew_into_positive_orthant(D, nparam, options->MaxRays);
 166     for (Polyhedron *P = D2; P; P = P->next) {
 167         assert(P->Dimension == D2->Dimension);
 168         gen_fun *P_gf;
 169
 170         P_gf = series(Polyhedron_Copy(P), P->Dimension, options);
 171         if (!gf)
 172             gf = P_gf;
 173         else {
 174             gf->add_union(P_gf, options);
 175             delete P_gf;
 176         }
 177     }
 178     /* we actually only need the convex union of the parameter space
 179      * but the reducer classes currently expect a polyhedron in
 180      * the combined space
 181      */
 182     Polyhedron_Free(gf->context);
 183     gf->context = DomainConvex(D2, options->MaxRays);
 184
 185     gf2 = gf->summate(D2->Dimension - nparam, options);
 186
 187     delete gf;
 188     if (D != D2)
 189         Domain_Free(D2);
 190     Domain_Free(D);
 191     return gf2;
 192 }
 193
 194 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
 195                                          unsigned MaxRays)
 196 {
 197     gen_fun *gf;
 198     barvinok_options *options = barvinok_options_new_with_defaults();
 199     options->MaxRays = MaxRays;
 200     gf = barvinok_enumerate_union_series_with_options(D, C, options);
 201     barvinok_options_free(options);
 202     return gf;
 203 }