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(emptyQ2(P) || 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         POL_ENSURE_VERTICES(P);
  64         if (P->NbEq)
  65             gf = series(Polyhedron_Copy(P), nparam, options);
  66         else {
  67             gf_base *red;
  68             red = gf_base::create(Polyhedron_Project(P, nparam),
  69                                   P->Dimension, nparam, options);
  70             red->start_gf(P, options);
  71             gf = red->gf;
  72             delete red;
  73         }
  74     }
  75     if (CP) {
  76         gf->substitute(CP);
  77         Matrix_Free(CP);
  78     }
  79     Polyhedron_Free(P);
  80     return gf;
  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 = series(P, nparam, options);
  95
  96     return gf;
  97 }
  98
  99 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C, unsigned MaxRays)
 100 {
 101     gen_fun *gf;
 102     barvinok_options *options = barvinok_options_new_with_defaults();
 103     options->MaxRays = MaxRays;
 104     gf = barvinok_series_with_options(P, C, options);
 105     barvinok_options_free(options);
 106     return gf;
 107 }
 108
 109 static Polyhedron *skew_into_positive_orthant(Polyhedron *D, unsigned nparam,
 110                                               unsigned MaxRays)
 111 {
 112     Matrix *M = NULL;
 113     Value tmp;
 114     value_init(tmp);
 115     for (Polyhedron *P = D; P; P = P->next) {
 116         POL_ENSURE_VERTICES(P);
 117         assert(!Polyhedron_is_unbounded(P, nparam, MaxRays));
 118         assert(P->NbBid == 0);
 119         assert(Polyhedron_has_positive_rays(P, nparam));
 120
 121         for (int r = 0; r < P->NbRays; ++r) {
 122             if (value_notzero_p(P->Ray[r][P->Dimension+1]))
 123                 continue;
 124             for (int i = 0; i < nparam; ++i) {
 125                 int j;
 126                 if (value_posz_p(P->Ray[r][i+1]))
 127                     continue;
 128                 if (!M) {
 129                     M = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 130                     for (int i = 0; i < D->Dimension+1; ++i)
 131                         value_set_si(M->p[i][i], 1);
 132                 } else {
 133                     Inner_Product(P->Ray[r]+1, M->p[i], D->Dimension+1, &tmp);
 134                     if (value_posz_p(tmp))
 135                         continue;
 136                 }
 137                 for (j = P->Dimension - nparam; j < P->Dimension; ++j)
 138                     if (value_pos_p(P->Ray[r][j+1]))
 139                         break;
 140                 assert(j < P->Dimension);
 141                 value_pdivision(tmp, P->Ray[r][j+1], P->Ray[r][i+1]);
 142                 value_subtract(M->p[i][j], M->p[i][j], tmp);
 143             }
 144         }
 145     }
 146     value_clear(tmp);
 147     if (M) {
 148         D = DomainImage(D, M, MaxRays);
 149         Matrix_Free(M);
 150     }
 151     return D;
 152 }
 153
 154 gen_fun* barvinok_enumerate_union_series_with_options(Polyhedron *D, Polyhedron* C,
 155                                                       barvinok_options *options)
 156 {
 157     Polyhedron *conv, *D2;
 158     Polyhedron *CA;
 159     gen_fun *gf = NULL, *gf2;
 160     unsigned nparam = C->Dimension;
 161     ZZ one, mone;
 162     one = 1;
 163     mone = -1;
 164
 165     CA = align_context(C, D->Dimension, options->MaxRays);
 166     D = DomainIntersection(D, CA, options->MaxRays);
 167     Polyhedron_Free(CA);
 168
 169     D2 = skew_into_positive_orthant(D, nparam, options->MaxRays);
 170     for (Polyhedron *P = D2; P; P = P->next) {
 171         assert(P->Dimension == D2->Dimension);
 172         gen_fun *P_gf;
 173
 174         P_gf = series(Polyhedron_Copy(P), P->Dimension, options);
 175         if (!gf)
 176             gf = P_gf;
 177         else {
 178             gf->add_union(P_gf, options);
 179             delete P_gf;
 180         }
 181     }
 182     /* we actually only need the convex union of the parameter space
 183      * but the reducer classes currently expect a polyhedron in
 184      * the combined space
 185      */
 186     Polyhedron_Free(gf->context);
 187     gf->context = DomainConvex(D2, options->MaxRays);
 188
 189     gf2 = gf->summate(D2->Dimension - nparam, options);
 190
 191     delete gf;
 192     if (D != D2)
 193         Domain_Free(D2);
 194     Domain_Free(D);
 195     return gf2;
 196 }
 197
 198 gen_fun* barvinok_enumerate_union_series(Polyhedron *D, Polyhedron* C,
 199                                          unsigned MaxRays)
 200 {
 201     gen_fun *gf;
 202     barvinok_options *options = barvinok_options_new_with_defaults();
 203     options->MaxRays = MaxRays;
 204     gf = barvinok_enumerate_union_series_with_options(D, C, options);
 205     barvinok_options_free(options);
 206     return gf;
 207 }