2 * Sum polynomial over integer points in polytope using local
3 * Euler-Maclaurin formula by Berline and Vergne.
6 #include <barvinok/options.h>
7 #include <barvinok/util.h>
9 #include "conversion.h"
10 #include "decomposer.h"
12 #include "lattice_point.h"
13 #include "param_util.h"
14 #include "reduce_domain.h"
16 /* Compute total degree in first nvar variables */
17 static unsigned total_degree(const evalue
*e
, unsigned nvar
)
22 if (value_notzero_p(e
->d
))
24 assert(e
->x
.p
->type
== polynomial
);
25 if (e
->x
.p
->pos
-1 >= nvar
)
28 max_degree
= total_degree(&e
->x
.p
->arr
[0], nvar
);
29 for (i
= 1; i
< e
->x
.p
->size
; ++i
) {
30 unsigned degree
= i
+ total_degree(&e
->x
.p
->arr
[i
], nvar
);
31 if (degree
> max_degree
)
44 Value
*factorial(unsigned n
)
50 int size
= 3*(n
+ 5)/2;
52 fact
.fact
= (Value
*)realloc(fact
.fact
, size
*sizeof(Value
));
55 for (int i
= fact
.n
; i
<= n
; ++i
) {
56 value_init(fact
.fact
[i
]);
58 value_set_si(fact
.fact
[0], 1);
60 mpz_mul_ui(fact
.fact
[i
], fact
.fact
[i
-1], i
);
73 Value
*binomial(unsigned n
, unsigned k
)
76 return &binom
.binom
[n
]->p
[k
];
78 if (n
>= binom
.size
) {
79 int size
= 3*(n
+ 5)/2;
81 binom
.binom
= (Vector
**)realloc(binom
.binom
, size
*sizeof(Vector
*));
84 for (int i
= binom
.n
; i
<= n
; ++i
) {
85 binom
.binom
[i
] = Vector_Alloc(i
+1);
87 value_set_si(binom
.binom
[0]->p
[0], 1);
89 value_set_si(binom
.binom
[i
]->p
[0], 1);
90 value_set_si(binom
.binom
[i
]->p
[i
], 1);
91 for (int j
= 1; j
< i
; ++j
)
92 value_addto(binom
.binom
[i
]->p
[j
],
93 binom
.binom
[i
-1]->p
[j
-1], binom
.binom
[i
-1]->p
[j
]);
97 return &binom
.binom
[n
]->p
[k
];
104 power(Value v
, int max
) {
105 powers
= Vector_Alloc(max
+1);
107 value_set_si(powers
->p
[0], 1);
109 value_assign(powers
->p
[1], v
);
115 Value
*operator[](int exp
) {
117 assert(exp
< powers
->Size
);
118 for (; n
<= exp
; ++n
)
119 value_multiply(powers
->p
[n
], powers
->p
[n
-1], powers
->p
[1]);
120 return &powers
->p
[exp
];
125 * Computes the coefficients of
127 * \mu(-t + R_+)(\xi) = \sum_{n=0)^\infty -b(n+1, t)/(n+1)! \xi^n
129 * where b(n, t) is the Bernoulli polynomial of degree n in t
130 * and t(p) is an expression (a fractional part) of the parameters p
131 * such that 0 <= t(p) < 1 for all values of the parameters.
132 * The coefficients are computed on demand up to (and including)
133 * the maximal degree max_degree.
137 evalue
**coefficients
;
140 mu_1d(unsigned max_degree
, evalue
*t
) : max_degree(max_degree
), t(t
) {
141 coefficients
= new evalue
*[max_degree
+1];
142 for (int i
= 0; i
< max_degree
+1; ++i
)
143 coefficients
[i
] = NULL
;
146 for (int i
= 0; i
< max_degree
+1; ++i
)
148 evalue_free(coefficients
[i
]);
149 delete [] coefficients
;
151 void compute_coefficient(unsigned n
);
152 const evalue
*coefficient(unsigned n
) {
153 if (!coefficients
[n
])
154 compute_coefficient(n
);
155 return coefficients
[n
];
159 void mu_1d::compute_coefficient(unsigned n
)
161 struct poly_list
*bernoulli
= bernoulli_compute(n
+1);
162 evalue
*c
= evalue_polynomial(bernoulli
->poly
[n
+1], t
);
164 evalue_div(c
, *factorial(n
+1));
170 * Computes the coefficients of
172 * \mu(a)(y) = \sum_{n_1} \sum_{n_2} c_{n_1,n_2} y^{n_1} y^{n_2}
174 * with c_{n1,n2} given by
176 * b(n1+1,t1)/(n1+1)! b(n2+1,t2)/(n2+1)!
177 * - b(n1+n2+2,t2)/(n1+n2+2)! (-c1)^{n1+1}
178 * - b(n1+n2+2,t1)/(n1+n2+2)! (-c2)^{n2+1}
180 * where b(n, t) is the Bernoulli polynomial of degree n in t,
181 * t1(p) and t2(p) are expressions (fractional parts) of the parameters p
182 * such that 0 <= t1(p), t2(p) < 1 for all values of the parameters
183 * and c1 = cn/c1d and c2 = cn/c2d.
184 * The coefficients are computed on demand up to (and including)
185 * the maximal degree (n1,n2) = (max_degree,max_degree).
187 * bernoulli_t[i][j] contains b(j+1, t_i)/(j+1)!
191 evalue
***coefficients
;
192 /* bernoulli_t[i][n] stores b(n+1, t_i)/(n+1)! */
193 evalue
**bernoulli_t
[2];
194 /* stores the powers of -cn */
195 struct power
*power_cn
;
196 struct power
*power_c1d
;
197 struct power
*power_c2d
;
200 mu_2d(unsigned max_degree
, evalue
*t1
, evalue
*t2
,
201 Value cn
, Value c1d
, Value c2d
) : max_degree(max_degree
) {
204 coefficients
= new evalue
**[max_degree
+1];
205 for (int i
= 0; i
< max_degree
+1; ++i
) {
206 coefficients
[i
] = new evalue
*[max_degree
+1];
207 for (int j
= 0; j
< max_degree
+1; ++j
)
208 coefficients
[i
][j
] = NULL
;
210 for (int i
= 0; i
< 2; ++i
) {
211 bernoulli_t
[i
] = new evalue
*[max_degree
+2];
212 for (int j
= 0; j
< max_degree
+2; ++j
)
213 bernoulli_t
[i
][j
] = NULL
;
215 value_oppose(cn
, cn
);
216 power_cn
= new struct power(cn
, max_degree
+1);
217 value_oppose(cn
, cn
);
218 power_c1d
= new struct power(c1d
, max_degree
+1);
219 power_c2d
= new struct power(c2d
, max_degree
+1);
222 for (int i
= 0; i
< max_degree
+1; ++i
) {
223 for (int j
= 0; j
< max_degree
+1; ++j
)
224 if (coefficients
[i
][j
])
225 evalue_free(coefficients
[i
][j
]);
226 delete [] coefficients
[i
];
228 delete [] coefficients
;
229 for (int i
= 0; i
< 2; ++i
)
230 for (int j
= 0; j
< max_degree
+2; ++j
)
231 if (bernoulli_t
[i
][j
])
232 evalue_free(bernoulli_t
[i
][j
]);
233 for (int i
= 0; i
< 2; ++i
)
234 delete [] bernoulli_t
[i
];
239 const evalue
*get_bernoulli(int n
, int i
);
241 void compute_coefficient(unsigned n1
, unsigned n2
);
242 const evalue
*coefficient(unsigned n1
, unsigned n2
) {
243 if (!coefficients
[n1
][n2
])
244 compute_coefficient(n1
, n2
);
245 return coefficients
[n1
][n2
];
250 * Returns b(n, t_i)/n!
252 const evalue
*mu_2d::get_bernoulli(int n
, int i
)
254 if (!bernoulli_t
[i
][n
-1]) {
255 struct poly_list
*bernoulli
= bernoulli_compute(n
);
256 bernoulli_t
[i
][n
-1] = evalue_polynomial(bernoulli
->poly
[n
], t
[i
]);
257 evalue_div(bernoulli_t
[i
][n
-1], *factorial(n
));
259 return bernoulli_t
[i
][n
-1];
262 void mu_2d::compute_coefficient(unsigned n1
, unsigned n2
)
264 evalue
*c
= evalue_dup(get_bernoulli(n1
+1, 0));
265 emul(get_bernoulli(n2
+1, 1), c
);
267 if (value_notzero_p(*(*power_cn
)[1])) {
271 value_init(neg_power
);
273 t
= evalue_dup(get_bernoulli(n1
+n2
+2, 1));
274 value_multiply(neg_power
,
275 *(*power_cn
)[n1
+1], *binomial(n1
+n2
+1, n1
+1));
276 value_oppose(neg_power
, neg_power
);
277 evalue_mul_div(t
, neg_power
, *(*power_c1d
)[n1
+1]);
281 t
= evalue_dup(get_bernoulli(n1
+n2
+2, 0));
282 value_multiply(neg_power
,
283 *(*power_cn
)[n2
+1], *binomial(n1
+n2
+1, n2
+1));
284 value_oppose(neg_power
, neg_power
);
285 evalue_mul_div(t
, neg_power
, *(*power_c2d
)[n2
+1]);
289 value_clear(neg_power
);
292 coefficients
[n1
][n2
] = c
;
295 /* Later: avoid recomputation of mu of faces that appear in
296 * more than one domain.
298 struct summator_2d
: public signed_cone_consumer
, public vertex_decomposer
{
299 const evalue
*polynomial
;
303 /* substitutions to use when result is 0-dimensional (only parameters) */
305 /* substitutions to use when result is 1-dimensional */
309 summator_2d(evalue
*e
, Param_Polyhedron
*PP
, Value
*inner
,
311 polynomial(e
), vertex_decomposer(PP
, *this),
312 inner(inner
), nparam(nparam
) {
315 subs_0d
= new evalue
*[2+nparam
];
316 subs_1d
= new evalue
*[2+nparam
];
321 for (int i
= 0; i
< nparam
; ++i
) {
322 subs_0d
[2+i
] = evalue_var(i
);
323 subs_1d
[2+i
] = evalue_var(1+i
);
327 for (int i
= 0; i
< nparam
; ++i
) {
328 evalue_free(subs_0d
[2+i
]);
329 evalue_free(subs_1d
[2+i
]);
334 evalue
*summate_over_pdomain(Polyhedron
*P
,
336 struct barvinok_options
*options
);
337 void handle_facet(Param_Polyhedron
*PP
, Param_Domain
*FD
, Value
*normal
);
338 void integrate(Polyhedron
*P
, Param_Domain
*PD
);
339 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
342 /* Replaces poly by its derivative along variable var */
343 static void evalue_derive(evalue
*poly
, int var
)
345 if (value_notzero_p(poly
->d
)) {
346 value_set_si(poly
->x
.n
, 0);
347 value_set_si(poly
->d
, 1);
350 assert(poly
->x
.p
->type
== polynomial
);
351 if (poly
->x
.p
->pos
-1 > var
) {
352 free_evalue_refs(poly
);
354 evalue_set_si(poly
, 0, 1);
358 if (poly
->x
.p
->pos
-1 < var
) {
359 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
)
360 evalue_derive(&poly
->x
.p
->arr
[i
], var
);
365 assert(poly
->x
.p
->size
> 1);
366 enode
*p
= poly
->x
.p
;
367 free_evalue_refs(&p
->arr
[0]);
369 value_clear(poly
->d
);
378 for (int i
= 0; i
< p
->size
; ++i
) {
379 value_set_si(factor
, i
+1);
380 p
->arr
[i
] = p
->arr
[i
+1];
381 evalue_mul(&p
->arr
[i
], factor
);
386 /* Check whether e is constant with respect to variable var */
387 static int evalue_is_constant(evalue
*e
, int var
)
391 if (value_notzero_p(e
->d
))
393 if (e
->x
.p
->type
== polynomial
&& e
->x
.p
->pos
-1 == var
)
395 assert(e
->x
.p
->type
== polynomial
||
396 e
->x
.p
->type
== fractional
||
397 e
->x
.p
->type
== flooring
);
398 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
399 if (!evalue_is_constant(&e
->x
.p
->arr
[i
], var
))
404 /* Replaces poly by its anti-derivative with constant 0 along variable var */
405 static void evalue_anti_derive(evalue
*poly
, int var
)
409 if (value_zero_p(poly
->d
) &&
410 poly
->x
.p
->type
== polynomial
&& poly
->x
.p
->pos
-1 < var
) {
411 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
)
412 evalue_anti_derive(&poly
->x
.p
->arr
[i
], var
);
417 if (evalue_is_constant(poly
, var
)) {
418 p
= new_enode(polynomial
, 2, 1+var
);
419 evalue_set_si(&p
->arr
[0], 0, 1);
420 value_clear(p
->arr
[1].d
);
426 assert(poly
->x
.p
->type
== polynomial
);
428 p
= new_enode(polynomial
, poly
->x
.p
->size
+1, 1+var
);
429 evalue_set_si(&p
->arr
[0], 0, 1);
430 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
) {
431 value_clear(p
->arr
[1+i
].d
);
432 p
->arr
[1+i
] = poly
->x
.p
->arr
[i
];
433 value_set_si(poly
->d
, 1+i
);
434 evalue_div(&p
->arr
[1+i
], poly
->d
);
438 value_set_si(poly
->d
, 0);
441 /* Computes offsets in the basis given by the rays of the cone
442 * to the integer point in the fundamental parallelepiped of
444 * The resulting evalues contain only the parameters as variables.
446 evalue
**offsets_to_integer_point(Matrix
*Rays
, Matrix
*vertex
)
448 unsigned nparam
= vertex
->NbColumns
-2;
449 evalue
**t
= new evalue
*[2];
451 if (value_one_p(vertex
->p
[0][nparam
+1])) {
452 t
[0] = evalue_zero();
453 t
[1] = evalue_zero();
455 Matrix
*R2
= Matrix_Copy(Rays
);
456 Matrix_Transposition(R2
);
457 Matrix
*inv
= Matrix_Alloc(Rays
->NbColumns
, Rays
->NbRows
);
458 int ok
= Matrix_Inverse(R2
, inv
);
462 /* We want the fractional part of the negative relative coordinates */
463 Vector_Oppose(inv
->p
[0], inv
->p
[0], inv
->NbColumns
);
464 Vector_Oppose(inv
->p
[1], inv
->p
[1], inv
->NbColumns
);
466 Matrix
*neg_rel
= Matrix_Alloc(2, nparam
+1);
468 Matrix_Product(inv
, vertex
, neg_rel
);
471 t
[0] = fractional_part(neg_rel
->p
[0], vertex
->p
[0][nparam
+1],
473 t
[1] = fractional_part(neg_rel
->p
[1], vertex
->p
[0][nparam
+1],
475 Matrix_Free(neg_rel
);
482 * Called from decompose_at_vertex.
484 * Handles a cone in the signed decomposition of the supporting
485 * cone of a vertex. The cone is assumed to be unimodular.
487 void summator_2d::handle(const signed_cone
& sc
, barvinok_options
*options
)
490 unsigned degree
= total_degree(polynomial
, 2);
492 subs_0d
[0] = affine2evalue(V
->Vertex
->p
[0],
493 V
->Vertex
->p
[0][nparam
+1], nparam
);
494 subs_0d
[1] = affine2evalue(V
->Vertex
->p
[1],
495 V
->Vertex
->p
[1][nparam
+1], nparam
);
499 assert(V
->Vertex
->NbRows
> 0);
500 Param_Vertex_Common_Denominator(V
);
502 Matrix
*Rays
= zz2matrix(sc
.rays
);
504 t
= offsets_to_integer_point(Rays
, V
->Vertex
);
506 Vector
*c
= Vector_Alloc(3);
507 Inner_Product(Rays
->p
[0], Rays
->p
[1], 2, &c
->p
[0]);
508 Inner_Product(Rays
->p
[0], Rays
->p
[0], 2, &c
->p
[1]);
509 Inner_Product(Rays
->p
[1], Rays
->p
[1], 2, &c
->p
[2]);
511 mu_2d
mu(degree
, t
[0], t
[1], c
->p
[0], c
->p
[1], c
->p
[2]);
514 struct power
power_r00(Rays
->p
[0][0], degree
);
515 struct power
power_r01(Rays
->p
[0][1], degree
);
516 struct power
power_r10(Rays
->p
[1][0], degree
);
517 struct power
power_r11(Rays
->p
[1][1], degree
);
519 Value factor
, tmp1
, tmp2
;
523 evalue
*res
= evalue_zero();
524 evalue
*dx1
= evalue_dup(polynomial
);
525 for (int i
= 0; !EVALUE_IS_ZERO(*dx1
); ++i
) {
526 evalue
*dx2
= evalue_dup(dx1
);
527 for (int j
= 0; !EVALUE_IS_ZERO(*dx2
); ++j
) {
528 evalue
*cij
= evalue_zero();
529 for (int n1
= 0; n1
<= i
+j
; ++n1
) {
531 value_set_si(factor
, 0);
532 for (int k
= max(0, i
-n2
); k
<= i
&& k
<= n1
; ++k
) {
534 value_multiply(tmp1
, *power_r00
[k
], *power_r01
[n1
-k
]);
535 value_multiply(tmp1
, tmp1
, *power_r10
[l
]);
536 value_multiply(tmp1
, tmp1
, *power_r11
[n2
-l
]);
537 value_multiply(tmp1
, tmp1
, *binomial(n1
, k
));
538 value_multiply(tmp1
, tmp1
, *binomial(n2
, l
));
539 value_addto(factor
, factor
, tmp1
);
541 if (value_zero_p(factor
))
544 evalue
*c
= evalue_dup(mu
.coefficient(n1
, n2
));
545 evalue_mul(c
, factor
);
549 evalue
*d
= evalue_dup(dx2
);
550 evalue_substitute(d
, subs_0d
);
555 evalue_derive(dx2
, 1);
558 evalue_derive(dx1
, 0);
561 for (int i
= 0; i
< 2; ++i
) {
562 evalue_free(subs_0d
[i
]);
578 evalue
*summator_2d::summate_over_pdomain(Polyhedron
*P
,
580 struct barvinok_options
*options
)
585 assert(PP
->V
->Vertex
->NbRows
== 2);
587 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, PP
) // _i internal counter
588 decompose_at_vertex(V
, _i
, options
);
589 END_FORALL_PVertex_in_ParamPolyhedron
;
591 Vector
*normal
= Vector_Alloc(2);
592 for (j
= P
->NbEq
; j
< P
->NbConstraints
; ++j
) {
594 Vector_Copy(P
->Constraint
[j
]+1, normal
->p
, 2);
595 if (value_zero_p(normal
->p
[0]) && value_zero_p(normal
->p
[1]))
598 FD
= Param_Polyhedron_Facet(PP
, PD
, P
, j
);
599 Vector_Normalize(normal
->p
, 2);
600 handle_facet(PP
, FD
, normal
->p
);
601 Param_Domain_Free(FD
);
610 void summator_2d::handle_facet(Param_Polyhedron
*PP
, Param_Domain
*FD
,
615 Param_Vertices
*vertex
[2];
617 unsigned degree
= total_degree(polynomial
, 2);
619 FORALL_PVertex_in_ParamPolyhedron(V
, FD
, PP
)
621 END_FORALL_PVertex_in_ParamPolyhedron
;
625 /* We can take either vertex[0] or vertex[1];
626 * the result should be the same
628 Param_Vertex_Common_Denominator(vertex
[0]);
630 /* The extra variable in front is the coordinate along the facet. */
631 Vector
*coef_normal
= Vector_Alloc(1 + nparam
+ 2);
632 Vector_Combine(vertex
[0]->Vertex
->p
[0], vertex
[0]->Vertex
->p
[1],
633 coef_normal
->p
+1, normal
[0], normal
[1], nparam
+1);
634 value_assign(coef_normal
->p
[1+nparam
+1], vertex
[0]->Vertex
->p
[0][nparam
+1]);
635 Vector_Normalize(coef_normal
->p
, coef_normal
->Size
);
637 Vector
*base
= Vector_Alloc(2);
638 value_oppose(base
->p
[0], normal
[1]);
639 value_assign(base
->p
[1], normal
[0]);
642 Inner_Product(normal
, normal
, 2, &det
);
644 Vector
*s
= Vector_Alloc(1+nparam
+2);
645 value_multiply(s
->p
[1+nparam
+1], coef_normal
->p
[1+nparam
+1], det
);
646 value_multiply(s
->p
[0], base
->p
[0], s
->p
[1+nparam
+1]);
647 Vector_Scale(coef_normal
->p
+1, s
->p
+1, normal
[0], nparam
+1);
648 subs_1d
[0] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
649 value_multiply(s
->p
[0], base
->p
[1], s
->p
[1+nparam
+1]);
650 Vector_Scale(coef_normal
->p
+1, s
->p
+1, normal
[1], nparam
+1);
651 subs_1d
[1] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
655 if (value_one_p(coef_normal
->p
[coef_normal
->Size
-1]))
658 Vector_Oppose(coef_normal
->p
+1, coef_normal
->p
+1, nparam
+1);
659 t
= fractional_part(coef_normal
->p
,
660 coef_normal
->p
[coef_normal
->Size
-1],
663 Vector_Free(coef_normal
);
667 struct power
power_normal0(normal
[0], degree
);
668 struct power
power_normal1(normal
[1], degree
);
669 struct power
power_det(det
, degree
);
673 evalue
*res
= evalue_zero();
674 evalue
*dx1
= evalue_dup(polynomial
);
675 for (int i
= 0; !EVALUE_IS_ZERO(*dx1
); ++i
) {
676 evalue
*dx2
= evalue_dup(dx1
);
677 for (int j
= 0; !EVALUE_IS_ZERO(*dx2
); ++j
) {
678 value_multiply(factor
, *power_normal0
[i
], *power_normal1
[j
]);
679 if (value_notzero_p(factor
)) {
680 value_multiply(factor
, factor
, *binomial(i
+j
, i
));
682 evalue
*c
= evalue_dup(mu
.coefficient(i
+j
));
683 evalue_mul_div(c
, factor
, *power_det
[i
+j
]);
685 evalue
*d
= evalue_dup(dx2
);
686 evalue_substitute(d
, subs_1d
);
692 evalue_derive(dx2
, 1);
695 evalue_derive(dx1
, 0);
699 for (int i
= 0; i
< 2; ++i
) {
700 evalue_free(subs_1d
[i
]);
704 evalue_anti_derive(res
, 0);
706 Matrix
*z
= Matrix_Alloc(2, nparam
+2);
707 Vector
*fixed_z
= Vector_Alloc(2);
708 for (int i
= 0; i
< 2; ++i
) {
709 Vector_Combine(vertex
[i
]->Vertex
->p
[0], vertex
[i
]->Vertex
->p
[1],
710 z
->p
[i
], base
->p
[0], base
->p
[1], nparam
+1);
711 value_multiply(z
->p
[i
][nparam
+1],
712 det
, vertex
[i
]->Vertex
->p
[0][nparam
+1]);
713 Inner_Product(z
->p
[i
], inner
, nparam
+1, &fixed_z
->p
[i
]);
716 /* Put on a common denominator */
717 value_multiply(fixed_z
->p
[0], fixed_z
->p
[0], z
->p
[1][nparam
+1]);
718 value_multiply(fixed_z
->p
[1], fixed_z
->p
[1], z
->p
[0][nparam
+1]);
719 /* Make sure z->p[0] is smaller than z->p[1] (for an internal
720 * point of the chamber and hence for all parameter values in
721 * the chamber), to ensure we integrate in the right direction.
723 if (value_lt(fixed_z
->p
[1], fixed_z
->p
[0]))
724 Vector_Exchange(z
->p
[0], z
->p
[1], nparam
+2);
725 Vector_Free(fixed_z
);
728 subs_0d
[1] = affine2evalue(z
->p
[1], z
->p
[1][nparam
+1], nparam
);
729 evalue
*up
= evalue_dup(res
);
730 evalue_substitute(up
, subs_0d
+1);
731 evalue_free(subs_0d
[1]);
733 subs_0d
[1] = affine2evalue(z
->p
[0], z
->p
[0][nparam
+1], nparam
);
734 evalue_substitute(res
, subs_0d
+1);
737 evalue_free(subs_0d
[1]);
747 /* Integrate the polynomial over the whole polygon using
748 * the Green-Stokes theorem.
750 void summator_2d::integrate(Polyhedron
*P
, Param_Domain
*PD
)
753 evalue
*res
= evalue_zero();
755 evalue
*I
= evalue_dup(polynomial
);
756 evalue_anti_derive(I
, 0);
759 Vector
*normal
= Vector_Alloc(2);
760 Vector
*dir
= Vector_Alloc(2);
761 Matrix
*v0v1
= Matrix_Alloc(2, nparam
+2);
762 Vector
*f_v0v1
= Vector_Alloc(2);
763 Vector
*s
= Vector_Alloc(1+nparam
+2);
764 for (int j
= P
->NbEq
; j
< P
->NbConstraints
; ++j
) {
767 Param_Vertices
*vertex
[2];
768 Vector_Copy(P
->Constraint
[j
]+1, normal
->p
, 2);
770 if (value_zero_p(normal
->p
[0]))
773 Vector_Normalize(normal
->p
, 2);
774 value_assign(dir
->p
[0], normal
->p
[1]);
775 value_oppose(dir
->p
[1], normal
->p
[0]);
777 FD
= Param_Polyhedron_Facet(PP
, PD
, P
, j
);
779 FORALL_PVertex_in_ParamPolyhedron(V
, FD
, PP
)
781 END_FORALL_PVertex_in_ParamPolyhedron
;
785 Param_Vertex_Common_Denominator(vertex
[0]);
786 Param_Vertex_Common_Denominator(vertex
[1]);
788 value_oppose(tmp
, vertex
[1]->Vertex
->p
[0][nparam
+1]);
789 for (int i
= 0; i
< 2; ++i
)
790 Vector_Combine(vertex
[1]->Vertex
->p
[i
],
791 vertex
[0]->Vertex
->p
[i
],
793 vertex
[0]->Vertex
->p
[0][nparam
+1], tmp
, nparam
+1);
794 value_multiply(v0v1
->p
[0][nparam
+1],
795 vertex
[0]->Vertex
->p
[0][nparam
+1],
796 vertex
[1]->Vertex
->p
[0][nparam
+1]);
797 value_assign(v0v1
->p
[1][nparam
+1], v0v1
->p
[0][nparam
+1]);
799 /* Order vertices to ensure we integrate in the right
800 * direction, i.e., with the polytope "on the left".
802 for (int i
= 0; i
< 2; ++i
)
803 Inner_Product(v0v1
->p
[i
], inner
, nparam
+1, &f_v0v1
->p
[i
]);
805 Inner_Product(dir
->p
, f_v0v1
->p
, 2, &tmp
);
806 if (value_neg_p(tmp
)) {
807 Param_Vertices
*PV
= vertex
[0];
808 vertex
[0] = vertex
[1];
810 for (int i
= 0; i
< 2; ++i
)
811 Vector_Oppose(v0v1
->p
[i
], v0v1
->p
[i
], nparam
+1);
813 value_oppose(tmp
, normal
->p
[0]);
814 if (value_neg_p(tmp
)) {
815 value_oppose(tmp
, tmp
);
816 Vector_Oppose(v0v1
->p
[1], v0v1
->p
[1], nparam
+1);
818 value_multiply(tmp
, tmp
, v0v1
->p
[1][nparam
+1]);
819 evalue
*top
= affine2evalue(v0v1
->p
[1], tmp
, nparam
);
821 value_multiply(s
->p
[0], normal
->p
[1], vertex
[0]->Vertex
->p
[0][nparam
+1]);
822 Vector_Copy(vertex
[0]->Vertex
->p
[0], s
->p
+1, nparam
+2);
823 subs_1d
[0] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
824 value_multiply(s
->p
[0], normal
->p
[0], vertex
[0]->Vertex
->p
[0][nparam
+1]);
825 value_oppose(s
->p
[0], s
->p
[0]);
826 Vector_Copy(vertex
[0]->Vertex
->p
[1], s
->p
+1, nparam
+2);
827 subs_1d
[1] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
829 evalue
*d
= evalue_dup(I
);
830 evalue_substitute(d
, subs_1d
);
831 evalue_anti_derive(d
, 0);
833 evalue_free(subs_1d
[0]);
834 evalue_free(subs_1d
[1]);
837 evalue_substitute(d
, subs_0d
+1);
838 evalue_mul(d
, dir
->p
[1]);
839 evalue_free(subs_0d
[1]);
844 Param_Domain_Free(FD
);
858 evalue
*summate_over_1d_pdomain(evalue
*e
,
859 Param_Polyhedron
*PP
, Param_Domain
*PD
,
860 Polyhedron
*P
, Value
*inner
,
861 struct barvinok_options
*options
)
865 Param_Vertices
*vertex
[2];
866 unsigned nparam
= PP
->V
->Vertex
->NbColumns
-2;
867 evalue
*subs_0d
[1+nparam
];
870 unsigned degree
= total_degree(e
, 1);
872 for (int i
= 0; i
< nparam
; ++i
)
873 subs_0d
[1+i
] = evalue_var(i
);
875 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, PP
)
877 END_FORALL_PVertex_in_ParamPolyhedron
;
880 Vector
*fixed
= Vector_Alloc(2);
881 for (int i
= 0; i
< 2; ++i
) {
882 Inner_Product(vertex
[i
]->Vertex
->p
[0], inner
, nparam
+1, &fixed
->p
[i
]);
883 value_multiply(fixed
->p
[i
],
884 fixed
->p
[i
], vertex
[1-i
]->Vertex
->p
[0][nparam
+1]);
886 if (value_lt(fixed
->p
[1], fixed
->p
[0])) {
888 vertex
[0] = vertex
[1];
893 Vector
*coef
= Vector_Alloc(nparam
+1);
894 for (int i
= 0; i
< 2; ++i
)
895 a
[i
] = affine2evalue(vertex
[i
]->Vertex
->p
[0],
896 vertex
[i
]->Vertex
->p
[0][nparam
+1], nparam
);
897 if (value_one_p(vertex
[0]->Vertex
->p
[0][nparam
+1]))
898 t
[0] = evalue_zero();
900 Vector_Oppose(vertex
[0]->Vertex
->p
[0], coef
->p
, nparam
+1);
901 t
[0] = fractional_part(coef
->p
, vertex
[0]->Vertex
->p
[0][nparam
+1],
904 if (value_one_p(vertex
[1]->Vertex
->p
[0][nparam
+1]))
905 t
[1] = evalue_zero();
907 Vector_Copy(vertex
[1]->Vertex
->p
[0], coef
->p
, nparam
+1);
908 t
[1] = fractional_part(coef
->p
, vertex
[1]->Vertex
->p
[0][nparam
+1],
913 evalue
*I
= evalue_dup(e
);
914 evalue_anti_derive(I
, 0);
915 evalue
*up
= evalue_dup(I
);
917 evalue_substitute(up
, subs_0d
);
920 evalue_substitute(I
, subs_0d
);
927 mu_1d
mu0(degree
, t
[0]);
928 mu_1d
mu1(degree
, t
[1]);
930 evalue
*dx
= evalue_dup(e
);
931 for (int n
= 0; !EVALUE_IS_ZERO(*dx
); ++n
) {
936 evalue_substitute(d
, subs_0d
);
937 emul(mu0
.coefficient(n
), d
);
943 evalue_substitute(d
, subs_0d
);
944 emul(mu1
.coefficient(n
), d
);
950 evalue_derive(dx
, 0);
954 for (int i
= 0; i
< nparam
; ++i
)
955 evalue_free(subs_0d
[1+i
]);
957 for (int i
= 0; i
< 2; ++i
) {
965 static evalue
*summate_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
966 struct barvinok_options
*options
)
969 Param_Polyhedron
*PP
;
973 struct evalue_section
*s
;
976 MaxRays
= options
->MaxRays
;
977 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
979 U
= Universe_Polyhedron(D
->Dimension
- nvar
);
980 PP
= Polyhedron2Param_Polyhedron(D
, U
, options
);
982 for (nd
= 0, PD
= PP
->D
; PD
; ++nd
, PD
= PD
->next
);
983 s
= ALLOCN(struct evalue_section
, nd
);
985 Polyhedron
*TC
= true_context(D
, U
, MaxRays
);
986 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, PD
, rVD
)
989 CA
= align_context(PD
->Domain
, D
->Dimension
, options
->MaxRays
);
990 F
= DomainIntersection(D
, CA
, options
->MaxRays
);
993 Vector
*inner
= inner_point(rVD
);
997 s
[i
].E
= summate_over_1d_pdomain(e
, PP
, PD
, F
, inner
->p
+1, options
);
998 } else if (nvar
== 2) {
999 summator_2d
s2d(e
, PP
, inner
->p
+1, rVD
->Dimension
);
1001 s
[i
].E
= s2d
.summate_over_pdomain(F
, PD
, options
);
1006 END_FORALL_REDUCED_DOMAIN
1007 Polyhedron_Free(TC
);
1009 Param_Polyhedron_Free(PP
);
1011 options
->MaxRays
= MaxRays
;
1013 res
= evalue_from_section_array(s
, nd
);
1019 evalue
*euler_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
1027 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
1028 return evalue_dup(e
);
1030 assert(value_zero_p(e
->d
));
1031 assert(e
->x
.p
->type
== partition
);
1033 res
= evalue_zero();
1035 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
1037 t
= summate_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
1038 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), options
);