summate.c

   1 #include <barvinok/options.h>
   2 #include "bernoulli.h"
   3 #include "euler.h"
   4 #include "laurent.h"
   5 #include "summate.h"
   6 #include "section_array.h"
   7
   8 #define ALLOC(type) (type*)malloc(sizeof(type))
   9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  10
  11 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
  12                                  unsigned nvar,
  13                                  struct evalue_section_array *sections,
  14                                  struct barvinok_options *options)
  15 {
  16     unsigned dim = P->Dimension;
  17     unsigned new_dim, new_nparam;
  18     Matrix *T = NULL, *CP = NULL;
  19     evalue **subs;
  20     evalue *sum;
  21     int j;
  22
  23     if (emptyQ(P))
  24         return evalue_zero();
  25
  26     assert(P->NbEq > 0);
  27
  28     remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
  29
  30     if (emptyQ(P)) {
  31         Polyhedron_Free(P);
  32         return evalue_zero();
  33     }
  34
  35     new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
  36     new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
  37
  38     /* We can avoid these substitutions if E is a constant */
  39     subs = ALLOCN(evalue *, dim);
  40     for (j = 0; j < nvar; ++j) {
  41         if (T)
  42             subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim],
  43                                     new_dim);
  44         else
  45             subs[j] = evalue_var(j);
  46     }
  47     for (j = 0; j < dim-nvar; ++j) {
  48         if (CP)
  49             subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam],
  50                                          new_nparam);
  51         else
  52             subs[nvar+j] = evalue_var(j);
  53         evalue_shift_variables(subs[nvar+j], 0, new_dim-new_nparam);
  54     }
  55
  56     E = evalue_dup(E);
  57     evalue_substitute(E, subs);
  58     reduce_evalue(E);
  59
  60     for (j = 0; j < dim; ++j)
  61         evalue_free(subs[j]);
  62     free(subs);
  63
  64     if (new_dim-new_nparam > 0) {
  65         sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam,
  66                                          sections, options);
  67         evalue_free(E);
  68         Polyhedron_Free(P);
  69     } else {
  70         sum = ALLOC(evalue);
  71         value_init(sum->d);
  72         sum->x.p = new_enode(partition, 2, new_dim);
  73         EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
  74         value_clear(sum->x.p->arr[1].d);
  75         sum->x.p->arr[1] = *E;
  76         free(E);
  77     }
  78
  79     if (CP) {
  80         evalue_backsubstitute(sum, CP, options->MaxRays);
  81         Matrix_Free(CP);
  82     }
  83
  84     if (T)
  85         Matrix_Free(T);
  86
  87     return sum;
  88 }
  89
  90 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
  91                                      struct evalue_section_array *sections,
  92                                      struct barvinok_options *options)
  93 {
  94     if (P->NbEq)
  95         return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
  96
  97     if (options->summation == BV_SUM_EULER)
  98         return euler_summate(P, E, nvar, options);
  99     else if (options->summation == BV_SUM_LAURENT)
 100         return laurent_summate(P, E, nvar, options);
 101     else if (options->summation == BV_SUM_BERNOULLI)
 102         return bernoulli_summate(P, E, nvar, sections, options);
 103     else
 104         return box_summate(P, E, nvar, options->MaxRays);
 105 }
 106
 107 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
 108 {
 109     int i;
 110     struct evalue_section_array sections;
 111     evalue *sum;
 112
 113     assert(nvar >= 0);
 114     if (nvar == 0 || EVALUE_IS_ZERO(*e))
 115         return evalue_dup(e);
 116
 117     assert(value_zero_p(e->d));
 118     assert(e->x.p->type == partition);
 119
 120     evalue_section_array_init(&sections);
 121     sum = evalue_zero();
 122
 123     for (i = 0; i < e->x.p->size/2; ++i) {
 124         Polyhedron *D;
 125         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 126             Polyhedron *next = D->next;
 127             evalue *tmp;
 128             D->next = NULL;
 129
 130             tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
 131                                              &sections, options);
 132             assert(tmp);
 133             eadd(tmp, sum);
 134             evalue_free(tmp);
 135
 136             D->next = next;
 137         }
 138     }
 139
 140     free(sections.s);
 141
 142     reduce_evalue(sum);
 143     return sum;
 144 }
 145
 146 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
 147 {
 148     evalue *sum;
 149     struct barvinok_options *options = barvinok_options_new_with_defaults();
 150     options->MaxRays = MaxRays;
 151     sum = barvinok_summate(E, nvar, options);
 152     barvinok_options_free(options);
 153     return sum;
 154 }
 155
 156 evalue *esum(evalue *e, int nvar)
 157 {
 158     evalue *sum;
 159     struct barvinok_options *options = barvinok_options_new_with_defaults();
 160     sum = barvinok_summate(e, nvar, options);
 161     barvinok_options_free(options);
 162     return sum;
 163 }
 164
 165 /* Turn unweighted counting problem into "weighted" counting problem
 166  * with weight equal to 1 and call barvinok_summate on this weighted problem.
 167  */
 168 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
 169                                     struct barvinok_options *options)
 170 {
 171     Polyhedron *CA, *D;
 172     evalue e;
 173     evalue *sum;
 174
 175     if (emptyQ(P) || emptyQ(C))
 176         return evalue_zero();
 177
 178     CA = align_context(C, P->Dimension, options->MaxRays);
 179     D = DomainIntersection(P, CA, options->MaxRays);
 180     Domain_Free(CA);
 181
 182     if (emptyQ(D)) {
 183         Domain_Free(D);
 184         return evalue_zero();
 185     }
 186
 187     value_init(e.d);
 188     e.x.p = new_enode(partition, 2, P->Dimension);
 189     EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
 190     evalue_set_si(&e.x.p->arr[1], 1, 1);
 191     sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
 192     free_evalue_refs(&e);
 193     return sum;
 194 }