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 extern evalue *evalue_outer_floor(evalue *e);
   9 extern int evalue_replace_floor(evalue *e, const evalue *floor, int var);
  10 extern void evalue_drop_floor(evalue *e, const evalue *floor);
  11
  12 #define ALLOC(type) (type*)malloc(sizeof(type))
  13 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  14
  15 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
  16                                  unsigned nvar,
  17                                  struct evalue_section_array *sections,
  18                                  struct barvinok_options *options)
  19 {
  20     unsigned dim = P->Dimension;
  21     unsigned new_dim, new_nparam;
  22     Matrix *T = NULL, *CP = NULL;
  23     evalue **subs;
  24     evalue *sum;
  25     int j;
  26
  27     if (emptyQ(P))
  28         return evalue_zero();
  29
  30     assert(P->NbEq > 0);
  31
  32     remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
  33
  34     if (emptyQ(P)) {
  35         Polyhedron_Free(P);
  36         return evalue_zero();
  37     }
  38
  39     new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
  40     new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
  41
  42     /* We can avoid these substitutions if E is a constant */
  43     subs = ALLOCN(evalue *, dim);
  44     for (j = 0; j < nvar; ++j) {
  45         if (T)
  46             subs[j] = affine2evalue(T->p[j], T->p[nvar+new_nparam][new_dim],
  47                                     new_dim);
  48         else
  49             subs[j] = evalue_var(j);
  50     }
  51     for (j = 0; j < dim-nvar; ++j) {
  52         if (CP)
  53             subs[nvar+j] = affine2evalue(CP->p[j], CP->p[dim-nvar][new_nparam],
  54                                          new_nparam);
  55         else
  56             subs[nvar+j] = evalue_var(j);
  57         evalue_shift_variables(subs[nvar+j], 0, new_dim-new_nparam);
  58     }
  59
  60     E = evalue_dup(E);
  61     evalue_substitute(E, subs);
  62     reduce_evalue(E);
  63
  64     for (j = 0; j < dim; ++j)
  65         evalue_free(subs[j]);
  66     free(subs);
  67
  68     if (new_dim-new_nparam > 0) {
  69         sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam,
  70                                          sections, options);
  71         evalue_free(E);
  72         Polyhedron_Free(P);
  73     } else {
  74         sum = ALLOC(evalue);
  75         value_init(sum->d);
  76         sum->x.p = new_enode(partition, 2, new_dim);
  77         EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
  78         value_clear(sum->x.p->arr[1].d);
  79         sum->x.p->arr[1] = *E;
  80         free(E);
  81     }
  82
  83     if (CP) {
  84         evalue_backsubstitute(sum, CP, options->MaxRays);
  85         Matrix_Free(CP);
  86     }
  87
  88     if (T)
  89         Matrix_Free(T);
  90
  91     return sum;
  92 }
  93
  94 /* Add two constraints corresponding to floor = floor(e/d),
  95  *
  96  *       e - d t       >= 0
  97  *      -e + d t + d-1 >= 0
  98  *
  99  * e is assumed to be an affine expression.
 100  */
 101 Polyhedron *add_floor_var(Polyhedron *P, unsigned nvar, const evalue *floor,
 102                                      struct barvinok_options *options)
 103 {
 104     int i;
 105     unsigned dim = P->Dimension+1;
 106     Matrix *M = Matrix_Alloc(P->NbConstraints+2, 2+dim);
 107     Polyhedron *CP;
 108     Value *d = &M->p[0][1+nvar];
 109     evalue_extract_affine(floor, M->p[0]+1, M->p[0]+1+dim, d);
 110     value_oppose(*d, *d);
 111     value_set_si(M->p[0][0], 1);
 112     value_set_si(M->p[1][0], 1);
 113     Vector_Oppose(M->p[0]+1, M->p[1]+1, M->NbColumns-1);
 114     value_subtract(M->p[1][1+dim], M->p[1][1+dim], *d);
 115     value_decrement(M->p[1][1+dim], M->p[1][1+dim]);
 116
 117     for (i = 0; i < P->NbConstraints; ++i) {
 118         Vector_Copy(P->Constraint[i], M->p[i+2], 1+nvar);
 119         Vector_Copy(P->Constraint[i]+1+nvar, M->p[i+2]+1+nvar+1, dim-nvar-1+1);
 120     }
 121
 122     CP = Constraints2Polyhedron(M, options->MaxRays);
 123     Matrix_Free(M);
 124     return CP;
 125 }
 126
 127 static evalue *evalue_add(evalue *a, evalue *b)
 128 {
 129     if (!a)
 130         return b;
 131     if (!b)
 132         return a;
 133     eadd(a, b);
 134     evalue_free(a);
 135     return b;
 136 }
 137
 138 /* Compute sum of a step-polynomial over a polytope by grouping
 139  * terms containing the same floor-expressions and introducing
 140  * new variables for each such expression.
 141  * In particular, while there is any floor-expression left,
 142  * the step-polynomial is split into a polynomial containing
 143  * the expression, which is then converted to a new variable,
 144  * and a polynomial not containing the expression.
 145  */
 146 static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar,
 147                                      struct barvinok_options *options)
 148 {
 149     evalue *floor;
 150     evalue *cur = E;
 151     evalue *sum = NULL;
 152     evalue *t;
 153
 154     while ((floor = evalue_outer_floor(cur))) {
 155         Polyhedron *CP;
 156         evalue *converted = evalue_dup(cur);
 157         evalue *converted_floor = evalue_dup(floor);
 158
 159         evalue_shift_variables(converted, nvar, 1);
 160         evalue_shift_variables(converted_floor, nvar, 1);
 161         evalue_replace_floor(converted, converted_floor, nvar);
 162         CP = add_floor_var(P, nvar, converted_floor, options);
 163         evalue_free(converted_floor);
 164         t = sum_step_polynomial(CP, converted, nvar+1, options);
 165         evalue_free(converted);
 166         Polyhedron_Free(CP);
 167         sum = evalue_add(t, sum);
 168
 169         if (cur == E)
 170             cur = evalue_dup(cur);
 171         evalue_drop_floor(cur, floor);
 172         evalue_free(floor);
 173     }
 174
 175     if (EVALUE_IS_ZERO(*cur))
 176         t = NULL;
 177     else if (options->summation == BV_SUM_EULER)
 178         t = euler_summate(P, cur, nvar, options);
 179     else if (options->summation == BV_SUM_LAURENT)
 180         t = laurent_summate(P, cur, nvar, options);
 181     else
 182         assert(0);
 183
 184     if (E != cur)
 185         evalue_free(cur);
 186
 187     return evalue_add(t, sum);
 188 }
 189
 190 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
 191                                      struct evalue_section_array *sections,
 192                                      struct barvinok_options *options)
 193 {
 194     if (P->NbEq)
 195         return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
 196
 197     if (options->summation == BV_SUM_BERNOULLI)
 198         return bernoulli_summate(P, E, nvar, sections, options);
 199     else if (options->summation == BV_SUM_BARVINOK)
 200         return box_summate(P, E, nvar, options->MaxRays);
 201
 202     evalue_frac2floor2(E, 0);
 203
 204     return sum_step_polynomial(P, E, nvar, options);
 205 }
 206
 207 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
 208 {
 209     int i;
 210     struct evalue_section_array sections;
 211     evalue *sum;
 212
 213     assert(nvar >= 0);
 214     if (nvar == 0 || EVALUE_IS_ZERO(*e))
 215         return evalue_dup(e);
 216
 217     assert(value_zero_p(e->d));
 218     assert(e->x.p->type == partition);
 219
 220     evalue_section_array_init(&sections);
 221     sum = evalue_zero();
 222
 223     for (i = 0; i < e->x.p->size/2; ++i) {
 224         Polyhedron *D;
 225         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 226             Polyhedron *next = D->next;
 227             evalue *tmp;
 228             D->next = NULL;
 229
 230             tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
 231                                              &sections, options);
 232             assert(tmp);
 233             eadd(tmp, sum);
 234             evalue_free(tmp);
 235
 236             D->next = next;
 237         }
 238     }
 239
 240     free(sections.s);
 241
 242     reduce_evalue(sum);
 243     return sum;
 244 }
 245
 246 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
 247 {
 248     evalue *sum;
 249     struct barvinok_options *options = barvinok_options_new_with_defaults();
 250     options->MaxRays = MaxRays;
 251     sum = barvinok_summate(E, nvar, options);
 252     barvinok_options_free(options);
 253     return sum;
 254 }
 255
 256 evalue *esum(evalue *e, int nvar)
 257 {
 258     evalue *sum;
 259     struct barvinok_options *options = barvinok_options_new_with_defaults();
 260     sum = barvinok_summate(e, nvar, options);
 261     barvinok_options_free(options);
 262     return sum;
 263 }
 264
 265 /* Turn unweighted counting problem into "weighted" counting problem
 266  * with weight equal to 1 and call barvinok_summate on this weighted problem.
 267  */
 268 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
 269                                     struct barvinok_options *options)
 270 {
 271     Polyhedron *CA, *D;
 272     evalue e;
 273     evalue *sum;
 274
 275     if (emptyQ(P) || emptyQ(C))
 276         return evalue_zero();
 277
 278     CA = align_context(C, P->Dimension, options->MaxRays);
 279     D = DomainIntersection(P, CA, options->MaxRays);
 280     Domain_Free(CA);
 281
 282     if (emptyQ(D)) {
 283         Domain_Free(D);
 284         return evalue_zero();
 285     }
 286
 287     value_init(e.d);
 288     e.x.p = new_enode(partition, 2, P->Dimension);
 289     EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
 290     evalue_set_si(&e.x.p->arr[1], 1, 1);
 291     sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
 292     free_evalue_refs(&e);
 293     return sum;
 294 }