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
)
147 if (coefficients
[i
]) {
148 free_evalue_refs(coefficients
[i
]);
149 free(coefficients
[i
]);
151 delete [] coefficients
;
153 void compute_coefficient(unsigned n
);
154 const evalue
*coefficient(unsigned n
) {
155 if (!coefficients
[n
])
156 compute_coefficient(n
);
157 return coefficients
[n
];
161 void mu_1d::compute_coefficient(unsigned n
)
163 struct poly_list
*bernoulli
= bernoulli_compute(n
+1);
164 evalue
*c
= evalue_polynomial(bernoulli
->poly
[n
+1], t
);
166 evalue_div(c
, *factorial(n
+1));
172 * Computes the coefficients of
174 * \mu(a)(y) = \sum_{n_1} \sum_{n_2} c_{n_1,n_2} y^{n_1} y^{n_2}
176 * with c_{n1,n2} given by
178 * b(n1+1,t1)/(n1+1)! b(n2+1,t2)/(n2+1)!
179 * - b(n1+n2+2,t2)/(n1+n2+2)! (-c1)^{n1+1}
180 * - b(n1+n2+2,t1)/(n1+n2+2)! (-c2)^{n2+1}
182 * where b(n, t) is the Bernoulli polynomial of degree n in t,
183 * t1(p) and t2(p) are expressions (fractional parts) of the parameters p
184 * such that 0 <= t1(p), t2(p) < 1 for all values of the parameters
185 * and c1 = cn/c1d and c2 = cn/c2d.
186 * The coefficients are computed on demand up to (and including)
187 * the maximal degree (n1,n2) = (max_degree,max_degree).
189 * bernoulli_t[i][j] contains b(j+1, t_i)/(j+1)!
193 evalue
***coefficients
;
194 /* bernoulli_t[i][n] stores b(n+1, t_i)/(n+1)! */
195 evalue
**bernoulli_t
[2];
196 /* stores the powers of -cn */
197 struct power
*power_cn
;
198 struct power
*power_c1d
;
199 struct power
*power_c2d
;
202 mu_2d(unsigned max_degree
, evalue
*t1
, evalue
*t2
,
203 Value cn
, Value c1d
, Value c2d
) : max_degree(max_degree
) {
206 coefficients
= new evalue
**[max_degree
+1];
207 for (int i
= 0; i
< max_degree
+1; ++i
) {
208 coefficients
[i
] = new evalue
*[max_degree
+1];
209 for (int j
= 0; j
< max_degree
+1; ++j
)
210 coefficients
[i
][j
] = NULL
;
212 for (int i
= 0; i
< 2; ++i
) {
213 bernoulli_t
[i
] = new evalue
*[max_degree
+2];
214 for (int j
= 0; j
< max_degree
+2; ++j
)
215 bernoulli_t
[i
][j
] = NULL
;
217 value_oppose(cn
, cn
);
218 power_cn
= new struct power(cn
, max_degree
+1);
219 value_oppose(cn
, cn
);
220 power_c1d
= new struct power(c1d
, max_degree
+1);
221 power_c2d
= new struct power(c2d
, max_degree
+1);
224 for (int i
= 0; i
< max_degree
+1; ++i
) {
225 for (int j
= 0; j
< max_degree
+1; ++j
)
226 if (coefficients
[i
][j
]) {
227 free_evalue_refs(coefficients
[i
][j
]);
228 free(coefficients
[i
][j
]);
230 delete [] coefficients
[i
];
232 delete [] coefficients
;
233 for (int i
= 0; i
< 2; ++i
)
234 for (int j
= 0; j
< max_degree
+2; ++j
)
235 if (bernoulli_t
[i
][j
]) {
236 free_evalue_refs(bernoulli_t
[i
][j
]);
237 free(bernoulli_t
[i
][j
]);
239 for (int i
= 0; i
< 2; ++i
)
240 delete [] bernoulli_t
[i
];
245 const evalue
*get_bernoulli(int n
, int i
);
247 void compute_coefficient(unsigned n1
, unsigned n2
);
248 const evalue
*coefficient(unsigned n1
, unsigned n2
) {
249 if (!coefficients
[n1
][n2
])
250 compute_coefficient(n1
, n2
);
251 return coefficients
[n1
][n2
];
256 * Returns b(n, t_i)/n!
258 const evalue
*mu_2d::get_bernoulli(int n
, int i
)
260 if (!bernoulli_t
[i
][n
-1]) {
261 struct poly_list
*bernoulli
= bernoulli_compute(n
);
262 bernoulli_t
[i
][n
-1] = evalue_polynomial(bernoulli
->poly
[n
], t
[i
]);
263 evalue_div(bernoulli_t
[i
][n
-1], *factorial(n
));
265 return bernoulli_t
[i
][n
-1];
268 void mu_2d::compute_coefficient(unsigned n1
, unsigned n2
)
270 evalue
*c
= evalue_dup(get_bernoulli(n1
+1, 0));
271 emul(get_bernoulli(n2
+1, 1), c
);
273 if (value_notzero_p(*(*power_cn
)[1])) {
277 value_init(neg_power
);
279 t
= evalue_dup(get_bernoulli(n1
+n2
+2, 1));
280 value_multiply(neg_power
,
281 *(*power_cn
)[n1
+1], *binomial(n1
+n2
+1, n1
+1));
282 value_oppose(neg_power
, neg_power
);
283 evalue_mul_div(t
, neg_power
, *(*power_c1d
)[n1
+1]);
288 t
= evalue_dup(get_bernoulli(n1
+n2
+2, 0));
289 value_multiply(neg_power
,
290 *(*power_cn
)[n2
+1], *binomial(n1
+n2
+1, n2
+1));
291 value_oppose(neg_power
, neg_power
);
292 evalue_mul_div(t
, neg_power
, *(*power_c2d
)[n2
+1]);
297 value_clear(neg_power
);
300 coefficients
[n1
][n2
] = c
;
303 /* returns an evalue that corresponds to
307 static evalue
*term(int param
)
309 evalue
*EP
= new evalue();
311 value_set_si(EP
->d
,0);
312 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
313 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
314 evalue_set_si(&EP
->x
.p
->arr
[1], 1, 1);
318 /* Later: avoid recomputation of mu of faces that appear in
319 * more than one domain.
321 struct summator_2d
: public signed_cone_consumer
, public vertex_decomposer
{
322 const evalue
*polynomial
;
326 /* substitutions to use when result is 0-dimensional (only parameters) */
328 /* substitutions to use when result is 1-dimensional */
332 summator_2d(evalue
*e
, Polyhedron
*P
, unsigned nbV
, Value
*inner
,
334 polynomial(e
), vertex_decomposer(P
, nbV
, *this),
335 inner(inner
), nparam(nparam
) {
338 subs_0d
= new evalue
*[2+nparam
];
339 subs_1d
= new evalue
*[2+nparam
];
344 for (int i
= 0; i
< nparam
; ++i
) {
345 subs_0d
[2+i
] = term(i
);
346 subs_1d
[2+i
] = term(1+i
);
350 for (int i
= 0; i
< nparam
; ++i
) {
351 free_evalue_refs(subs_0d
[2+i
]);
353 free_evalue_refs(subs_1d
[2+i
]);
359 evalue
*summate_over_pdomain(Param_Polyhedron
*PP
,
361 struct barvinok_options
*options
);
362 void handle_facet(Param_Polyhedron
*PP
, Param_Domain
*FD
, Value
*normal
);
363 void integrate(Param_Polyhedron
*PP
, Param_Domain
*PD
);
364 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
367 /* Replaces poly by its derivative along variable var */
368 static void evalue_derive(evalue
*poly
, int var
)
370 if (value_notzero_p(poly
->d
)) {
371 value_set_si(poly
->x
.n
, 0);
372 value_set_si(poly
->d
, 1);
375 assert(poly
->x
.p
->type
== polynomial
);
376 if (poly
->x
.p
->pos
-1 > var
) {
377 free_evalue_refs(poly
);
379 evalue_set_si(poly
, 0, 1);
383 if (poly
->x
.p
->pos
-1 < var
) {
384 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
)
385 evalue_derive(&poly
->x
.p
->arr
[i
], var
);
390 assert(poly
->x
.p
->size
> 1);
391 enode
*p
= poly
->x
.p
;
392 free_evalue_refs(&p
->arr
[0]);
394 value_clear(poly
->d
);
403 for (int i
= 0; i
< p
->size
; ++i
) {
404 value_set_si(factor
, i
+1);
405 p
->arr
[i
] = p
->arr
[i
+1];
406 evalue_mul(&p
->arr
[i
], factor
);
411 /* Check whether e is constant with respect to variable var */
412 static int evalue_is_constant(evalue
*e
, int var
)
416 if (value_notzero_p(e
->d
))
418 if (e
->x
.p
->type
== polynomial
&& e
->x
.p
->pos
-1 == var
)
420 assert(e
->x
.p
->type
== polynomial
||
421 e
->x
.p
->type
== fractional
||
422 e
->x
.p
->type
== flooring
);
423 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
424 if (!evalue_is_constant(&e
->x
.p
->arr
[i
], var
))
429 /* Replaces poly by its anti-derivative with constant 0 along variable var */
430 static void evalue_anti_derive(evalue
*poly
, int var
)
434 if (value_zero_p(poly
->d
) &&
435 poly
->x
.p
->type
== polynomial
&& poly
->x
.p
->pos
-1 < var
) {
436 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
)
437 evalue_anti_derive(&poly
->x
.p
->arr
[i
], var
);
442 if (evalue_is_constant(poly
, var
)) {
443 p
= new_enode(polynomial
, 2, 1+var
);
444 evalue_set_si(&p
->arr
[0], 0, 1);
445 value_clear(p
->arr
[1].d
);
451 assert(poly
->x
.p
->type
== polynomial
);
453 p
= new_enode(polynomial
, poly
->x
.p
->size
+1, 1+var
);
454 evalue_set_si(&p
->arr
[0], 0, 1);
455 for (int i
= 0; i
< poly
->x
.p
->size
; ++i
) {
456 value_clear(p
->arr
[1+i
].d
);
457 p
->arr
[1+i
] = poly
->x
.p
->arr
[i
];
458 value_set_si(poly
->d
, 1+i
);
459 evalue_div(&p
->arr
[1+i
], poly
->d
);
463 value_set_si(poly
->d
, 0);
466 /* Computes offsets in the basis given by the rays of the cone
467 * to the integer point in the fundamental parallelepiped of
469 * The resulting evalues contain only the parameters as variables.
471 evalue
**offsets_to_integer_point(Matrix
*Rays
, Matrix
*vertex
)
473 unsigned nparam
= vertex
->NbColumns
-2;
474 evalue
**t
= new evalue
*[2];
476 if (value_one_p(vertex
->p
[0][nparam
+1])) {
477 t
[0] = evalue_zero();
478 t
[1] = evalue_zero();
480 Matrix
*R2
= Matrix_Copy(Rays
);
481 Matrix_Transposition(R2
);
482 Matrix
*inv
= Matrix_Alloc(Rays
->NbColumns
, Rays
->NbRows
);
483 int ok
= Matrix_Inverse(R2
, inv
);
487 /* We want the fractional part of the negative relative coordinates */
488 Vector_Oppose(inv
->p
[0], inv
->p
[0], inv
->NbColumns
);
489 Vector_Oppose(inv
->p
[1], inv
->p
[1], inv
->NbColumns
);
491 Matrix
*neg_rel
= Matrix_Alloc(2, nparam
+1);
493 Matrix_Product(inv
, vertex
, neg_rel
);
496 t
[0] = fractional_part(neg_rel
->p
[0], vertex
->p
[0][nparam
+1],
498 t
[1] = fractional_part(neg_rel
->p
[1], vertex
->p
[0][nparam
+1],
500 Matrix_Free(neg_rel
);
507 * Called from decompose_at_vertex.
509 * Handles a cone in the signed decomposition of the supporting
510 * cone of a vertex. The cone is assumed to be unimodular.
512 void summator_2d::handle(const signed_cone
& sc
, barvinok_options
*options
)
515 unsigned degree
= total_degree(polynomial
, 2);
517 subs_0d
[0] = affine2evalue(V
->Vertex
->p
[0],
518 V
->Vertex
->p
[0][nparam
+1], nparam
);
519 subs_0d
[1] = affine2evalue(V
->Vertex
->p
[1],
520 V
->Vertex
->p
[1][nparam
+1], nparam
);
524 assert(V
->Vertex
->NbRows
> 0);
525 Param_Vertex_Common_Denominator(V
);
527 Matrix
*Rays
= zz2matrix(sc
.rays
);
529 t
= offsets_to_integer_point(Rays
, V
->Vertex
);
531 Vector
*c
= Vector_Alloc(3);
532 Inner_Product(Rays
->p
[0], Rays
->p
[1], 2, &c
->p
[0]);
533 Inner_Product(Rays
->p
[0], Rays
->p
[0], 2, &c
->p
[1]);
534 Inner_Product(Rays
->p
[1], Rays
->p
[1], 2, &c
->p
[2]);
536 mu_2d
mu(degree
, t
[0], t
[1], c
->p
[0], c
->p
[1], c
->p
[2]);
539 struct power
power_r00(Rays
->p
[0][0], degree
);
540 struct power
power_r01(Rays
->p
[0][1], degree
);
541 struct power
power_r10(Rays
->p
[1][0], degree
);
542 struct power
power_r11(Rays
->p
[1][1], degree
);
544 Value factor
, tmp1
, tmp2
;
548 evalue
*res
= evalue_zero();
549 evalue
*dx1
= evalue_dup(polynomial
);
550 for (int i
= 0; !EVALUE_IS_ZERO(*dx1
); ++i
) {
551 evalue
*dx2
= evalue_dup(dx1
);
552 for (int j
= 0; !EVALUE_IS_ZERO(*dx2
); ++j
) {
553 evalue
*cij
= evalue_zero();
554 for (int n1
= 0; n1
<= i
+j
; ++n1
) {
556 value_set_si(factor
, 0);
557 for (int k
= max(0, i
-n2
); k
<= i
&& k
<= n1
; ++k
) {
559 value_multiply(tmp1
, *power_r00
[k
], *power_r01
[n1
-k
]);
560 value_multiply(tmp1
, tmp1
, *power_r10
[l
]);
561 value_multiply(tmp1
, tmp1
, *power_r11
[n2
-l
]);
562 value_multiply(tmp1
, tmp1
, *binomial(n1
, k
));
563 value_multiply(tmp1
, tmp1
, *binomial(n2
, l
));
564 value_addto(factor
, factor
, tmp1
);
566 if (value_zero_p(factor
))
569 evalue
*c
= evalue_dup(mu
.coefficient(n1
, n2
));
570 evalue_mul(c
, factor
);
575 evalue
*d
= evalue_dup(dx2
);
576 evalue_substitute(d
, subs_0d
);
578 free_evalue_refs(cij
);
583 evalue_derive(dx2
, 1);
585 free_evalue_refs(dx2
);
587 evalue_derive(dx1
, 0);
589 free_evalue_refs(dx1
);
591 for (int i
= 0; i
< 2; ++i
) {
592 free_evalue_refs(subs_0d
[i
]);
595 free_evalue_refs(t
[i
]);
607 free_evalue_refs(res
);
611 evalue
*summator_2d::summate_over_pdomain(Param_Polyhedron
*PP
,
613 struct barvinok_options
*options
)
618 assert(PP
->V
->Vertex
->NbRows
== 2);
620 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, PP
) // _i internal counter
621 decompose_at_vertex(V
, _i
, options
);
622 END_FORALL_PVertex_in_ParamPolyhedron
;
624 Vector
*normal
= Vector_Alloc(2);
625 for (j
= P
->NbEq
; j
< P
->NbConstraints
; ++j
) {
627 Vector_Copy(P
->Constraint
[j
]+1, normal
->p
, 2);
628 if (value_zero_p(normal
->p
[0]) && value_zero_p(normal
->p
[1]))
631 FD
= Param_Polyhedron_Facet(PP
, PD
, P
, j
);
632 Vector_Normalize(normal
->p
, 2);
633 handle_facet(PP
, FD
, normal
->p
);
634 Param_Domain_Free(FD
);
643 void summator_2d::handle_facet(Param_Polyhedron
*PP
, Param_Domain
*FD
,
648 Param_Vertices
*vertex
[2];
650 unsigned degree
= total_degree(polynomial
, 2);
652 FORALL_PVertex_in_ParamPolyhedron(V
, FD
, PP
)
654 END_FORALL_PVertex_in_ParamPolyhedron
;
658 /* We can take either vertex[0] or vertex[1];
659 * the result should be the same
661 Param_Vertex_Common_Denominator(vertex
[0]);
663 /* The extra variable in front is the coordinate along the facet. */
664 Vector
*coef_normal
= Vector_Alloc(1 + nparam
+ 2);
665 Vector_Combine(vertex
[0]->Vertex
->p
[0], vertex
[0]->Vertex
->p
[1],
666 coef_normal
->p
+1, normal
[0], normal
[1], nparam
+1);
667 value_assign(coef_normal
->p
[1+nparam
+1], vertex
[0]->Vertex
->p
[0][nparam
+1]);
668 Vector_Normalize(coef_normal
->p
, coef_normal
->Size
);
670 Vector
*base
= Vector_Alloc(2);
671 value_oppose(base
->p
[0], normal
[1]);
672 value_assign(base
->p
[1], normal
[0]);
675 Inner_Product(normal
, normal
, 2, &det
);
677 Vector
*s
= Vector_Alloc(1+nparam
+2);
678 value_multiply(s
->p
[1+nparam
+1], coef_normal
->p
[1+nparam
+1], det
);
679 value_multiply(s
->p
[0], base
->p
[0], s
->p
[1+nparam
+1]);
680 Vector_Scale(coef_normal
->p
+1, s
->p
+1, normal
[0], nparam
+1);
681 subs_1d
[0] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
682 value_multiply(s
->p
[0], base
->p
[1], s
->p
[1+nparam
+1]);
683 Vector_Scale(coef_normal
->p
+1, s
->p
+1, normal
[1], nparam
+1);
684 subs_1d
[1] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
688 if (value_one_p(coef_normal
->p
[coef_normal
->Size
-1]))
691 Vector_Oppose(coef_normal
->p
+1, coef_normal
->p
+1, nparam
+1);
692 t
= fractional_part(coef_normal
->p
,
693 coef_normal
->p
[coef_normal
->Size
-1],
696 Vector_Free(coef_normal
);
700 struct power
power_normal0(normal
[0], degree
);
701 struct power
power_normal1(normal
[1], degree
);
702 struct power
power_det(det
, degree
);
706 evalue
*res
= evalue_zero();
707 evalue
*dx1
= evalue_dup(polynomial
);
708 for (int i
= 0; !EVALUE_IS_ZERO(*dx1
); ++i
) {
709 evalue
*dx2
= evalue_dup(dx1
);
710 for (int j
= 0; !EVALUE_IS_ZERO(*dx2
); ++j
) {
711 value_multiply(factor
, *power_normal0
[i
], *power_normal1
[j
]);
712 if (value_notzero_p(factor
)) {
713 value_multiply(factor
, factor
, *binomial(i
+j
, i
));
715 evalue
*c
= evalue_dup(mu
.coefficient(i
+j
));
716 evalue_mul_div(c
, factor
, *power_det
[i
+j
]);
718 evalue
*d
= evalue_dup(dx2
);
719 evalue_substitute(d
, subs_1d
);
727 evalue_derive(dx2
, 1);
729 free_evalue_refs(dx2
);
731 evalue_derive(dx1
, 0);
733 free_evalue_refs(dx1
);
737 for (int i
= 0; i
< 2; ++i
) {
738 free_evalue_refs(subs_1d
[i
]);
743 evalue_anti_derive(res
, 0);
745 Matrix
*z
= Matrix_Alloc(2, nparam
+2);
746 Vector
*fixed_z
= Vector_Alloc(2);
747 for (int i
= 0; i
< 2; ++i
) {
748 Vector_Combine(vertex
[i
]->Vertex
->p
[0], vertex
[i
]->Vertex
->p
[1],
749 z
->p
[i
], base
->p
[0], base
->p
[1], nparam
+1);
750 value_multiply(z
->p
[i
][nparam
+1],
751 det
, vertex
[i
]->Vertex
->p
[0][nparam
+1]);
752 Inner_Product(z
->p
[i
], inner
, nparam
+1, &fixed_z
->p
[i
]);
755 /* Put on a common denominator */
756 value_multiply(fixed_z
->p
[0], fixed_z
->p
[0], z
->p
[1][nparam
+1]);
757 value_multiply(fixed_z
->p
[1], fixed_z
->p
[1], z
->p
[0][nparam
+1]);
758 /* Make sure z->p[0] is smaller than z->p[1] (for an internal
759 * point of the chamber and hence for all parameter values in
760 * the chamber), to ensure we integrate in the right direction.
762 if (value_lt(fixed_z
->p
[1], fixed_z
->p
[0]))
763 Vector_Exchange(z
->p
[0], z
->p
[1], nparam
+2);
764 Vector_Free(fixed_z
);
767 subs_0d
[1] = affine2evalue(z
->p
[1], z
->p
[1][nparam
+1], nparam
);
768 evalue
*up
= evalue_dup(res
);
769 evalue_substitute(up
, subs_0d
+1);
770 free_evalue_refs(subs_0d
[1]);
773 subs_0d
[1] = affine2evalue(z
->p
[0], z
->p
[0][nparam
+1], nparam
);
774 evalue_substitute(res
, subs_0d
+1);
777 free_evalue_refs(subs_0d
[1]);
780 free_evalue_refs(up
);
786 free_evalue_refs(res
);
790 /* Integrate the polynomial over the whole polygon using
791 * the Green-Stokes theorem.
793 void summator_2d::integrate(Param_Polyhedron
*PP
, Param_Domain
*PD
)
796 evalue
*res
= evalue_zero();
798 evalue
*I
= evalue_dup(polynomial
);
799 evalue_anti_derive(I
, 0);
802 Vector
*normal
= Vector_Alloc(2);
803 Vector
*dir
= Vector_Alloc(2);
804 Matrix
*v0v1
= Matrix_Alloc(2, nparam
+2);
805 Vector
*f_v0v1
= Vector_Alloc(2);
806 Vector
*s
= Vector_Alloc(1+nparam
+2);
807 for (int j
= P
->NbEq
; j
< P
->NbConstraints
; ++j
) {
810 Param_Vertices
*vertex
[2];
811 Vector_Copy(P
->Constraint
[j
]+1, normal
->p
, 2);
813 if (value_zero_p(normal
->p
[0]))
816 Vector_Normalize(normal
->p
, 2);
817 value_assign(dir
->p
[0], normal
->p
[1]);
818 value_oppose(dir
->p
[1], normal
->p
[0]);
820 FD
= Param_Polyhedron_Facet(PP
, PD
, P
, j
);
822 FORALL_PVertex_in_ParamPolyhedron(V
, FD
, PP
)
824 END_FORALL_PVertex_in_ParamPolyhedron
;
828 Param_Vertex_Common_Denominator(vertex
[0]);
829 Param_Vertex_Common_Denominator(vertex
[1]);
831 value_oppose(tmp
, vertex
[1]->Vertex
->p
[0][nparam
+1]);
832 for (int i
= 0; i
< 2; ++i
)
833 Vector_Combine(vertex
[1]->Vertex
->p
[i
],
834 vertex
[0]->Vertex
->p
[i
],
836 vertex
[0]->Vertex
->p
[0][nparam
+1], tmp
, nparam
+1);
837 value_multiply(v0v1
->p
[0][nparam
+1],
838 vertex
[0]->Vertex
->p
[0][nparam
+1],
839 vertex
[1]->Vertex
->p
[0][nparam
+1]);
840 value_assign(v0v1
->p
[1][nparam
+1], v0v1
->p
[0][nparam
+1]);
842 /* Order vertices to ensure we integrate in the right
843 * direction, i.e., with the polytope "on the left".
845 for (int i
= 0; i
< 2; ++i
)
846 Inner_Product(v0v1
->p
[i
], inner
, nparam
+1, &f_v0v1
->p
[i
]);
848 Inner_Product(dir
->p
, f_v0v1
->p
, 2, &tmp
);
849 if (value_neg_p(tmp
)) {
850 Param_Vertices
*PV
= vertex
[0];
851 vertex
[0] = vertex
[1];
853 for (int i
= 0; i
< 2; ++i
)
854 Vector_Oppose(v0v1
->p
[i
], v0v1
->p
[i
], nparam
+1);
856 value_oppose(tmp
, normal
->p
[0]);
857 if (value_neg_p(tmp
)) {
858 value_oppose(tmp
, tmp
);
859 Vector_Oppose(v0v1
->p
[1], v0v1
->p
[1], nparam
+1);
861 value_multiply(tmp
, tmp
, v0v1
->p
[1][nparam
+1]);
862 evalue
*top
= affine2evalue(v0v1
->p
[1], tmp
, nparam
);
864 value_multiply(s
->p
[0], normal
->p
[1], vertex
[0]->Vertex
->p
[0][nparam
+1]);
865 Vector_Copy(vertex
[0]->Vertex
->p
[0], s
->p
+1, nparam
+2);
866 subs_1d
[0] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
867 value_multiply(s
->p
[0], normal
->p
[0], vertex
[0]->Vertex
->p
[0][nparam
+1]);
868 value_oppose(s
->p
[0], s
->p
[0]);
869 Vector_Copy(vertex
[0]->Vertex
->p
[1], s
->p
+1, nparam
+2);
870 subs_1d
[1] = affine2evalue(s
->p
, s
->p
[1+nparam
+1], 1+nparam
);
872 evalue
*d
= evalue_dup(I
);
873 evalue_substitute(d
, subs_1d
);
874 evalue_anti_derive(d
, 0);
876 free_evalue_refs(subs_1d
[0]);
878 free_evalue_refs(subs_1d
[1]);
882 evalue_substitute(d
, subs_0d
+1);
883 evalue_mul(d
, dir
->p
[1]);
884 free_evalue_refs(subs_0d
[1]);
891 Param_Domain_Free(FD
);
903 free_evalue_refs(res
);
907 evalue
*summate_over_1d_pdomain(evalue
*e
,
908 Param_Polyhedron
*PP
, Param_Domain
*PD
,
909 Polyhedron
*P
, Value
*inner
,
910 struct barvinok_options
*options
)
914 Param_Vertices
*vertex
[2];
915 unsigned nparam
= PP
->V
->Vertex
->NbColumns
-2;
916 evalue
*subs_0d
[1+nparam
];
919 unsigned degree
= total_degree(e
, 1);
921 for (int i
= 0; i
< nparam
; ++i
)
922 subs_0d
[1+i
] = term(i
);
924 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, PP
)
926 END_FORALL_PVertex_in_ParamPolyhedron
;
929 Vector
*fixed
= Vector_Alloc(2);
930 for (int i
= 0; i
< 2; ++i
) {
931 Inner_Product(vertex
[i
]->Vertex
->p
[0], inner
, nparam
+1, &fixed
->p
[i
]);
932 value_multiply(fixed
->p
[i
],
933 fixed
->p
[i
], vertex
[1-i
]->Vertex
->p
[0][nparam
+1]);
935 if (value_lt(fixed
->p
[1], fixed
->p
[0])) {
937 vertex
[0] = vertex
[1];
942 Vector
*coef
= Vector_Alloc(nparam
+1);
943 for (int i
= 0; i
< 2; ++i
)
944 a
[i
] = affine2evalue(vertex
[i
]->Vertex
->p
[0],
945 vertex
[i
]->Vertex
->p
[0][nparam
+1], nparam
);
946 if (value_one_p(vertex
[0]->Vertex
->p
[0][nparam
+1]))
947 t
[0] = evalue_zero();
949 Vector_Oppose(vertex
[0]->Vertex
->p
[0], coef
->p
, nparam
+1);
950 t
[0] = fractional_part(coef
->p
, vertex
[0]->Vertex
->p
[0][nparam
+1],
953 if (value_one_p(vertex
[1]->Vertex
->p
[0][nparam
+1]))
954 t
[1] = evalue_zero();
956 Vector_Copy(vertex
[1]->Vertex
->p
[0], coef
->p
, nparam
+1);
957 t
[1] = fractional_part(coef
->p
, vertex
[1]->Vertex
->p
[0][nparam
+1],
962 evalue
*I
= evalue_dup(e
);
963 evalue_anti_derive(I
, 0);
964 evalue
*up
= evalue_dup(I
);
966 evalue_substitute(up
, subs_0d
);
969 evalue_substitute(I
, subs_0d
);
972 free_evalue_refs(up
);
977 mu_1d
mu0(degree
, t
[0]);
978 mu_1d
mu1(degree
, t
[1]);
980 evalue
*dx
= evalue_dup(e
);
981 for (int n
= 0; !EVALUE_IS_ZERO(*dx
); ++n
) {
986 evalue_substitute(d
, subs_0d
);
987 emul(mu0
.coefficient(n
), d
);
994 evalue_substitute(d
, subs_0d
);
995 emul(mu1
.coefficient(n
), d
);
1002 evalue_derive(dx
, 0);
1004 free_evalue_refs(dx
);
1007 for (int i
= 0; i
< nparam
; ++i
) {
1008 free_evalue_refs(subs_0d
[1+i
]);
1009 delete subs_0d
[1+i
];
1012 for (int i
= 0; i
< 2; ++i
) {
1013 free_evalue_refs(a
[i
]);
1015 free_evalue_refs(t
[i
]);
1022 static evalue
*summate_over_domain(evalue
*e
, int nvar
, Polyhedron
*D
,
1023 struct barvinok_options
*options
)
1026 Param_Polyhedron
*PP
;
1030 struct evalue_section
*s
;
1033 MaxRays
= options
->MaxRays
;
1034 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
1036 U
= Universe_Polyhedron(D
->Dimension
- nvar
);
1037 PP
= Polyhedron2Param_Polyhedron(D
, U
, options
);
1039 for (nd
= 0, PD
= PP
->D
; PD
; ++nd
, PD
= PD
->next
);
1040 s
= ALLOCN(struct evalue_section
, nd
);
1042 Polyhedron
*TC
= true_context(D
, U
, MaxRays
);
1043 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, PD
, rVD
)
1046 CA
= align_context(PD
->Domain
, D
->Dimension
, options
->MaxRays
);
1047 F
= DomainIntersection(D
, CA
, options
->MaxRays
);
1050 Vector
*inner
= inner_point(rVD
);
1054 s
[i
].E
= summate_over_1d_pdomain(e
, PP
, PD
, F
, inner
->p
+1, options
);
1055 } else if (nvar
== 2) {
1056 summator_2d
s2d(e
, F
, PP
->nbV
, inner
->p
+1, rVD
->Dimension
);
1058 s
[i
].E
= s2d
.summate_over_pdomain(PP
, PD
, options
);
1063 END_FORALL_REDUCED_DOMAIN
1064 Polyhedron_Free(TC
);
1066 Param_Polyhedron_Free(PP
);
1068 options
->MaxRays
= MaxRays
;
1070 res
= evalue_from_section_array(s
, nd
);
1076 evalue
*euler_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
1084 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
1085 return evalue_dup(e
);
1087 assert(value_zero_p(e
->d
));
1088 assert(e
->x
.p
->type
== partition
);
1090 res
= evalue_zero();
1092 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
1094 t
= summate_over_domain(&e
->x
.p
->arr
[2*i
+1], nvar
,
1095 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), options
);
1097 free_evalue_refs(t
);