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 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
   9                                  struct evalue_section_array *sections,
  10                                  struct barvinok_options *options)
  11 {
  12     if (options->summation == BV_SUM_EULER)
  13         return euler_summate(P, E, nvar, options);
  14     else if (options->summation == BV_SUM_LAURENT)
  15         return laurent_summate(P, E, nvar, options);
  16     else if (options->summation == BV_SUM_BERNOULLI)
  17         return bernoulli_summate(P, E, nvar, sections, options);
  18     else
  19         return box_summate(P, E, nvar, options->MaxRays);
  20 }
  21
  22 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
  23 {
  24     int i;
  25     struct evalue_section_array sections;
  26     evalue *sum;
  27
  28     assert(nvar >= 0);
  29     if (nvar == 0 || EVALUE_IS_ZERO(*e))
  30         return evalue_dup(e);
  31
  32     assert(value_zero_p(e->d));
  33     assert(e->x.p->type == partition);
  34
  35     evalue_section_array_init(&sections);
  36     sum = evalue_zero();
  37
  38     for (i = 0; i < e->x.p->size/2; ++i) {
  39         Polyhedron *D;
  40         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
  41             Polyhedron *next = D->next;
  42             evalue *tmp;
  43             D->next = NULL;
  44
  45             tmp = sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
  46                                         &sections, options);
  47             assert(tmp);
  48             eadd(tmp, sum);
  49             evalue_free(tmp);
  50
  51             D->next = next;
  52         }
  53     }
  54
  55     free(sections.s);
  56
  57     reduce_evalue(sum);
  58     return sum;
  59 }
  60
  61 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
  62 {
  63     evalue *sum;
  64     struct barvinok_options *options = barvinok_options_new_with_defaults();
  65     options->MaxRays = MaxRays;
  66     sum = barvinok_summate(E, nvar, options);
  67     barvinok_options_free(options);
  68     return sum;
  69 }
  70
  71 evalue *esum(evalue *e, int nvar)
  72 {
  73     evalue *sum;
  74     struct barvinok_options *options = barvinok_options_new_with_defaults();
  75     sum = barvinok_summate(e, nvar, options);
  76     barvinok_options_free(options);
  77     return sum;
  78 }
  79
  80 /* Turn unweighted counting problem into "weighted" counting problem
  81  * with weight equal to 1 and call barvinok_summate on this weighted problem.
  82  */
  83 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
  84                                     struct barvinok_options *options)
  85 {
  86     Polyhedron *CA, *D;
  87     evalue e;
  88     evalue *sum;
  89
  90     if (emptyQ(P) || emptyQ(C))
  91         return evalue_zero();
  92
  93     CA = align_context(C, P->Dimension, options->MaxRays);
  94     D = DomainIntersection(P, CA, options->MaxRays);
  95     Domain_Free(CA);
  96
  97     if (emptyQ(D)) {
  98         Domain_Free(D);
  99         return evalue_zero();
 100     }
 101
 102     value_init(e.d);
 103     e.x.p = new_enode(partition, 2, P->Dimension);
 104     EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
 105     evalue_set_si(&e.x.p->arr[1], 1, 1);
 106     sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
 107     free_evalue_refs(&e);
 108     return sum;
 109 }