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 /* returns an evalue that corresponds to
299 static evalue
*term(int param
)
301 evalue
*EP
= new evalue();
303 value_set_si(EP
->d
,0);
304 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
305 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
306 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
310 /* Later: avoid recomputation of mu of faces that appear in
311 * more than one domain.
313 struct summator_2d
: public signed_cone_consumer
, public vertex_decomposer
{
314 const evalue
*polynomial
;
318 /* substitutions to use when result is 0-dimensional (only parameters) */
320 /* substitutions to use when result is 1-dimensional */
324 summator_2d(evalue
*e
, Polyhedron
*P
, unsigned nbV
, Value
*inner
,
326 polynomial(e
), vertex_decomposer(P
, nbV
, *this),
327 inner(inner
), nparam(nparam
) {
330 subs_0d
= new evalue
*[2+nparam
];
331 subs_1d
= new evalue
*[2+nparam
];
336 for (int i
= 0; i
< nparam
; ++i
) {
337 subs_0d
[2+i
] = term(i
);
338 subs_1d
[2+i
] = term(1+i
);
342 for (int i
= 0; i
< nparam
; ++i
) {
343 free_evalue_refs(subs_0d
[2+i
]);
345 free_evalue_refs(subs_1d
[2+i
]);
351 evalue
*summate_over_pdomain(Param_Polyhedron
*PP
,
353 struct barvinok_options
*options
);
354 void handle_facet(Param_Polyhedron
*PP
, Param_Domain
*FD
, Value
*normal
);
355 void integrate(Param_Polyhedron
*PP
, Param_Domain
*PD
);
356 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
359 /* Replaces poly by its derivative along variable var */
360 static void evalue_derive(evalue
*poly
, int var
)
362 if (value_notzero_p(poly
->d
)) {
363 value_set_si(poly
->x
.n
, 0);
364 value_set_si(poly
->d
, 1);
367 assert(poly
->x
.p
->type
== polynomial
);
368 if (poly
->x
.p
->pos
-1 > var
) {
369 free_evalue_refs(poly
);
371 evalue_set_si(poly
, 0, 1);
375 if (poly
->x
.p
->pos
-1 < var
) {
376 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
)
377 evalue_derive(&poly
->x
.p
->arr
[i
], var
);
382 assert(poly
->x
.p
->size
> 1);
383 enode
*p
= poly
->x
.p
;
384 free_evalue_refs(&p
->arr
[0]);
386 value_clear(poly
->d
);
395 for (int i
= 0; i
< p
->size
; ++i
) {
396 value_set_si(factor
, i
+1);
397 p
->arr
[i
] = p
->arr
[i
+1];
398 evalue_mul(&p
->arr
[i
], factor
);
403 /* Check whether e is constant with respect to variable var */
404 static int evalue_is_constant(evalue
*e
, int var
)
408 if (value_notzero_p(e
->d
))
410 if (e
->x
.p
->type
== polynomial
&& e
->x
.p
->pos
-1 == var
)
412 assert(e
->x
.p
->type
== polynomial
||
413 e
->x
.p
->type
== fractional
||
414 e
->x
.p
->type
== flooring
);
415 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
416 if (!evalue_is_constant(&e
->x
.p
->arr
[i
], var
))
421 /* Replaces poly by its anti-derivative with constant 0 along variable var */
422 static void evalue_anti_derive(evalue
*poly
, int var
)
426 if (value_zero_p(poly
->d
) &&
427 poly
->x
.p
->type
== polynomial
&& poly
->x
.p
->pos
-1 < var
) {
428 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
)
429 evalue_anti_derive(&poly
->x
.p
->arr
[i
], var
);
434 if (evalue_is_constant(poly
, var
)) {
435 p
= new_enode(polynomial
, 2, 1+var
);
436 evalue_set_si(&p
->arr
[0], 0, 1);
437 value_clear(p
->arr
[1].d
);
443 assert(poly
->x
.p
->type
== polynomial
);
445 p
= new_enode(polynomial
, poly
->x
.p
->size
+1, 1+var
);
446 evalue_set_si(&p
->arr
[0], 0, 1);
447 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
) {
448 value_clear(p
->arr
[1+i
].d
);
449 p
->arr
[1+i
] = poly
->x
.p
->arr
[i
];
450 value_set_si(poly
->d
, 1+i
);
451 evalue_div(&p
->arr
[1+i
], poly
->d
);
455 value_set_si(poly
->d
, 0);
458 /* Computes offsets in the basis given by the rays of the cone
459 * to the integer point in the fundamental parallelepiped of
461 * The resulting evalues contain only the parameters as variables.
463 evalue
**offsets_to_integer_point(Matrix
*Rays
, Matrix
*vertex
)
465 unsigned nparam
= vertex
->NbColumns
-2;
466 evalue
**t
= new evalue
*[2];
468 if (value_one_p(vertex
->p
[0][nparam
+1])) {
469 t
[0] = evalue_zero();
470 t
[1] = evalue_zero();
472 Matrix
*R2
= Matrix_Copy(Rays
);
473 Matrix_Transposition(R2
);
474 Matrix
*inv
= Matrix_Alloc(Rays
->NbColumns
, Rays
->NbRows
);
475 int ok
= Matrix_Inverse(R2
, inv
);
479 /* We want the fractional part of the negative relative coordinates */
480 Vector_Oppose(inv
->p
[0], inv
->p
[0], inv
->NbColumns
);
481 Vector_Oppose(inv
->p
[1], inv
->p
[1], inv
->NbColumns
);
483 Matrix
*neg_rel
= Matrix_Alloc(2, nparam
+1);
485 Matrix_Product(inv
, vertex
, neg_rel
);
488 t
[0] = fractional_part(neg_rel
->p
[0], vertex
->p
[0][nparam
+1],
490 t
[1] = fractional_part(neg_rel
->p
[1], vertex
->p
[0][nparam
+1],
492 Matrix_Free(neg_rel
);
499 * Called from decompose_at_vertex.
501 * Handles a cone in the signed decomposition of the supporting
502 * cone of a vertex. The cone is assumed to be unimodular.
504 void summator_2d::handle(const signed_cone
& sc
, barvinok_options
*options
)
507 unsigned degree
= total_degree(polynomial
, 2);
509 subs_0d
[0] = affine2evalue(V
->Vertex
->p
[0],
510 V
->Vertex
->p
[0][nparam
+1], nparam
);
511 subs_0d
[1] = affine2evalue(V
->Vertex
->p
[1],
512 V
->Vertex
->p
[1][nparam
+1], nparam
);
516 assert(V
->Vertex
->NbRows
> 0);
517 Param_Vertex_Common_Denominator(V
);
519 Matrix
*Rays
= zz2matrix(sc
.rays
);
521 t
= offsets_to_integer_point(Rays
, V
->Vertex
);
523 Vector
*c
= Vector_Alloc(3);
524 Inner_Product(Rays
->p
[0], Rays
->p
[1], 2, &c
->p
[0]);
525 Inner_Product(Rays
->p
[0], Rays
->p
[0], 2, &c
->p
[1]);
526 Inner_Product(Rays
->p
[1], Rays
->p
[1], 2, &c
->p
[2]);
528 mu_2d
mu(degree
, t
[0], t
[1], c
->p
[0], c
->p
[1], c
->p
[2]);
531 struct power
power_r00(Rays
->p
[0][0], degree
);
532 struct power
power_r01(Rays
->p
[0][1], degree
);
533 struct power
power_r10(Rays
->p
[1][0], degree
);
534 struct power
power_r11(Rays
->p
[1][1], degree
);
536 Value factor
, tmp1
, tmp2
;
540 evalue
*res
= evalue_zero();
541 evalue
*dx1
= evalue_dup(polynomial
);
542 for (int i
= 0; !EVALUE_IS_ZERO(*dx1
); ++i
) {
543 evalue
*dx2
= evalue_dup(dx1
);
544 for (int j
= 0; !EVALUE_IS_ZERO(*dx2
); ++j
) {
545 evalue
*cij
= evalue_zero();
546 for (int n1
= 0; n1
<= i
+j
; ++n1
) {
548 value_set_si(factor
, 0);
549 for (int k
= max(0, i
-n2
); k
<= i
&& k
<= n1
; ++k
) {
551 value_multiply(tmp1
, *power_r00
[k
], *power_r01
[n1
-k
]);
552 value_multiply(tmp1
, tmp1
, *power_r10
[l
]);
553 value_multiply(tmp1
, tmp1
, *power_r11
[n2
-l
]);
554 value_multiply(tmp1
, tmp1
, *binomial(n1
, k
));
555 value_multiply(tmp1
, tmp1
, *binomial(n2
, l
));
556 value_addto(factor
, factor
, tmp1
);
558 if (value_zero_p(factor
))
561 evalue
*c
= evalue_dup(mu
.coefficient(n1
, n2
));
562 evalue_mul(c
, factor
);
566 evalue
*d
= evalue_dup(dx2
);
567 evalue_substitute(d
, subs_0d
);
572 evalue_derive(dx2
, 1);
575 evalue_derive(dx1
, 0);
578 for (int i
= 0; i
< 2; ++i
) {
579 evalue_free(subs_0d
[i
]);
595 evalue
*summator_2d::summate_over_pdomain(Param_Polyhedron
*PP
,
597 struct barvinok_options
*options
)
602 assert(PP
->V
->Vertex
->NbRows
== 2);
604 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, PP
) // _i internal counter
605 decompose_at_vertex(V
, _i
, options
);
606 END_FORALL_PVertex_in_ParamPolyhedron
;
608 Vector
*normal
= Vector_Alloc(2);
609 for (j
= P
->NbEq
; j
< P
->NbConstraints
; ++j
) {
611 Vector_Copy(P
->Constraint
[j
]+1, normal
->p
, 2);
612 if (value_zero_p(normal
->p
[0]) && value_zero_p(normal
->p
[1]))
615 FD
= Param_Polyhedron_Facet(PP
, PD
, P
, j
);
616 Vector_Normalize(normal
->p
, 2);
617 handle_facet(PP
, FD
, normal
->p
);
618 Param_Domain_Free(FD
);
627 void summator_2d::handle_facet(Param_Polyhedron
*PP
, Param_Domain
*FD
,
632 Param_Vertices
*vertex
[2];
634 unsigned degree
= total_degree(polynomial
, 2);
636 FORALL_PVertex_in_ParamPolyhedron(V
, FD
, PP
)
638 END_FORALL_PVertex_in_ParamPolyhedron
;
642 /* We can take either vertex[0] or vertex[1];
643 * the result should be the same
645 Param_Vertex_Common_Denominator(vertex
[0]);
647 /* The extra variable in front is the coordinate along the facet. */
648 Vector
*coef_normal
= Vector_Alloc(1 + nparam
+ 2);
649 Vector_Combine(vertex
[0]->Vertex
->p
[0], vertex
[0]->Vertex
->p
[1],
650 coef_normal
->p
+1, normal
[0], normal
[1], nparam
+1);
651 value_assign(coef_normal
->p
[1+nparam
+1], vertex
[0]->Vertex
->p
[0][nparam
+1]);
652 Vector_Normalize(coef_normal
->p
, coef_normal
->Size
);
654 Vector
*base
= Vector_Alloc(2);
655 value_oppose(base
->p
[0], normal
[1]);
656 value_assign(base
->p
[1], normal
[0]);
659 Inner_Product(normal
, normal
, 2, &det
);
661 Vector
*s
= Vector_Alloc(1+nparam
+2);
662 value_multiply(s
->p
[1+nparam
+1], coef_normal
->p
[1+nparam
+1], det
);
663 value_multiply(s
->p
[0], base
->p
[0], s
->p
[1+nparam
+1]);
664 Vector_Scale(coef_normal
->p
+1, s
->p
+1, normal
[0], nparam
+1);
665 subs_1d
[0] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
666 value_multiply(s
->p
[0], base
->p
[1], s
->p
[1+nparam
+1]);
667 Vector_Scale(coef_normal
->p
+1, s
->p
+1, normal
[1], nparam
+1);
668 subs_1d
[1] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
672 if (value_one_p(coef_normal
->p
[coef_normal
->Size
-1]))
675 Vector_Oppose(coef_normal
->p
+1, coef_normal
->p
+1, nparam
+1);
676 t
= fractional_part(coef_normal
->p
,
677 coef_normal
->p
[coef_normal
->Size
-1],
680 Vector_Free(coef_normal
);
684 struct power
power_normal0(normal
[0], degree
);
685 struct power
power_normal1(normal
[1], degree
);
686 struct power
power_det(det
, degree
);
690 evalue
*res
= evalue_zero();
691 evalue
*dx1
= evalue_dup(polynomial
);
692 for (int i
= 0; !EVALUE_IS_ZERO(*dx1
); ++i
) {
693 evalue
*dx2
= evalue_dup(dx1
);
694 for (int j
= 0; !EVALUE_IS_ZERO(*dx2
); ++j
) {
695 value_multiply(factor
, *power_normal0
[i
], *power_normal1
[j
]);
696 if (value_notzero_p(factor
)) {
697 value_multiply(factor
, factor
, *binomial(i
+j
, i
));
699 evalue
*c
= evalue_dup(mu
.coefficient(i
+j
));
700 evalue_mul_div(c
, factor
, *power_det
[i
+j
]);
702 evalue
*d
= evalue_dup(dx2
);
703 evalue_substitute(d
, subs_1d
);
709 evalue_derive(dx2
, 1);
712 evalue_derive(dx1
, 0);
716 for (int i
= 0; i
< 2; ++i
) {
717 evalue_free(subs_1d
[i
]);
721 evalue_anti_derive(res
, 0);
723 Matrix
*z
= Matrix_Alloc(2, nparam
+2);
724 Vector
*fixed_z
= Vector_Alloc(2);
725 for (int i
= 0; i
< 2; ++i
) {
726 Vector_Combine(vertex
[i
]->Vertex
->p
[0], vertex
[i
]->Vertex
->p
[1],
727 z
->p
[i
], base
->p
[0], base
->p
[1], nparam
+1);
728 value_multiply(z
->p
[i
][nparam
+1],
729 det
, vertex
[i
]->Vertex
->p
[0][nparam
+1]);
730 Inner_Product(z
->p
[i
], inner
, nparam
+1, &fixed_z
->p
[i
]);
733 /* Put on a common denominator */
734 value_multiply(fixed_z
->p
[0], fixed_z
->p
[0], z
->p
[1][nparam
+1]);
735 value_multiply(fixed_z
->p
[1], fixed_z
->p
[1], z
->p
[0][nparam
+1]);
736 /* Make sure z->p[0] is smaller than z->p[1] (for an internal
737 * point of the chamber and hence for all parameter values in
738 * the chamber), to ensure we integrate in the right direction.
740 if (value_lt(fixed_z
->p
[1], fixed_z
->p
[0]))
741 Vector_Exchange(z
->p
[0], z
->p
[1], nparam
+2);
742 Vector_Free(fixed_z
);
745 subs_0d
[1] = affine2evalue(z
->p
[1], z
->p
[1][nparam
+1], nparam
);
746 evalue
*up
= evalue_dup(res
);
747 evalue_substitute(up
, subs_0d
+1);
748 evalue_free(subs_0d
[1]);
750 subs_0d
[1] = affine2evalue(z
->p
[0], z
->p
[0][nparam
+1], nparam
);
751 evalue_substitute(res
, subs_0d
+1);
754 evalue_free(subs_0d
[1]);
764 /* Integrate the polynomial over the whole polygon using
765 * the Green-Stokes theorem.
767 void summator_2d::integrate(Param_Polyhedron
*PP
, Param_Domain
*PD
)
770 evalue
*res
= evalue_zero();
772 evalue
*I
= evalue_dup(polynomial
);
773 evalue_anti_derive(I
, 0);
776 Vector
*normal
= Vector_Alloc(2);
777 Vector
*dir
= Vector_Alloc(2);
778 Matrix
*v0v1
= Matrix_Alloc(2, nparam
+2);
779 Vector
*f_v0v1
= Vector_Alloc(2);
780 Vector
*s
= Vector_Alloc(1+nparam
+2);
781 for (int j
= P
->NbEq
; j
< P
->NbConstraints
; ++j
) {
784 Param_Vertices
*vertex
[2];
785 Vector_Copy(P
->Constraint
[j
]+1, normal
->p
, 2);
787 if (value_zero_p(normal
->p
[0]))
790 Vector_Normalize(normal
->p
, 2);
791 value_assign(dir
->p
[0], normal
->p
[1]);
792 value_oppose(dir
->p
[1], normal
->p
[0]);
794 FD
= Param_Polyhedron_Facet(PP
, PD
, P
, j
);
796 FORALL_PVertex_in_ParamPolyhedron(V
, FD
, PP
)
798 END_FORALL_PVertex_in_ParamPolyhedron
;
802 Param_Vertex_Common_Denominator(vertex
[0]);
803 Param_Vertex_Common_Denominator(vertex
[1]);
805 value_oppose(tmp
, vertex
[1]->Vertex
->p
[0][nparam
+1]);
806 for (int i
= 0; i
< 2; ++i
)
807 Vector_Combine(vertex
[1]->Vertex
->p
[i
],
808 vertex
[0]->Vertex
->p
[i
],
810 vertex
[0]->Vertex
->p
[0][nparam
+1], tmp
, nparam
+1);
811 value_multiply(v0v1
->p
[0][nparam
+1],
812 vertex
[0]->Vertex
->p
[0][nparam
+1],
813 vertex
[1]->Vertex
->p
[0][nparam
+1]);
814 value_assign(v0v1
->p
[1][nparam
+1], v0v1
->p
[0][nparam
+1]);
816 /* Order vertices to ensure we integrate in the right
817 * direction, i.e., with the polytope "on the left".
819 for (int i
= 0; i
< 2; ++i
)
820 Inner_Product(v0v1
->p
[i
], inner
, nparam
+1, &f_v0v1
->p
[i
]);
822 Inner_Product(dir
->p
, f_v0v1
->p
, 2, &tmp
);
823 if (value_neg_p(tmp
)) {
824 Param_Vertices
*PV
= vertex
[0];
825 vertex
[0] = vertex
[1];
827 for (int i
= 0; i
< 2; ++i
)
828 Vector_Oppose(v0v1
->p
[i
], v0v1
->p
[i
], nparam
+1);
830 value_oppose(tmp
, normal
->p
[0]);
831 if (value_neg_p(tmp
)) {
832 value_oppose(tmp
, tmp
);
833 Vector_Oppose(v0v1
->p
[1], v0v1
->p
[1], nparam
+1);
835 value_multiply(tmp
, tmp
, v0v1
->p
[1][nparam
+1]);
836 evalue
*top
= affine2evalue(v0v1
->p
[1], tmp
, nparam
);
838 value_multiply(s
->p
[0], normal
->p
[1], vertex
[0]->Vertex
->p
[0][nparam
+1]);
839 Vector_Copy(vertex
[0]->Vertex
->p
[0], s
->p
+1, nparam
+2);
840 subs_1d
[0] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
841 value_multiply(s
->p
[0], normal
->p
[0], vertex
[0]->Vertex
->p
[0][nparam
+1]);
842 value_oppose(s
->p
[0], s
->p
[0]);
843 Vector_Copy(vertex
[0]->Vertex
->p
[1], s
->p
+1, nparam
+2);
844 subs_1d
[1] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
846 evalue
*d
= evalue_dup(I
);
847 evalue_substitute(d
, subs_1d
);
848 evalue_anti_derive(d
, 0);
850 evalue_free(subs_1d
[0]);
851 evalue_free(subs_1d
[1]);
854 evalue_substitute(d
, subs_0d
+1);
855 evalue_mul(d
, dir
->p
[1]);
856 evalue_free(subs_0d
[1]);
861 Param_Domain_Free(FD
);
875 evalue
*summate_over_1d_pdomain(evalue
*e
,
876 Param_Polyhedron
*PP
, Param_Domain
*PD
,
877 Polyhedron
*P
, Value
*inner
,
878 struct barvinok_options
*options
)
882 Param_Vertices
*vertex
[2];
883 unsigned nparam
= PP
->V
->Vertex
->NbColumns
-2;
884 evalue
*subs_0d
[1+nparam
];
887 unsigned degree
= total_degree(e
, 1);
889 for (int i
= 0; i
< nparam
; ++i
)
890 subs_0d
[1+i
] = term(i
);
892 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, PP
)
894 END_FORALL_PVertex_in_ParamPolyhedron
;
897 Vector
*fixed
= Vector_Alloc(2);
898 for (int i
= 0; i
< 2; ++i
) {
899 Inner_Product(vertex
[i
]->Vertex
->p
[0], inner
, nparam
+1, &fixed
->p
[i
]);
900 value_multiply(fixed
->p
[i
],
901 fixed
->p
[i
], vertex
[1-i
]->Vertex
->p
[0][nparam
+1]);
903 if (value_lt(fixed
->p
[1], fixed
->p
[0])) {
905 vertex
[0] = vertex
[1];
910 Vector
*coef
= Vector_Alloc(nparam
+1);
911 for (int i
= 0; i
< 2; ++i
)
912 a
[i
] = affine2evalue(vertex
[i
]->Vertex
->p
[0],
913 vertex
[i
]->Vertex
->p
[0][nparam
+1], nparam
);
914 if (value_one_p(vertex
[0]->Vertex
->p
[0][nparam
+1]))
915 t
[0] = evalue_zero();
917 Vector_Oppose(vertex
[0]->Vertex
->p
[0], coef
->p
, nparam
+1);
918 t
[0] = fractional_part(coef
->p
, vertex
[0]->Vertex
->p
[0][nparam
+1],
921 if (value_one_p(vertex
[1]->Vertex
->p
[0][nparam
+1]))
922 t
[1] = evalue_zero();
924 Vector_Copy(vertex
[1]->Vertex
->p
[0], coef
->p
, nparam
+1);
925 t
[1] = fractional_part(coef
->p
, vertex
[1]->Vertex
->p
[0][nparam
+1],
930 evalue
*I
= evalue_dup(e
);
931 evalue_anti_derive(I
, 0);
932 evalue
*up
= evalue_dup(I
);
934 evalue_substitute(up
, subs_0d
);
937 evalue_substitute(I
, subs_0d
);
944 mu_1d
mu0(degree
, t
[0]);
945 mu_1d
mu1(degree
, t
[1]);
947 evalue
*dx
= evalue_dup(e
);
948 for (int n
= 0; !EVALUE_IS_ZERO(*dx
); ++n
) {
953 evalue_substitute(d
, subs_0d
);
954 emul(mu0
.coefficient(n
), d
);
960 evalue_substitute(d
, subs_0d
);
961 emul(mu1
.coefficient(n
), d
);
967 evalue_derive(dx
, 0);
971 for (int i
= 0; i
< nparam
; ++i
) {
972 free_evalue_refs(subs_0d
[1+i
]);
976 for (int i
= 0; i
< 2; ++i
) {
984 static evalue
*summate_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
985 struct barvinok_options
*options
)
988 Param_Polyhedron
*PP
;
992 struct evalue_section
*s
;
995 MaxRays
= options
->MaxRays
;
996 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
998 U
= Universe_Polyhedron(D
->Dimension
- nvar
);
999 PP
= Polyhedron2Param_Polyhedron(D
, U
, options
);
1001 for (nd
= 0, PD
= PP
->D
; PD
; ++nd
, PD
= PD
->next
);
1002 s
= ALLOCN(struct evalue_section
, nd
);
1004 Polyhedron
*TC
= true_context(D
, U
, MaxRays
);
1005 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, PD
, rVD
)
1008 CA
= align_context(PD
->Domain
, D
->Dimension
, options
->MaxRays
);
1009 F
= DomainIntersection(D
, CA
, options
->MaxRays
);
1012 Vector
*inner
= inner_point(rVD
);
1016 s
[i
].E
= summate_over_1d_pdomain(e
, PP
, PD
, F
, inner
->p
+1, options
);
1017 } else if (nvar
== 2) {
1018 summator_2d
s2d(e
, F
, PP
->nbV
, inner
->p
+1, rVD
->Dimension
);
1020 s
[i
].E
= s2d
.summate_over_pdomain(PP
, PD
, options
);
1025 END_FORALL_REDUCED_DOMAIN
1026 Polyhedron_Free(TC
);
1028 Param_Polyhedron_Free(PP
);
1030 options
->MaxRays
= MaxRays
;
1032 res
= evalue_from_section_array(s
, nd
);
1038 evalue
*euler_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
1046 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
1047 return evalue_dup(e
);
1049 assert(value_zero_p(e
->d
));
1050 assert(e
->x
.p
->type
== partition
);
1052 res
= evalue_zero();
1054 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
1056 t
= summate_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
1057 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), options
);