iscc: test isl_stream for eof rather than the underlying FILE
[barvinok.git] / euler.cc
blobcb982f1516afb1ad4c34a04049c8aeca2dffc8cc
1 /*
2 * Sum polynomial over integer points in polytope using local
3 * Euler-Maclaurin formula by Berline and Vergne.
4 */
6 #include <assert.h>
7 #include <barvinok/options.h>
8 #include <barvinok/util.h>
9 #include "bernoulli.h"
10 #include "binomial.h"
11 #include "conversion.h"
12 #include "decomposer.h"
13 #include "euler.h"
14 #include "lattice_point.h"
15 #include "param_util.h"
16 #include "power.h"
17 #include "reduce_domain.h"
19 /* Compute total degree in first nvar variables */
20 static unsigned total_degree(const evalue *e, unsigned nvar)
22 int i;
23 unsigned max_degree;
25 if (value_notzero_p(e->d))
26 return 0;
27 assert(e->x.p->type == polynomial);
28 if (e->x.p->pos-1 >= nvar)
29 return 0;
31 max_degree = total_degree(&e->x.p->arr[0], nvar);
32 for (i = 1; i < e->x.p->size; ++i) {
33 unsigned degree = i + total_degree(&e->x.p->arr[i], nvar);
34 if (degree > max_degree)
35 max_degree = degree;
37 return max_degree;
41 * Computes the coefficients of
43 * \mu(-t + R_+)(\xi) = \sum_{n=0)^\infty -b(n+1, t)/(n+1)! \xi^n
45 * where b(n, t) is the Bernoulli polynomial of degree n in t
46 * and t(p) is an expression (a fractional part) of the parameters p
47 * such that 0 <= t(p) < 1 for all values of the parameters.
48 * The coefficients are computed on demand up to (and including)
49 * the maximal degree max_degree.
51 struct mu_1d {
52 unsigned max_degree;
53 evalue **coefficients;
54 const evalue *t;
56 mu_1d(unsigned max_degree, evalue *t) : max_degree(max_degree), t(t) {
57 coefficients = new evalue *[max_degree+1];
58 for (int i = 0; i < max_degree+1; ++i)
59 coefficients[i] = NULL;
61 ~mu_1d() {
62 for (int i = 0; i < max_degree+1; ++i)
63 if (coefficients[i])
64 evalue_free(coefficients[i]);
65 delete [] coefficients;
67 void compute_coefficient(unsigned n);
68 const evalue *coefficient(unsigned n) {
69 if (!coefficients[n])
70 compute_coefficient(n);
71 return coefficients[n];
75 void mu_1d::compute_coefficient(unsigned n)
77 struct poly_list *bernoulli = bernoulli_compute(n+1);
78 evalue *c = evalue_polynomial(bernoulli->poly[n+1], t);
79 evalue_negate(c);
80 evalue_div(c, *factorial(n+1));
82 coefficients[n] = c;
86 * Computes the coefficients of
88 * \mu(a)(y) = \sum_{n_1} \sum_{n_2} c_{n_1,n_2} y^{n_1} y^{n_2}
90 * with c_{n1,n2} given by
92 * b(n1+1,t1)/(n1+1)! b(n2+1,t2)/(n2+1)!
93 * - b(n1+n2+2,t2)/(n1+n2+2)! (-c1)^{n1+1}
94 * - b(n1+n2+2,t1)/(n1+n2+2)! (-c2)^{n2+1}
96 * where b(n, t) is the Bernoulli polynomial of degree n in t,
97 * t1(p) and t2(p) are expressions (fractional parts) of the parameters p
98 * such that 0 <= t1(p), t2(p) < 1 for all values of the parameters
99 * and c1 = cn/c1d and c2 = cn/c2d.
100 * The coefficients are computed on demand up to (and including)
101 * the maximal degree (n1,n2) = (max_degree,max_degree).
103 * bernoulli_t[i][j] contains b(j+1, t_i)/(j+1)!
105 struct mu_2d {
106 unsigned max_degree;
107 evalue ***coefficients;
108 /* bernoulli_t[i][n] stores b(n+1, t_i)/(n+1)! */
109 evalue **bernoulli_t[2];
110 /* stores the powers of -cn */
111 struct power *power_cn;
112 struct power *power_c1d;
113 struct power *power_c2d;
114 const evalue *t[2];
116 mu_2d(unsigned max_degree, evalue *t1, evalue *t2,
117 Value cn, Value c1d, Value c2d) : max_degree(max_degree) {
118 t[0] = t1;
119 t[1] = t2;
120 coefficients = new evalue **[max_degree+1];
121 for (int i = 0; i < max_degree+1; ++i) {
122 coefficients[i] = new evalue *[max_degree+1];
123 for (int j = 0; j < max_degree+1; ++j)
124 coefficients[i][j] = NULL;
126 for (int i = 0; i < 2; ++i) {
127 bernoulli_t[i] = new evalue *[max_degree+2];
128 for (int j = 0; j < max_degree+2; ++j)
129 bernoulli_t[i][j] = NULL;
131 value_oppose(cn, cn);
132 power_cn = new struct power(cn, max_degree+1);
133 value_oppose(cn, cn);
134 power_c1d = new struct power(c1d, max_degree+1);
135 power_c2d = new struct power(c2d, max_degree+1);
137 ~mu_2d() {
138 for (int i = 0; i < max_degree+1; ++i) {
139 for (int j = 0; j < max_degree+1; ++j)
140 if (coefficients[i][j])
141 evalue_free(coefficients[i][j]);
142 delete [] coefficients[i];
144 delete [] coefficients;
145 for (int i = 0; i < 2; ++i)
146 for (int j = 0; j < max_degree+2; ++j)
147 if (bernoulli_t[i][j])
148 evalue_free(bernoulli_t[i][j]);
149 for (int i = 0; i < 2; ++i)
150 delete [] bernoulli_t[i];
151 delete power_cn;
152 delete power_c1d;
153 delete power_c2d;
155 const evalue *get_bernoulli(int n, int i);
157 void compute_coefficient(unsigned n1, unsigned n2);
158 const evalue *coefficient(unsigned n1, unsigned n2) {
159 if (!coefficients[n1][n2])
160 compute_coefficient(n1, n2);
161 return coefficients[n1][n2];
166 * Returns b(n, t_i)/n!
168 const evalue *mu_2d::get_bernoulli(int n, int i)
170 if (!bernoulli_t[i][n-1]) {
171 struct poly_list *bernoulli = bernoulli_compute(n);
172 bernoulli_t[i][n-1] = evalue_polynomial(bernoulli->poly[n], t[i]);
173 evalue_div(bernoulli_t[i][n-1], *factorial(n));
175 return bernoulli_t[i][n-1];
178 void mu_2d::compute_coefficient(unsigned n1, unsigned n2)
180 evalue *c = evalue_dup(get_bernoulli(n1+1, 0));
181 emul(get_bernoulli(n2+1, 1), c);
183 if (value_notzero_p(*(*power_cn)[1])) {
184 evalue *t;
185 Value neg_power;
187 value_init(neg_power);
189 t = evalue_dup(get_bernoulli(n1+n2+2, 1));
190 value_multiply(neg_power,
191 *(*power_cn)[n1+1], *binomial(n1+n2+1, n1+1));
192 value_oppose(neg_power, neg_power);
193 evalue_mul_div(t, neg_power, *(*power_c1d)[n1+1]);
194 eadd(t, c);
195 evalue_free(t);
197 t = evalue_dup(get_bernoulli(n1+n2+2, 0));
198 value_multiply(neg_power,
199 *(*power_cn)[n2+1], *binomial(n1+n2+1, n2+1));
200 value_oppose(neg_power, neg_power);
201 evalue_mul_div(t, neg_power, *(*power_c2d)[n2+1]);
202 eadd(t, c);
203 evalue_free(t);
205 value_clear(neg_power);
208 coefficients[n1][n2] = c;
211 /* Later: avoid recomputation of mu of faces that appear in
212 * more than one domain.
214 struct summator_2d : public signed_cone_consumer, public vertex_decomposer {
215 const evalue *polynomial;
216 Value *inner;
217 unsigned nparam;
219 /* substitutions to use when result is 0-dimensional (only parameters) */
220 evalue **subs_0d;
221 /* substitutions to use when result is 1-dimensional */
222 evalue **subs_1d;
223 evalue *sum;
225 summator_2d(evalue *e, Param_Polyhedron *PP, Value *inner,
226 unsigned nparam) :
227 polynomial(e), vertex_decomposer(PP, *this),
228 inner(inner), nparam(nparam) {
229 sum = evalue_zero();
231 subs_0d = new evalue *[2+nparam];
232 subs_1d = new evalue *[2+nparam];
233 subs_0d[0] = NULL;
234 subs_0d[1] = NULL;
235 subs_1d[0] = NULL;
236 subs_1d[1] = NULL;
237 for (int i = 0; i < nparam; ++i) {
238 subs_0d[2+i] = evalue_var(i);
239 subs_1d[2+i] = evalue_var(1+i);
242 ~summator_2d() {
243 for (int i = 0; i < nparam; ++i) {
244 evalue_free(subs_0d[2+i]);
245 evalue_free(subs_1d[2+i]);
247 delete [] subs_0d;
248 delete [] subs_1d;
250 evalue *summate_over_pdomain(Param_Polyhedron *PP, unsigned *facets,
251 Param_Domain *PD,
252 struct barvinok_options *options);
253 void handle_facet(Param_Polyhedron *PP, Param_Domain *FD, Value *normal);
254 void integrate(Param_Polyhedron *PP, unsigned *facets, Param_Domain *PD);
255 virtual void handle(const signed_cone& sc, barvinok_options *options);
258 /* Replaces poly by its derivative along variable var */
259 static void evalue_derive(evalue *poly, int var)
261 if (value_notzero_p(poly->d)) {
262 value_set_si(poly->x.n, 0);
263 value_set_si(poly->d, 1);
264 return;
266 assert(poly->x.p->type == polynomial);
267 if (poly->x.p->pos-1 > var) {
268 free_evalue_refs(poly);
269 value_init(poly->d);
270 evalue_set_si(poly, 0, 1);
271 return;
274 if (poly->x.p->pos-1 < var) {
275 for (int i = 0; i < poly->x.p->size; ++i)
276 evalue_derive(&poly->x.p->arr[i], var);
277 reduce_evalue(poly);
278 return;
281 assert(poly->x.p->size >= 1);
282 enode *p = poly->x.p;
283 free_evalue_refs(&p->arr[0]);
284 if (p->size == 1) {
285 free(p);
286 evalue_set_si(poly, 0, 1);
287 return;
289 if (p->size == 2) {
290 value_clear(poly->d);
291 *poly = p->arr[1];
292 free(p);
293 return;
296 p->size--;
297 Value factor;
298 value_init(factor);
299 for (int i = 0; i < p->size; ++i) {
300 value_set_si(factor, i+1);
301 p->arr[i] = p->arr[i+1];
302 evalue_mul(&p->arr[i], factor);
304 value_clear(factor);
307 /* Check whether e is constant with respect to variable var */
308 static int evalue_is_constant(evalue *e, int var)
310 int i;
312 if (value_notzero_p(e->d))
313 return 1;
314 if (e->x.p->type == polynomial && e->x.p->pos-1 == var)
315 return 0;
316 assert(e->x.p->type == polynomial ||
317 e->x.p->type == fractional ||
318 e->x.p->type == flooring);
319 for (i = 0; i < e->x.p->size; ++i)
320 if (!evalue_is_constant(&e->x.p->arr[i], var))
321 return 0;
322 return 1;
325 /* Replaces poly by its anti-derivative with constant 0 along variable var */
326 static void evalue_anti_derive(evalue *poly, int var)
328 enode *p;
330 if (value_zero_p(poly->d) &&
331 poly->x.p->type == polynomial && poly->x.p->pos-1 < var) {
332 for (int i = 0; i < poly->x.p->size; ++i)
333 evalue_anti_derive(&poly->x.p->arr[i], var);
334 reduce_evalue(poly);
335 return;
338 if (evalue_is_constant(poly, var)) {
339 p = new_enode(polynomial, 2, 1+var);
340 evalue_set_si(&p->arr[0], 0, 1);
341 value_clear(p->arr[1].d);
342 p->arr[1] = *poly;
343 value_init(poly->d);
344 poly->x.p = p;
345 return;
347 assert(poly->x.p->type == polynomial);
349 p = new_enode(polynomial, poly->x.p->size+1, 1+var);
350 evalue_set_si(&p->arr[0], 0, 1);
351 for (int i = 0; i < poly->x.p->size; ++i) {
352 value_clear(p->arr[1+i].d);
353 p->arr[1+i] = poly->x.p->arr[i];
354 value_set_si(poly->d, 1+i);
355 evalue_div(&p->arr[1+i], poly->d);
357 free(poly->x.p);
358 poly->x.p = p;
359 value_set_si(poly->d, 0);
362 /* Computes offsets in the basis given by the rays of the cone
363 * to the integer point in the fundamental parallelepiped of
364 * the cone.
365 * The resulting evalues contain only the parameters as variables.
367 evalue **offsets_to_integer_point(Matrix *Rays, Matrix *vertex)
369 unsigned nparam = vertex->NbColumns-2;
370 evalue **t = new evalue *[2];
372 if (value_one_p(vertex->p[0][nparam+1])) {
373 t[0] = evalue_zero();
374 t[1] = evalue_zero();
375 } else {
376 Matrix *R2 = Matrix_Copy(Rays);
377 Matrix_Transposition(R2);
378 Matrix *inv = Matrix_Alloc(Rays->NbColumns, Rays->NbRows);
379 int ok = Matrix_Inverse(R2, inv);
380 assert(ok);
381 Matrix_Free(R2);
383 /* We want the fractional part of the negative relative coordinates */
384 Vector_Oppose(inv->p[0], inv->p[0], inv->NbColumns);
385 Vector_Oppose(inv->p[1], inv->p[1], inv->NbColumns);
387 Matrix *neg_rel = Matrix_Alloc(2, nparam+1);
388 vertex->NbColumns--;
389 Matrix_Product(inv, vertex, neg_rel);
390 Matrix_Free(inv);
391 vertex->NbColumns++;
392 t[0] = fractional_part(neg_rel->p[0], vertex->p[0][nparam+1],
393 nparam, NULL);
394 t[1] = fractional_part(neg_rel->p[1], vertex->p[0][nparam+1],
395 nparam, NULL);
396 Matrix_Free(neg_rel);
399 return t;
403 * Called from decompose_at_vertex.
405 * Handles a cone in the signed decomposition of the supporting
406 * cone of a vertex. The cone is assumed to be unimodular.
408 void summator_2d::handle(const signed_cone& sc, barvinok_options *options)
410 evalue **t;
411 unsigned degree = total_degree(polynomial, 2);
413 subs_0d[0] = affine2evalue(V->Vertex->p[0],
414 V->Vertex->p[0][nparam+1], nparam);
415 subs_0d[1] = affine2evalue(V->Vertex->p[1],
416 V->Vertex->p[1][nparam+1], nparam);
418 assert(sc.det == 1);
420 assert(V->Vertex->NbRows > 0);
421 Param_Vertex_Common_Denominator(V);
423 Matrix *Rays = zz2matrix(sc.rays);
425 t = offsets_to_integer_point(Rays, V->Vertex);
427 Vector *c = Vector_Alloc(3);
428 Inner_Product(Rays->p[0], Rays->p[1], 2, &c->p[0]);
429 Inner_Product(Rays->p[0], Rays->p[0], 2, &c->p[1]);
430 Inner_Product(Rays->p[1], Rays->p[1], 2, &c->p[2]);
432 mu_2d mu(degree, t[0], t[1], c->p[0], c->p[1], c->p[2]);
433 Vector_Free(c);
435 struct power power_r00(Rays->p[0][0], degree);
436 struct power power_r01(Rays->p[0][1], degree);
437 struct power power_r10(Rays->p[1][0], degree);
438 struct power power_r11(Rays->p[1][1], degree);
440 Value factor, tmp1, tmp2;
441 value_init(factor);
442 value_init(tmp1);
443 value_init(tmp2);
444 evalue *res = evalue_zero();
445 evalue *dx1 = evalue_dup(polynomial);
446 for (int i = 0; !EVALUE_IS_ZERO(*dx1); ++i) {
447 evalue *dx2 = evalue_dup(dx1);
448 for (int j = 0; !EVALUE_IS_ZERO(*dx2); ++j) {
449 evalue *cij = evalue_zero();
450 for (int n1 = 0; n1 <= i+j; ++n1) {
451 int n2 = i+j-n1;
452 value_set_si(factor, 0);
453 for (int k = max(0, i-n2); k <= i && k <= n1; ++k) {
454 int l = i-k;
455 value_multiply(tmp1, *power_r00[k], *power_r01[n1-k]);
456 value_multiply(tmp1, tmp1, *power_r10[l]);
457 value_multiply(tmp1, tmp1, *power_r11[n2-l]);
458 value_multiply(tmp1, tmp1, *binomial(n1, k));
459 value_multiply(tmp1, tmp1, *binomial(n2, l));
460 value_addto(factor, factor, tmp1);
462 if (value_zero_p(factor))
463 continue;
465 evalue *c = evalue_dup(mu.coefficient(n1, n2));
466 evalue_mul(c, factor);
467 eadd(c, cij);
468 evalue_free(c);
470 evalue *d = evalue_dup(dx2);
471 evalue_substitute(d, subs_0d);
472 emul(cij, d);
473 evalue_free(cij);
474 eadd(d, res);
475 evalue_free(d);
476 evalue_derive(dx2, 1);
478 evalue_free(dx2);
479 evalue_derive(dx1, 0);
481 evalue_free(dx1);
482 for (int i = 0; i < 2; ++i) {
483 evalue_free(subs_0d[i]);
484 subs_0d[i] = NULL;
485 evalue_free(t[i]);
487 delete [] t;
488 value_clear(factor);
489 value_clear(tmp1);
490 value_clear(tmp2);
491 Matrix_Free(Rays);
493 if (sc.sign < 0)
494 evalue_negate(res);
495 eadd(res, sum);
496 evalue_free(res);
499 evalue *summator_2d::summate_over_pdomain(Param_Polyhedron *PP,
500 unsigned *facets, Param_Domain *PD,
501 struct barvinok_options *options)
503 Param_Vertices *V;
504 int i, ix;
505 unsigned bx;
507 assert(PP->V->Vertex->NbRows == 2);
509 FORALL_PVertex_in_ParamPolyhedron(V, PD, PP) // _i internal counter
510 decompose_at_vertex(V, _i, options);
511 END_FORALL_PVertex_in_ParamPolyhedron;
513 Vector *normal = Vector_Alloc(2);
514 for (i = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
515 Param_Domain *FD;
517 if (!(facets[ix] & bx)) {
518 NEXT(ix, bx);
519 continue;
522 Vector_Copy(PP->Constraints->p[i]+1, normal->p, 2);
523 if (value_zero_p(normal->p[0]) && value_zero_p(normal->p[1]))
524 continue;
526 FD = Param_Polyhedron_Facet(PP, PD, PP->Constraints->p[i]);
527 Vector_Normalize(normal->p, 2);
528 handle_facet(PP, FD, normal->p);
529 Param_Domain_Free(FD);
530 NEXT(ix, bx);
532 Vector_Free(normal);
534 integrate(PP, facets, PD);
536 return sum;
539 void summator_2d::handle_facet(Param_Polyhedron *PP, Param_Domain *FD,
540 Value *normal)
542 Param_Vertices *V;
543 int nbV = 0;
544 Param_Vertices *vertex[2];
545 Value det;
546 unsigned degree = total_degree(polynomial, 2);
548 FORALL_PVertex_in_ParamPolyhedron(V, FD, PP)
549 vertex[nbV++] = V;
550 END_FORALL_PVertex_in_ParamPolyhedron;
552 assert(nbV == 2);
554 /* We can take either vertex[0] or vertex[1];
555 * the result should be the same
557 Param_Vertex_Common_Denominator(vertex[0]);
559 /* The extra variable in front is the coordinate along the facet. */
560 Vector *coef_normal = Vector_Alloc(1 + nparam + 2);
561 Vector_Combine(vertex[0]->Vertex->p[0], vertex[0]->Vertex->p[1],
562 coef_normal->p+1, normal[0], normal[1], nparam+1);
563 value_assign(coef_normal->p[1+nparam+1], vertex[0]->Vertex->p[0][nparam+1]);
564 Vector_Normalize(coef_normal->p, coef_normal->Size);
566 Vector *base = Vector_Alloc(2);
567 value_oppose(base->p[0], normal[1]);
568 value_assign(base->p[1], normal[0]);
570 value_init(det);
571 Inner_Product(normal, normal, 2, &det);
573 Vector *s = Vector_Alloc(1+nparam+2);
574 value_multiply(s->p[1+nparam+1], coef_normal->p[1+nparam+1], det);
575 value_multiply(s->p[0], base->p[0], s->p[1+nparam+1]);
576 Vector_Scale(coef_normal->p+1, s->p+1, normal[0], nparam+1);
577 subs_1d[0] = affine2evalue(s->p, s->p[1+nparam+1], 1+nparam);
578 value_multiply(s->p[0], base->p[1], s->p[1+nparam+1]);
579 Vector_Scale(coef_normal->p+1, s->p+1, normal[1], nparam+1);
580 subs_1d[1] = affine2evalue(s->p, s->p[1+nparam+1], 1+nparam);
581 Vector_Free(s);
583 evalue *t;
584 if (value_one_p(coef_normal->p[coef_normal->Size-1]))
585 t = evalue_zero();
586 else {
587 Vector_Oppose(coef_normal->p+1, coef_normal->p+1, nparam+1);
588 t = fractional_part(coef_normal->p,
589 coef_normal->p[coef_normal->Size-1],
590 1+nparam, NULL);
592 Vector_Free(coef_normal);
594 mu_1d mu(degree, t);
596 struct power power_normal0(normal[0], degree);
597 struct power power_normal1(normal[1], degree);
598 struct power power_det(det, degree);
600 Value(factor);
601 value_init(factor);
602 evalue *res = evalue_zero();
603 evalue *dx1 = evalue_dup(polynomial);
604 for (int i = 0; !EVALUE_IS_ZERO(*dx1); ++i) {
605 evalue *dx2 = evalue_dup(dx1);
606 for (int j = 0; !EVALUE_IS_ZERO(*dx2); ++j) {
607 value_multiply(factor, *power_normal0[i], *power_normal1[j]);
608 if (value_notzero_p(factor)) {
609 value_multiply(factor, factor, *binomial(i+j, i));
611 evalue *c = evalue_dup(mu.coefficient(i+j));
612 evalue_mul_div(c, factor, *power_det[i+j]);
614 evalue *d = evalue_dup(dx2);
615 evalue_substitute(d, subs_1d);
616 emul(c, d);
617 evalue_free(c);
618 eadd(d, res);
619 evalue_free(d);
621 evalue_derive(dx2, 1);
623 evalue_free(dx2);
624 evalue_derive(dx1, 0);
626 evalue_free(dx1);
627 evalue_free(t);
628 for (int i = 0; i < 2; ++i) {
629 evalue_free(subs_1d[i]);
630 subs_1d[i] = NULL;
632 value_clear(factor);
633 evalue_anti_derive(res, 0);
635 Matrix *z = Matrix_Alloc(2, nparam+2);
636 Vector *fixed_z = Vector_Alloc(2);
637 for (int i = 0; i < 2; ++i) {
638 Vector_Combine(vertex[i]->Vertex->p[0], vertex[i]->Vertex->p[1],
639 z->p[i], base->p[0], base->p[1], nparam+1);
640 value_multiply(z->p[i][nparam+1],
641 det, vertex[i]->Vertex->p[0][nparam+1]);
642 Inner_Product(z->p[i], inner, nparam+1, &fixed_z->p[i]);
644 Vector_Free(base);
645 /* Put on a common denominator */
646 value_multiply(fixed_z->p[0], fixed_z->p[0], z->p[1][nparam+1]);
647 value_multiply(fixed_z->p[1], fixed_z->p[1], z->p[0][nparam+1]);
648 /* Make sure z->p[0] is smaller than z->p[1] (for an internal
649 * point of the chamber and hence for all parameter values in
650 * the chamber), to ensure we integrate in the right direction.
652 if (value_lt(fixed_z->p[1], fixed_z->p[0]))
653 Vector_Exchange(z->p[0], z->p[1], nparam+2);
654 Vector_Free(fixed_z);
655 value_clear(det);
657 subs_0d[1] = affine2evalue(z->p[1], z->p[1][nparam+1], nparam);
658 evalue *up = evalue_dup(res);
659 evalue_substitute(up, subs_0d+1);
660 evalue_free(subs_0d[1]);
662 subs_0d[1] = affine2evalue(z->p[0], z->p[0][nparam+1], nparam);
663 evalue_substitute(res, subs_0d+1);
664 evalue_negate(res);
665 eadd(up, res);
666 evalue_free(subs_0d[1]);
667 subs_0d[1] = NULL;
668 evalue_free(up);
670 Matrix_Free(z);
672 eadd(res, sum);
673 evalue_free(res);
676 /* Integrate the polynomial over the whole polygon using
677 * the Green-Stokes theorem.
679 void summator_2d::integrate(Param_Polyhedron *PP, unsigned *facets,
680 Param_Domain *PD)
682 Value tmp;
683 evalue *res = evalue_zero();
684 int i, ix;
685 unsigned bx;
687 evalue *I = evalue_dup(polynomial);
688 evalue_anti_derive(I, 0);
690 value_init(tmp);
691 Vector *normal = Vector_Alloc(2);
692 Vector *dir = Vector_Alloc(2);
693 Matrix *v0v1 = Matrix_Alloc(2, nparam+2);
694 Vector *f_v0v1 = Vector_Alloc(2);
695 Vector *s = Vector_Alloc(1+nparam+2);
696 for (i = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
697 Param_Domain *FD;
698 int nbV = 0;
699 Param_Vertices *vertex[2];
701 if (!(facets[ix] & bx)) {
702 NEXT(ix, bx);
703 continue;
706 Vector_Copy(PP->Constraints->p[i]+1, normal->p, 2);
708 if (value_zero_p(normal->p[0]))
709 continue;
711 Vector_Normalize(normal->p, 2);
712 value_assign(dir->p[0], normal->p[1]);
713 value_oppose(dir->p[1], normal->p[0]);
715 FD = Param_Polyhedron_Facet(PP, PD, PP->Constraints->p[i]);
717 FORALL_PVertex_in_ParamPolyhedron(V, FD, PP)
718 vertex[nbV++] = V;
719 END_FORALL_PVertex_in_ParamPolyhedron;
721 assert(nbV == 2);
723 Param_Vertex_Common_Denominator(vertex[0]);
724 Param_Vertex_Common_Denominator(vertex[1]);
726 value_oppose(tmp, vertex[1]->Vertex->p[0][nparam+1]);
727 for (int i = 0; i < 2; ++i)
728 Vector_Combine(vertex[1]->Vertex->p[i],
729 vertex[0]->Vertex->p[i],
730 v0v1->p[i],
731 vertex[0]->Vertex->p[0][nparam+1], tmp, nparam+1);
732 value_multiply(v0v1->p[0][nparam+1],
733 vertex[0]->Vertex->p[0][nparam+1],
734 vertex[1]->Vertex->p[0][nparam+1]);
735 value_assign(v0v1->p[1][nparam+1], v0v1->p[0][nparam+1]);
737 /* Order vertices to ensure we integrate in the right
738 * direction, i.e., with the polytope "on the left".
740 for (int i = 0; i < 2; ++i)
741 Inner_Product(v0v1->p[i], inner, nparam+1, &f_v0v1->p[i]);
743 Inner_Product(dir->p, f_v0v1->p, 2, &tmp);
744 if (value_neg_p(tmp)) {
745 Param_Vertices *PV = vertex[0];
746 vertex[0] = vertex[1];
747 vertex[1] = PV;
748 for (int i = 0; i < 2; ++i)
749 Vector_Oppose(v0v1->p[i], v0v1->p[i], nparam+1);
751 value_oppose(tmp, normal->p[0]);
752 if (value_neg_p(tmp)) {
753 value_oppose(tmp, tmp);
754 Vector_Oppose(v0v1->p[1], v0v1->p[1], nparam+1);
756 value_multiply(tmp, tmp, v0v1->p[1][nparam+1]);
757 evalue *top = affine2evalue(v0v1->p[1], tmp, nparam);
759 value_multiply(s->p[0], normal->p[1], vertex[0]->Vertex->p[0][nparam+1]);
760 Vector_Copy(vertex[0]->Vertex->p[0], s->p+1, nparam+2);
761 subs_1d[0] = affine2evalue(s->p, s->p[1+nparam+1], 1+nparam);
762 value_multiply(s->p[0], normal->p[0], vertex[0]->Vertex->p[0][nparam+1]);
763 value_oppose(s->p[0], s->p[0]);
764 Vector_Copy(vertex[0]->Vertex->p[1], s->p+1, nparam+2);
765 subs_1d[1] = affine2evalue(s->p, s->p[1+nparam+1], 1+nparam);
767 evalue *d = evalue_dup(I);
768 evalue_substitute(d, subs_1d);
769 evalue_anti_derive(d, 0);
771 evalue_free(subs_1d[0]);
772 evalue_free(subs_1d[1]);
774 subs_0d[1] = top;
775 evalue_substitute(d, subs_0d+1);
776 evalue_mul(d, dir->p[1]);
777 evalue_free(subs_0d[1]);
779 eadd(d, res);
780 evalue_free(d);
782 Param_Domain_Free(FD);
783 NEXT(ix, bx);
785 Vector_Free(s);
786 Vector_Free(f_v0v1);
787 Matrix_Free(v0v1);
788 Vector_Free(normal);
789 Vector_Free(dir);
790 value_clear(tmp);
791 evalue_free(I);
793 eadd(res, sum);
794 evalue_free(res);
797 evalue *summate_over_1d_pdomain(evalue *e,
798 Param_Polyhedron *PP, Param_Domain *PD,
799 Value *inner,
800 struct barvinok_options *options)
802 Param_Vertices *V;
803 int nbV = 0;
804 Param_Vertices *vertex[2];
805 unsigned nparam = PP->V->Vertex->NbColumns-2;
806 evalue **subs_0d = new evalue *[1+nparam];
807 evalue *a[2];
808 evalue *t[2];
809 unsigned degree = total_degree(e, 1);
811 for (int i = 0; i < nparam; ++i)
812 subs_0d[1+i] = evalue_var(i);
814 FORALL_PVertex_in_ParamPolyhedron(V, PD, PP)
815 vertex[nbV++] = V;
816 END_FORALL_PVertex_in_ParamPolyhedron;
817 assert(nbV == 2);
819 Vector *fixed = Vector_Alloc(2);
820 for (int i = 0; i < 2; ++i) {
821 Inner_Product(vertex[i]->Vertex->p[0], inner, nparam+1, &fixed->p[i]);
822 value_multiply(fixed->p[i],
823 fixed->p[i], vertex[1-i]->Vertex->p[0][nparam+1]);
825 if (value_lt(fixed->p[1], fixed->p[0])) {
826 V = vertex[0];
827 vertex[0] = vertex[1];
828 vertex[1] = V;
830 Vector_Free(fixed);
832 Vector *coef = Vector_Alloc(nparam+1);
833 for (int i = 0; i < 2; ++i)
834 a[i] = affine2evalue(vertex[i]->Vertex->p[0],
835 vertex[i]->Vertex->p[0][nparam+1], nparam);
836 if (value_one_p(vertex[0]->Vertex->p[0][nparam+1]))
837 t[0] = evalue_zero();
838 else {
839 Vector_Oppose(vertex[0]->Vertex->p[0], coef->p, nparam+1);
840 t[0] = fractional_part(coef->p, vertex[0]->Vertex->p[0][nparam+1],
841 nparam, NULL);
843 if (value_one_p(vertex[1]->Vertex->p[0][nparam+1]))
844 t[1] = evalue_zero();
845 else {
846 Vector_Copy(vertex[1]->Vertex->p[0], coef->p, nparam+1);
847 t[1] = fractional_part(coef->p, vertex[1]->Vertex->p[0][nparam+1],
848 nparam, NULL);
850 Vector_Free(coef);
852 evalue *I = evalue_dup(e);
853 evalue_anti_derive(I, 0);
854 evalue *up = evalue_dup(I);
855 subs_0d[0] = a[1];
856 evalue_substitute(up, subs_0d);
858 subs_0d[0] = a[0];
859 evalue_substitute(I, subs_0d);
860 evalue_negate(I);
861 eadd(up, I);
862 evalue_free(up);
864 evalue *res = I;
866 mu_1d mu0(degree, t[0]);
867 mu_1d mu1(degree, t[1]);
869 evalue *dx = evalue_dup(e);
870 for (int n = 0; !EVALUE_IS_ZERO(*dx); ++n) {
871 evalue *d;
873 d = evalue_dup(dx);
874 subs_0d[0] = a[0];
875 evalue_substitute(d, subs_0d);
876 emul(mu0.coefficient(n), d);
877 eadd(d, res);
878 evalue_free(d);
880 d = evalue_dup(dx);
881 subs_0d[0] = a[1];
882 evalue_substitute(d, subs_0d);
883 emul(mu1.coefficient(n), d);
884 if (n % 2)
885 evalue_negate(d);
886 eadd(d, res);
887 evalue_free(d);
889 evalue_derive(dx, 0);
891 evalue_free(dx);
893 for (int i = 0; i < nparam; ++i)
894 evalue_free(subs_0d[1+i]);
895 delete [] subs_0d;
897 for (int i = 0; i < 2; ++i) {
898 evalue_free(a[i]);
899 evalue_free(t[i]);
902 return res;
905 #define INT_BITS (sizeof(unsigned) * 8)
907 static unsigned *active_constraints(Param_Polyhedron *PP, Param_Domain *D)
909 int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
910 unsigned *facets = (unsigned *)calloc(len, sizeof(unsigned));
911 Param_Vertices *V;
913 FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
914 if (!V->Facets)
915 Param_Vertex_Set_Facets(PP, V);
916 for (int i = 0; i < len; ++i)
917 facets[i] |= V->Facets[i];
918 END_FORALL_PVertex_in_ParamPolyhedron;
920 return facets;
923 evalue *euler_summate(Polyhedron *P, evalue *e, unsigned nvar,
924 struct barvinok_options *options)
926 Polyhedron *U;
927 Param_Polyhedron *PP;
928 evalue *res;
929 int nd = -1;
930 unsigned MaxRays;
931 struct evalue_section *s;
932 Param_Domain *PD;
934 assert(nvar <= 2);
936 MaxRays = options->MaxRays;
937 POL_UNSET(options->MaxRays, POL_INTEGER);
939 U = Universe_Polyhedron(P->Dimension - nvar);
940 PP = Polyhedron2Param_Polyhedron(P, U, options);
942 for (nd = 0, PD = PP->D; PD; ++nd, PD = PD->next);
943 s = ALLOCN(struct evalue_section, nd);
945 Polyhedron *TC = true_context(P, U, MaxRays);
946 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, PD, rVD)
947 unsigned *facets;
949 facets = active_constraints(PP, PD);
951 Vector *inner = inner_point(rVD);
952 s[i].D = rVD;
954 if (nvar == 1) {
955 s[i].E = summate_over_1d_pdomain(e, PP, PD, inner->p+1, options);
956 } else if (nvar == 2) {
957 summator_2d s2d(e, PP, inner->p+1, rVD->Dimension);
959 s[i].E = s2d.summate_over_pdomain(PP, facets, PD, options);
962 Vector_Free(inner);
963 free(facets);
964 END_FORALL_REDUCED_DOMAIN
965 Polyhedron_Free(TC);
966 Polyhedron_Free(U);
967 Param_Polyhedron_Free(PP);
969 options->MaxRays = MaxRays;
971 res = evalue_from_section_array(s, nd);
972 free(s);
974 return res;