doc: add experimental comparison between old and new Laurent expansion method
[barvinok.git] / summate.c
blobbf15e1a66ef55f334a219db9f2e622dcb747ff44
1 #include <barvinok/options.h>
2 #include <barvinok/util.h>
3 #include "bernoulli.h"
4 #include "euler.h"
5 #include "laurent.h"
6 #include "summate.h"
7 #include "section_array.h"
8 #include "remove_equalities.h"
10 extern evalue *evalue_outer_floor(evalue *e);
11 extern int evalue_replace_floor(evalue *e, const evalue *floor, int var);
12 extern void evalue_drop_floor(evalue *e, const evalue *floor);
14 #define ALLOC(type) (type*)malloc(sizeof(type))
15 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
17 /* Apply the variable transformation specified by T and CP on
18 * the polynomial e. T expresses the old variables in terms
19 * of the new variables (and optionally also the new parameters),
20 * while CP expresses the old parameters in terms of the new
21 * parameters.
23 static void transform_polynomial(evalue *E, Matrix *T, Matrix *CP,
24 unsigned nvar, unsigned nparam,
25 unsigned new_nvar, unsigned new_nparam)
27 int j;
28 evalue **subs;
30 subs = ALLOCN(evalue *, nvar+nparam);
32 for (j = 0; j < nvar; ++j) {
33 if (T)
34 subs[j] = affine2evalue(T->p[j], T->p[T->NbRows-1][T->NbColumns-1],
35 T->NbColumns-1);
36 else
37 subs[j] = evalue_var(j);
39 for (j = 0; j < nparam; ++j) {
40 if (CP)
41 subs[nvar+j] = affine2evalue(CP->p[j], CP->p[nparam][new_nparam],
42 new_nparam);
43 else
44 subs[nvar+j] = evalue_var(j);
45 evalue_shift_variables(subs[nvar+j], 0, new_nvar);
48 evalue_substitute(E, subs);
49 reduce_evalue(E);
51 for (j = 0; j < nvar+nparam; ++j)
52 evalue_free(subs[j]);
53 free(subs);
56 static evalue *sum_over_polytope_with_equalities(Polyhedron *P, evalue *E,
57 unsigned nvar,
58 struct evalue_section_array *sections,
59 struct barvinok_options *options)
61 unsigned dim = P->Dimension;
62 unsigned new_dim, new_nparam;
63 Matrix *T = NULL, *CP = NULL;
64 evalue *sum;
66 if (emptyQ(P))
67 return evalue_zero();
69 assert(P->NbEq > 0);
71 remove_all_equalities(&P, NULL, &CP, &T, dim-nvar, options->MaxRays);
73 if (emptyQ(P)) {
74 Polyhedron_Free(P);
75 return evalue_zero();
78 new_nparam = CP ? CP->NbColumns-1 : dim - nvar;
79 new_dim = T ? T->NbColumns-1 : nvar + new_nparam;
81 /* We can avoid these substitutions if E is a constant */
82 E = evalue_dup(E);
83 transform_polynomial(E, T, CP, nvar, dim-nvar,
84 new_dim-new_nparam, new_nparam);
86 if (new_dim-new_nparam > 0) {
87 sum = barvinok_sum_over_polytope(P, E, new_dim-new_nparam,
88 sections, options);
89 evalue_free(E);
90 Polyhedron_Free(P);
91 } else {
92 sum = ALLOC(evalue);
93 value_init(sum->d);
94 sum->x.p = new_enode(partition, 2, new_dim);
95 EVALUE_SET_DOMAIN(sum->x.p->arr[0], P);
96 value_clear(sum->x.p->arr[1].d);
97 sum->x.p->arr[1] = *E;
98 free(E);
101 if (CP) {
102 evalue_backsubstitute(sum, CP, options->MaxRays);
103 Matrix_Free(CP);
106 if (T)
107 Matrix_Free(T);
109 return sum;
112 static evalue *sum_base(Polyhedron *P, evalue *E, unsigned nvar,
113 struct barvinok_options *options)
115 if (options->summation == BV_SUM_EULER)
116 return euler_summate(P, E, nvar, options);
117 else if (options->summation == BV_SUM_LAURENT)
118 return laurent_summate(P, E, nvar, options);
119 else if (options->summation == BV_SUM_LAURENT_OLD)
120 return laurent_summate_old(P, E, nvar, options);
121 assert(0);
124 /* Count the number of non-zero terms in e when viewed as a polynomial
125 * in only the first nvar variables. "count" is the number counted
126 * so far.
128 static int evalue_count_terms(const evalue *e, unsigned nvar, int count)
130 int i;
132 if (EVALUE_IS_ZERO(*e))
133 return count;
135 if (value_zero_p(e->d))
136 assert(e->x.p->type == polynomial);
137 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1)
138 return count+1;
140 for (i = 0; i < e->x.p->size; ++i)
141 count = evalue_count_terms(&e->x.p->arr[i], nvar, count);
143 return count;
146 /* Create placeholder structure for unzipping.
147 * A "polynomial" is created with size terms in variable pos,
148 * with each term having itself as coefficient.
150 static evalue *create_placeholder(int size, int pos)
152 int i, j;
153 evalue *E = ALLOC(evalue);
154 value_init(E->d);
155 E->x.p = new_enode(polynomial, size, pos+1);
156 for (i = 0; i < size; ++i) {
157 E->x.p->arr[i].x.p = new_enode(polynomial, i+1, pos+1);
158 for (j = 0; j < i; ++j)
159 evalue_set_si(&E->x.p->arr[i].x.p->arr[j], 0, 1);
160 evalue_set_si(&E->x.p->arr[i].x.p->arr[i], 1, 1);
162 return E;
165 /* Interchange each non-zero term in e (when viewed as a polynomial
166 * in only the first nvar variables) with a placeholder in ph (created
167 * by create_placeholder), resulting in two polynomials in the
168 * placeholder variable such that for each non-zero term in e
169 * there is a power of the placeholder variable such that the factors
170 * in the first nvar variables form the coefficient of that power in
171 * the first polynomial (e) and the factors in the remaining variables
172 * form the coefficient of that power in the second polynomial (ph).
174 static int evalue_unzip_terms(evalue *e, evalue *ph, unsigned nvar, int count)
176 int i;
178 if (EVALUE_IS_ZERO(*e))
179 return count;
181 if (value_zero_p(e->d))
182 assert(e->x.p->type == polynomial);
183 if (value_notzero_p(e->d) || e->x.p->pos >= nvar+1) {
184 evalue t = *e;
185 *e = ph->x.p->arr[count];
186 ph->x.p->arr[count] = t;
187 return count+1;
190 for (i = 0; i < e->x.p->size; ++i)
191 count = evalue_unzip_terms(&e->x.p->arr[i], ph, nvar, count);
193 return count;
196 /* Remove n variables at pos (0-based) from the polyhedron P.
197 * Each of these variables is assumed to be completely free,
198 * i.e., there is a line in the polyhedron corresponding to
199 * each of these variables.
201 static Polyhedron *Polyhedron_Remove_Columns(Polyhedron *P, unsigned pos,
202 unsigned n)
204 int i, j;
205 unsigned NbConstraints = 0;
206 unsigned NbRays = 0;
207 Polyhedron *Q;
209 if (n == 0)
210 return P;
212 assert(pos <= P->Dimension);
214 if (POL_HAS(P, POL_INEQUALITIES))
215 NbConstraints = P->NbConstraints;
216 if (POL_HAS(P, POL_POINTS))
217 NbRays = P->NbRays - n;
219 Q = Polyhedron_Alloc(P->Dimension - n, NbConstraints, NbRays);
220 if (POL_HAS(P, POL_INEQUALITIES)) {
221 Q->NbEq = P->NbEq;
222 for (i = 0; i < P->NbConstraints; ++i) {
223 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
224 Vector_Copy(P->Constraint[i]+1+pos+n, Q->Constraint[i]+1+pos,
225 Q->Dimension-pos+1);
228 if (POL_HAS(P, POL_POINTS)) {
229 Q->NbBid = P->NbBid - n;
230 for (i = 0; i < n; ++i)
231 value_set_si(Q->Ray[i][1+pos+i], 1);
232 for (i = 0, j = 0; i < P->NbRays; ++i) {
233 int line = First_Non_Zero(P->Ray[i], 1+P->Dimension+1);
234 assert(line != -1);
235 if (line-1 >= pos && line-1 < pos+n) {
236 ++j;
237 assert(j <= n);
238 continue;
240 assert(i-j < Q->NbRays);
241 Vector_Copy(P->Ray[i], Q->Ray[i-j], 1+pos);
242 Vector_Copy(P->Ray[i]+1+pos+n, Q->Ray[i-j]+1+pos,
243 Q->Dimension-pos+1);
246 POL_SET(Q, POL_VALID);
247 if (POL_HAS(P, POL_INEQUALITIES))
248 POL_SET(Q, POL_INEQUALITIES);
249 if (POL_HAS(P, POL_POINTS))
250 POL_SET(Q, POL_POINTS);
251 if (POL_HAS(P, POL_VERTICES))
252 POL_SET(Q, POL_VERTICES);
253 return Q;
256 /* Remove n variables at pos (0-based) from the union of polyhedra P.
257 * Each of these variables is assumed to be completely free,
258 * i.e., there is a line in the polyhedron corresponding to
259 * each of these variables.
261 static Polyhedron *Domain_Remove_Columns(Polyhedron *P, unsigned pos,
262 unsigned n)
264 Polyhedron *R;
265 Polyhedron **next = &R;
267 for (; P; P = P->next) {
268 *next = Polyhedron_Remove_Columns(P, pos, n);
269 next = &(*next)->next;
271 return R;
274 /* Drop n parameters starting at first from partition evalue e */
275 static void drop_parameters(evalue *e, int first, int n)
277 int i;
279 if (EVALUE_IS_ZERO(*e))
280 return;
282 assert(value_zero_p(e->d) && e->x.p->type == partition);
283 for (i = 0; i < e->x.p->size/2; ++i) {
284 Polyhedron *P = EVALUE_DOMAIN(e->x.p->arr[2*i]);
285 Polyhedron *Q = Domain_Remove_Columns(P, first, n);
286 EVALUE_SET_DOMAIN(e->x.p->arr[2*i], Q);
287 Domain_Free(P);
288 evalue_shift_variables(&e->x.p->arr[2*i+1], first, -n);
290 e->x.p->pos -= n;
293 static void extract_term_into(const evalue *src, int var, int exp, evalue *dst)
295 int i;
297 if (value_notzero_p(src->d) ||
298 src->x.p->type != polynomial ||
299 src->x.p->pos > var+1) {
300 if (exp == 0)
301 evalue_copy(dst, src);
302 else
303 evalue_set_si(dst, 0, 1);
304 return;
307 if (src->x.p->pos == var+1) {
308 if (src->x.p->size > exp)
309 evalue_copy(dst, &src->x.p->arr[exp]);
310 else
311 evalue_set_si(dst, 0, 1);
312 return;
315 dst->x.p = new_enode(polynomial, src->x.p->size, src->x.p->pos);
316 for (i = 0; i < src->x.p->size; ++i)
317 extract_term_into(&src->x.p->arr[i], var, exp,
318 &dst->x.p->arr[i]);
321 /* Extract the coefficient of var^exp.
323 static evalue *extract_term(const evalue *e, int var, int exp)
325 int i;
326 evalue *res;
328 if (EVALUE_IS_ZERO(*e))
329 return evalue_zero();
331 assert(value_zero_p(e->d) && e->x.p->type == partition);
332 res = ALLOC(evalue);
333 value_init(res->d);
334 res->x.p = new_enode(partition, e->x.p->size, e->x.p->pos);
335 for (i = 0; i < e->x.p->size/2; ++i) {
336 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
337 Domain_Copy(EVALUE_DOMAIN(e->x.p->arr[2*i])));
338 extract_term_into(&e->x.p->arr[2*i+1], var, exp,
339 &res->x.p->arr[2*i+1]);
340 reduce_evalue(&res->x.p->arr[2*i+1]);
342 return res;
345 /* Insert n free variables at pos (0-based) in the polyhedron P.
347 static Polyhedron *Polyhedron_Insert_Columns(Polyhedron *P, unsigned pos,
348 unsigned n)
350 int i;
351 unsigned NbConstraints = 0;
352 unsigned NbRays = 0;
353 Polyhedron *Q;
355 if (!P)
356 return P;
357 if (n == 0)
358 return P;
360 assert(pos <= P->Dimension);
362 if (POL_HAS(P, POL_INEQUALITIES))
363 NbConstraints = P->NbConstraints;
364 if (POL_HAS(P, POL_POINTS))
365 NbRays = P->NbRays + n;
367 Q = Polyhedron_Alloc(P->Dimension+n, NbConstraints, NbRays);
368 if (POL_HAS(P, POL_INEQUALITIES)) {
369 Q->NbEq = P->NbEq;
370 for (i = 0; i < P->NbConstraints; ++i) {
371 Vector_Copy(P->Constraint[i], Q->Constraint[i], 1+pos);
372 Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[i]+1+pos+n,
373 P->Dimension-pos+1);
376 if (POL_HAS(P, POL_POINTS)) {
377 Q->NbBid = P->NbBid + n;
378 for (i = 0; i < n; ++i)
379 value_set_si(Q->Ray[i][1+pos+i], 1);
380 for (i = 0; i < P->NbRays; ++i) {
381 Vector_Copy(P->Constraint[i], Q->Constraint[n+i], 1+pos);
382 Vector_Copy(P->Constraint[i]+1+pos, Q->Constraint[n+i]+1+pos+n,
383 P->Dimension-pos+1);
386 POL_SET(Q, POL_VALID);
387 if (POL_HAS(P, POL_INEQUALITIES))
388 POL_SET(Q, POL_INEQUALITIES);
389 if (POL_HAS(P, POL_POINTS))
390 POL_SET(Q, POL_POINTS);
391 if (POL_HAS(P, POL_VERTICES))
392 POL_SET(Q, POL_VERTICES);
393 return Q;
396 /* Perform summation of e over a list of 1 or more factors F, with context C.
397 * nvar is the total number of variables in the remaining factors.
398 * extra is the number of placeholder parameters introduced in e,
399 * but not (yet) in F or C.
401 * If there is only one factor left, F is intersected with the
402 * context C, the placeholder variables are added, and then
403 * e is summed over the resulting parametric polytope.
405 * If there is more than one factor left, we create two polynomials
406 * in a new placeholder variable (which is placed after the regular
407 * parameters, but before any previously introduced placeholder
408 * variables) that has the factors of the variables in the first
409 * factor of F and the factor of the remaining variables of
410 * each term as its coefficients.
411 * These two polynomials are then summed over their domains
412 * and afterwards the results are combined and the placeholder
413 * variable is removed again.
415 static evalue *sum_factors(Polyhedron *F, Polyhedron *C, evalue *e,
416 unsigned nvar, unsigned extra,
417 struct barvinok_options *options)
419 Polyhedron *P = F;
420 unsigned nparam = C->Dimension;
421 unsigned F_var = F->Dimension - C->Dimension;
422 int i, n;
423 evalue *s1, *s2, *s;
424 evalue *ph;
426 if (!F->next) {
427 Polyhedron *CA = align_context(C, nvar+nparam, options->MaxRays);
428 Polyhedron *P = DomainIntersection(F, CA, options->MaxRays);
429 Polyhedron *Q = Polyhedron_Insert_Columns(P, nvar+nparam, extra);
430 Polyhedron_Free(CA);
431 Polyhedron_Free(F);
432 Polyhedron_Free(P);
433 evalue *sum = sum_base(Q, e, nvar, options);
434 Polyhedron_Free(Q);
435 return sum;
438 n = evalue_count_terms(e, F_var, 0);
439 ph = create_placeholder(n, nvar+nparam);
440 evalue_shift_variables(e, nvar+nparam, 1);
441 evalue_unzip_terms(e, ph, F_var, 0);
442 evalue_shift_variables(e, nvar, -(nvar-F_var));
443 evalue_reorder_terms(ph);
444 evalue_shift_variables(ph, 0, -F_var);
446 s2 = sum_factors(F->next, C, ph, nvar-F_var, extra+1, options);
447 evalue_free(ph);
448 F->next = NULL;
449 s1 = sum_factors(F, C, e, F_var, extra+1, options);
451 if (n == 1) {
452 /* remove placeholder "polynomial" */
453 reduce_evalue(s1);
454 emul(s1, s2);
455 evalue_free(s1);
456 drop_parameters(s2, nparam, 1);
457 return s2;
460 s = evalue_zero();
461 for (i = 0; i < n; ++i) {
462 evalue *t1, *t2;
463 t1 = extract_term(s1, nparam, i);
464 t2 = extract_term(s2, nparam, i);
465 emul(t1, t2);
466 eadd(t2, s);
467 evalue_free(t1);
468 evalue_free(t2);
470 evalue_free(s1);
471 evalue_free(s2);
473 drop_parameters(s, nparam, 1);
474 return s;
477 /* Perform summation over a product of factors F, obtained using
478 * variable transformation T from the original problem specification.
480 * We first perform the corresponding transformation on the polynomial E,
481 * compute the common context over all factors and then perform
482 * the actual summation over the factors.
484 static evalue *sum_product(Polyhedron *F, evalue *E, Matrix *T, unsigned nparam,
485 struct barvinok_options *options)
487 int i;
488 Matrix *T2;
489 unsigned nvar = T->NbRows;
490 Polyhedron *C;
491 evalue *sum;
493 assert(nvar == T->NbColumns);
494 T2 = Matrix_Alloc(nvar+1, nvar+1);
495 for (i = 0; i < nvar; ++i)
496 Vector_Copy(T->p[i], T2->p[i], nvar);
497 value_set_si(T2->p[nvar][nvar], 1);
499 transform_polynomial(E, T2, NULL, nvar, nparam, nvar, nparam);
501 C = Factor_Context(F, nparam, options->MaxRays);
502 if (F->Dimension == nparam) {
503 Polyhedron *T = F;
504 F = F->next;
505 T->next = NULL;
506 Polyhedron_Free(T);
508 sum = sum_factors(F, C, E, nvar, 0, options);
510 Polyhedron_Free(C);
511 Matrix_Free(T2);
512 Matrix_Free(T);
513 return sum;
516 /* Add two constraints corresponding to floor = floor(e/d),
518 * e - d t >= 0
519 * -e + d t + d-1 >= 0
521 * e is assumed to be an affine expression.
523 Polyhedron *add_floor_var(Polyhedron *P, unsigned nvar, const evalue *floor,
524 struct barvinok_options *options)
526 int i;
527 unsigned dim = P->Dimension+1;
528 Matrix *M = Matrix_Alloc(P->NbConstraints+2, 2+dim);
529 Polyhedron *CP;
530 Value *d = &M->p[0][1+nvar];
531 evalue_extract_affine(floor, M->p[0]+1, M->p[0]+1+dim, d);
532 value_oppose(*d, *d);
533 value_set_si(M->p[0][0], 1);
534 value_set_si(M->p[1][0], 1);
535 Vector_Oppose(M->p[0]+1, M->p[1]+1, M->NbColumns-1);
536 value_subtract(M->p[1][1+dim], M->p[1][1+dim], *d);
537 value_decrement(M->p[1][1+dim], M->p[1][1+dim]);
539 for (i = 0; i < P->NbConstraints; ++i) {
540 Vector_Copy(P->Constraint[i], M->p[i+2], 1+nvar);
541 Vector_Copy(P->Constraint[i]+1+nvar, M->p[i+2]+1+nvar+1, dim-nvar-1+1);
544 CP = Constraints2Polyhedron(M, options->MaxRays);
545 Matrix_Free(M);
546 return CP;
549 static evalue *evalue_add(evalue *a, evalue *b)
551 if (!a)
552 return b;
553 if (!b)
554 return a;
555 eadd(a, b);
556 evalue_free(a);
557 return b;
560 /* Compute sum of a step-polynomial over a polytope by grouping
561 * terms containing the same floor-expressions and introducing
562 * new variables for each such expression.
563 * In particular, while there is any floor-expression left,
564 * the step-polynomial is split into a polynomial containing
565 * the expression, which is then converted to a new variable,
566 * and a polynomial not containing the expression.
568 static evalue *sum_step_polynomial(Polyhedron *P, evalue *E, unsigned nvar,
569 struct barvinok_options *options)
571 evalue *floor;
572 evalue *cur = E;
573 evalue *sum = NULL;
574 evalue *t;
576 while ((floor = evalue_outer_floor(cur))) {
577 Polyhedron *CP;
578 evalue *converted;
579 evalue *converted_floor;
581 /* Ignore floors that do not depend on variables. */
582 if (value_notzero_p(floor->d) || floor->x.p->pos >= nvar+1)
583 break;
585 converted = evalue_dup(cur);
586 converted_floor = evalue_dup(floor);
587 evalue_shift_variables(converted, nvar, 1);
588 evalue_shift_variables(converted_floor, nvar, 1);
589 evalue_replace_floor(converted, converted_floor, nvar);
590 CP = add_floor_var(P, nvar, converted_floor, options);
591 evalue_free(converted_floor);
592 t = sum_step_polynomial(CP, converted, nvar+1, options);
593 evalue_free(converted);
594 Polyhedron_Free(CP);
595 sum = evalue_add(t, sum);
597 if (cur == E)
598 cur = evalue_dup(cur);
599 evalue_drop_floor(cur, floor);
600 evalue_free(floor);
602 if (floor) {
603 evalue_floor2frac(cur);
604 evalue_free(floor);
607 if (EVALUE_IS_ZERO(*cur))
608 t = NULL;
609 else {
610 Matrix *T;
611 unsigned nparam = P->Dimension - nvar;
612 Polyhedron *F = Polyhedron_Factor(P, nparam, &T, options->MaxRays);
613 if (!F)
614 t = sum_base(P, cur, nvar, options);
615 else {
616 if (cur == E)
617 cur = evalue_dup(cur);
618 t = sum_product(F, cur, T, nparam, options);
622 if (E != cur)
623 evalue_free(cur);
625 return evalue_add(t, sum);
628 evalue *barvinok_sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
629 struct evalue_section_array *sections,
630 struct barvinok_options *options)
632 if (P->NbEq)
633 return sum_over_polytope_with_equalities(P, E, nvar, sections, options);
635 if (options->summation == BV_SUM_BERNOULLI)
636 return bernoulli_summate(P, E, nvar, sections, options);
637 else if (options->summation == BV_SUM_BOX)
638 return box_summate(P, E, nvar, options->MaxRays);
640 evalue_frac2floor2(E, 0);
642 return sum_step_polynomial(P, E, nvar, options);
645 evalue *barvinok_summate(evalue *e, int nvar, struct barvinok_options *options)
647 int i;
648 struct evalue_section_array sections;
649 evalue *sum;
651 assert(nvar >= 0);
652 if (nvar == 0 || EVALUE_IS_ZERO(*e))
653 return evalue_dup(e);
655 assert(value_zero_p(e->d));
656 assert(e->x.p->type == partition);
658 evalue_section_array_init(&sections);
659 sum = evalue_zero();
661 for (i = 0; i < e->x.p->size/2; ++i) {
662 Polyhedron *D;
663 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
664 Polyhedron *next = D->next;
665 evalue *tmp;
666 D->next = NULL;
668 tmp = barvinok_sum_over_polytope(D, &e->x.p->arr[2*i+1], nvar,
669 &sections, options);
670 assert(tmp);
671 eadd(tmp, sum);
672 evalue_free(tmp);
674 D->next = next;
678 free(sections.s);
680 reduce_evalue(sum);
681 return sum;
684 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays)
686 evalue *sum;
687 struct barvinok_options *options = barvinok_options_new_with_defaults();
688 options->MaxRays = MaxRays;
689 sum = barvinok_summate(E, nvar, options);
690 barvinok_options_free(options);
691 return sum;
694 evalue *esum(evalue *e, int nvar)
696 evalue *sum;
697 struct barvinok_options *options = barvinok_options_new_with_defaults();
698 sum = barvinok_summate(e, nvar, options);
699 barvinok_options_free(options);
700 return sum;
703 /* Turn unweighted counting problem into "weighted" counting problem
704 * with weight equal to 1 and call barvinok_summate on this weighted problem.
706 evalue *barvinok_summate_unweighted(Polyhedron *P, Polyhedron *C,
707 struct barvinok_options *options)
709 Polyhedron *CA, *D;
710 evalue e;
711 evalue *sum;
713 if (emptyQ(P) || emptyQ(C))
714 return evalue_zero();
716 CA = align_context(C, P->Dimension, options->MaxRays);
717 D = DomainIntersection(P, CA, options->MaxRays);
718 Domain_Free(CA);
720 if (emptyQ(D)) {
721 Domain_Free(D);
722 return evalue_zero();
725 value_init(e.d);
726 e.x.p = new_enode(partition, 2, P->Dimension);
727 EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
728 evalue_set_si(&e.x.p->arr[1], 1, 1);
729 sum = barvinok_summate(&e, P->Dimension - C->Dimension, options);
730 free_evalue_refs(&e);
731 return sum;