2 #include <barvinok/options.h>
3 #include <barvinok/util.h>
6 #define ALLOC(type) (type*)malloc(sizeof(type))
7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
8 #define REALLOCN(ptr,type,n) (type*)realloc(ptr, (n) * sizeof(type))
10 static struct bernoulli_coef bernoulli_coef
;
11 static struct poly_list bernoulli
;
12 static struct poly_list faulhaber
;
14 struct bernoulli_coef
*bernoulli_coef_compute(int n
)
19 if (n
< bernoulli_coef
.n
)
20 return &bernoulli_coef
;
22 if (n
>= bernoulli_coef
.size
) {
23 int size
= 3*(n
+ 5)/2;
26 b
= Vector_Alloc(size
);
27 if (bernoulli_coef
.n
) {
28 Vector_Copy(bernoulli_coef
.num
->p
, b
->p
, bernoulli_coef
.n
);
29 Vector_Free(bernoulli_coef
.num
);
31 bernoulli_coef
.num
= b
;
32 b
= Vector_Alloc(size
);
33 if (bernoulli_coef
.n
) {
34 Vector_Copy(bernoulli_coef
.den
->p
, b
->p
, bernoulli_coef
.n
);
35 Vector_Free(bernoulli_coef
.den
);
37 bernoulli_coef
.den
= b
;
38 b
= Vector_Alloc(size
);
39 if (bernoulli_coef
.n
) {
40 Vector_Copy(bernoulli_coef
.lcm
->p
, b
->p
, bernoulli_coef
.n
);
41 Vector_Free(bernoulli_coef
.lcm
);
43 bernoulli_coef
.lcm
= b
;
45 bernoulli_coef
.size
= size
;
49 for (i
= bernoulli_coef
.n
; i
<= n
; ++i
) {
51 value_set_si(bernoulli_coef
.num
->p
[0], 1);
52 value_set_si(bernoulli_coef
.den
->p
[0], 1);
53 value_set_si(bernoulli_coef
.lcm
->p
[0], 1);
56 value_set_si(bernoulli_coef
.num
->p
[i
], 0);
57 value_set_si(factor
, -(i
+1));
58 for (j
= i
-1; j
>= 0; --j
) {
59 mpz_mul_ui(factor
, factor
, j
+1);
60 mpz_divexact_ui(factor
, factor
, i
+1-j
);
61 value_division(tmp
, bernoulli_coef
.lcm
->p
[i
-1],
62 bernoulli_coef
.den
->p
[j
]);
63 value_multiply(tmp
, tmp
, bernoulli_coef
.num
->p
[j
]);
64 value_multiply(tmp
, tmp
, factor
);
65 value_addto(bernoulli_coef
.num
->p
[i
], bernoulli_coef
.num
->p
[i
], tmp
);
67 mpz_mul_ui(bernoulli_coef
.den
->p
[i
], bernoulli_coef
.lcm
->p
[i
-1], i
+1);
68 value_gcd(tmp
, bernoulli_coef
.num
->p
[i
], bernoulli_coef
.den
->p
[i
]);
69 if (value_notone_p(tmp
)) {
70 value_division(bernoulli_coef
.num
->p
[i
],
71 bernoulli_coef
.num
->p
[i
], tmp
);
72 value_division(bernoulli_coef
.den
->p
[i
],
73 bernoulli_coef
.den
->p
[i
], tmp
);
75 value_lcm(bernoulli_coef
.lcm
->p
[i
],
76 bernoulli_coef
.lcm
->p
[i
-1], bernoulli_coef
.den
->p
[i
]);
78 bernoulli_coef
.n
= n
+1;
82 return &bernoulli_coef
;
86 * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
88 * B_n = sum_{k=0}^n { n \choose k } b_k x^{n-k}
89 * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
91 static struct poly_list
*bernoulli_faulhaber_compute(int n
, struct poly_list
*pl
,
96 struct bernoulli_coef
*bc
;
102 int size
= 3*(n
+ 5)/2;
105 poly
= ALLOCN(Vector
*, size
);
106 for (i
= 0; i
< pl
->n
; ++i
)
107 poly
[i
] = pl
->poly
[i
];
114 bc
= bernoulli_coef_compute(n
);
117 for (i
= pl
->n
; i
<= n
; ++i
) {
118 pl
->poly
[i
] = Vector_Alloc(i
+faulhaber
+2);
119 value_assign(pl
->poly
[i
]->p
[i
+faulhaber
], bc
->lcm
->p
[i
]);
121 mpz_mul_ui(pl
->poly
[i
]->p
[i
+2], bc
->lcm
->p
[i
], i
+1);
123 value_assign(pl
->poly
[i
]->p
[i
+1], bc
->lcm
->p
[i
]);
124 value_set_si(factor
, 1);
125 for (j
= 1; j
<= i
; ++j
) {
126 mpz_mul_ui(factor
, factor
, i
+faulhaber
+1-j
);
127 mpz_divexact_ui(factor
, factor
, j
);
128 value_division(pl
->poly
[i
]->p
[i
+faulhaber
-j
],
129 bc
->lcm
->p
[i
], bc
->den
->p
[j
]);
130 value_multiply(pl
->poly
[i
]->p
[i
+faulhaber
-j
],
131 pl
->poly
[i
]->p
[i
+faulhaber
-j
], bc
->num
->p
[j
]);
132 value_multiply(pl
->poly
[i
]->p
[i
+faulhaber
-j
],
133 pl
->poly
[i
]->p
[i
+faulhaber
-j
], factor
);
135 Vector_Normalize(pl
->poly
[i
]->p
, i
+faulhaber
+2);
143 struct poly_list
*bernoulli_compute(int n
)
145 return bernoulli_faulhaber_compute(n
, &bernoulli
, 0);
148 struct poly_list
*faulhaber_compute(int n
)
150 return bernoulli_faulhaber_compute(n
, &faulhaber
, 1);
153 /* shift variables in polynomial one down */
154 static void shift(evalue
*e
)
157 if (value_notzero_p(e
->d
))
159 assert(e
->x
.p
->type
== polynomial
);
160 assert(e
->x
.p
->pos
> 1);
162 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
163 shift(&e
->x
.p
->arr
[i
]);
166 /* shift variables in polynomial n up */
167 static void unshift(evalue
*e
, unsigned n
)
170 if (value_notzero_p(e
->d
))
172 assert(e
->x
.p
->type
== polynomial
|| e
->x
.p
->type
== fractional
);
173 if (e
->x
.p
->type
== polynomial
)
175 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
176 unshift(&e
->x
.p
->arr
[i
], n
);
179 static evalue
*shifted_copy(evalue
*src
)
181 evalue
*e
= ALLOC(evalue
);
188 static evalue
*power_sums(struct poly_list
*faulhaber
, evalue
*poly
,
189 Vector
*arg
, Value denom
, int neg
, int alt_neg
)
192 evalue
*base
= affine2evalue(arg
->p
, denom
, arg
->Size
-1);
193 evalue
*sum
= evalue_zero();
195 for (i
= 1; i
< poly
->x
.p
->size
; ++i
) {
196 evalue
*term
= evalue_polynomial(faulhaber
->poly
[i
], base
);
197 evalue
*factor
= shifted_copy(&poly
->x
.p
->arr
[i
]);
199 if (alt_neg
&& (i
% 2))
212 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
213 * +- (b y + c) +- a >=,> 0
216 * sign_affine sign_cst
218 static void bound_constraint(Value
*c
, unsigned dim
,
220 int sign_affine
, int sign_cst
, int strict
)
222 if (sign_affine
== -1)
223 Vector_Oppose(cst_affine
+1, c
, dim
+1);
225 Vector_Copy(cst_affine
+1, c
, dim
+1);
228 value_subtract(c
[dim
], c
[dim
], cst_affine
[0]);
229 else if (sign_cst
== 1)
230 value_addto(c
[dim
], c
[dim
], cst_affine
[0]);
233 value_decrement(c
[dim
], c
[dim
]);
236 struct Bernoulli_data
{
238 struct evalue_section
*s
;
244 static evalue
*compute_poly_u(evalue
*poly_u
, Value
*upper
, Vector
*row
,
245 unsigned dim
, Value tmp
,
246 struct poly_list
*faulhaber
,
247 struct Bernoulli_data
*data
)
251 Vector_Copy(upper
+2, row
->p
, dim
+1);
252 value_oppose(tmp
, upper
[1]);
253 value_addto(row
->p
[dim
], row
->p
[dim
], tmp
);
254 return power_sums(faulhaber
, data
->e
, row
, tmp
, 0, 0);
257 static evalue
*compute_poly_l(evalue
*poly_l
, Value
*lower
, Vector
*row
,
259 struct poly_list
*faulhaber
,
260 struct Bernoulli_data
*data
)
264 Vector_Copy(lower
+2, row
->p
, dim
+1);
265 value_addto(row
->p
[dim
], row
->p
[dim
], lower
[1]);
266 return power_sums(faulhaber
, data
->e
, row
, lower
[1], 0, 1);
269 static void Bernoulli_init(unsigned n
, void *cb_data
)
271 struct Bernoulli_data
*data
= (struct Bernoulli_data
*)cb_data
;
274 if (cases
* n
<= data
->size
)
277 data
->size
= cases
* (n
+ 16);
278 data
->s
= REALLOCN(data
->s
, struct evalue_section
, data
->size
);
281 static void Bernoulli_cb(Matrix
*M
, Value
*lower
, Value
*upper
, void *cb_data
)
283 struct Bernoulli_data
*data
= (struct Bernoulli_data
*)cb_data
;
286 evalue
*factor
= NULL
;
287 evalue
*linear
= NULL
;
290 unsigned dim
= M
->NbColumns
-2;
296 assert(data
->ns
+ cases
<= data
->size
);
299 T
= Constraints2Polyhedron(M2
, data
->MaxRays
);
302 POL_ENSURE_VERTICES(T
);
308 assert(lower
!= upper
);
310 row
= Vector_Alloc(dim
+1);
312 if (value_notzero_p(data
->e
->d
)) {
316 assert(data
->e
->x
.p
->type
== polynomial
);
317 if (data
->e
->x
.p
->pos
> 1) {
318 factor
= shifted_copy(data
->e
);
321 factor
= shifted_copy(&data
->e
->x
.p
->arr
[0]);
323 if (!EVALUE_IS_ZERO(*factor
)) {
324 value_absolute(tmp
, upper
[1]);
326 Vector_Combine(lower
+2, upper
+2, row
->p
, tmp
, lower
[1], dim
+1);
327 value_multiply(tmp
, tmp
, lower
[1]);
328 /* upper - lower + 1 */
329 value_addto(row
->p
[dim
], row
->p
[dim
], tmp
);
331 linear
= affine2evalue(row
->p
, tmp
, dim
);
332 emul(factor
, linear
);
334 linear
= evalue_zero();
337 data
->s
[data
->ns
].E
= linear
;
338 data
->s
[data
->ns
].D
= T
;
341 evalue
*poly_u
= NULL
, *poly_l
= NULL
;
343 struct poly_list
*faulhaber
;
344 assert(data
->e
->x
.p
->type
== polynomial
);
345 assert(data
->e
->x
.p
->pos
== 1);
346 faulhaber
= faulhaber_compute(data
->e
->x
.p
->size
-1);
347 /* lower is the constraint
348 * b i - l' >= 0 i >= l'/b = l
349 * upper is the constraint
350 * -a i + u' >= 0 i <= u'/a = u
352 M2
= Matrix_Alloc(3, 2+T
->Dimension
);
353 value_set_si(M2
->p
[0][0], 1);
354 value_set_si(M2
->p
[1][0], 1);
355 value_set_si(M2
->p
[2][0], 1);
359 bound_constraint(M2
->p
[0]+1, T
->Dimension
, lower
+1, -1, -1, 0);
360 D
= AddConstraints(M2
->p_Init
, 1, T
, data
->MaxRays
);
361 POL_ENSURE_VERTICES(D
);
366 poly_u
= compute_poly_u(poly_u
, upper
, row
, dim
, tmp
,
368 Vector_Oppose(lower
+2, row
->p
, dim
+1);
369 extra
= power_sums(faulhaber
, data
->e
, row
, lower
[1], 1, 0);
373 data
->s
[data
->ns
].E
= extra
;
374 data
->s
[data
->ns
].D
= D
;
379 * 1 <= -u -u' - a >= 0
381 bound_constraint(M2
->p
[0]+1, T
->Dimension
, upper
+1, -1, 1, 0);
382 D
= AddConstraints(M2
->p_Init
, 1, T
, data
->MaxRays
);
383 POL_ENSURE_VERTICES(D
);
388 poly_l
= compute_poly_l(poly_l
, lower
, row
, dim
, faulhaber
, data
);
389 Vector_Oppose(upper
+2, row
->p
, dim
+1);
390 value_oppose(tmp
, upper
[1]);
391 extra
= power_sums(faulhaber
, data
->e
, row
, tmp
, 1, 1);
395 data
->s
[data
->ns
].E
= extra
;
396 data
->s
[data
->ns
].D
= D
;
404 bound_constraint(M2
->p
[0]+1, T
->Dimension
, upper
+1, 1, 0, 0);
405 bound_constraint(M2
->p
[1]+1, T
->Dimension
, lower
+1, 1, 0, 0);
406 D
= AddConstraints(M2
->p_Init
, 2, T
, data
->MaxRays
);
407 POL_ENSURE_VERTICES(D
);
411 poly_l
= compute_poly_l(poly_l
, lower
, row
, dim
, faulhaber
, data
);
412 poly_u
= compute_poly_u(poly_u
, upper
, row
, dim
, tmp
,
415 data
->s
[data
->ns
].E
= ALLOC(evalue
);
416 value_init(data
->s
[data
->ns
].E
->d
);
417 evalue_copy(data
->s
[data
->ns
].E
, poly_u
);
418 eadd(poly_l
, data
->s
[data
->ns
].E
);
419 eadd(linear
, data
->s
[data
->ns
].E
);
420 data
->s
[data
->ns
].D
= D
;
425 * l < 1 -l' + b - 1 >= 0
429 bound_constraint(M2
->p
[0]+1, T
->Dimension
, lower
+1, 1, 1, 1);
430 bound_constraint(M2
->p
[1]+1, T
->Dimension
, lower
+1, -1, 0, 1);
431 bound_constraint(M2
->p
[2]+1, T
->Dimension
, upper
+1, 1, 1, 0);
432 D
= AddConstraints(M2
->p_Init
, 3, T
, data
->MaxRays
);
433 POL_ENSURE_VERTICES(D
);
437 poly_u
= compute_poly_u(poly_u
, upper
, row
, dim
, tmp
,
440 eadd(linear
, poly_u
);
441 data
->s
[data
->ns
].E
= poly_u
;
443 data
->s
[data
->ns
].D
= D
;
448 * -u < 1 u' + a - 1 >= 0
449 * 0 < -u -u' - 1 >= 0
450 * l <= -1 -l' - b >= 0
452 bound_constraint(M2
->p
[0]+1, T
->Dimension
, upper
+1, 1, -1, 1);
453 bound_constraint(M2
->p
[1]+1, T
->Dimension
, upper
+1, -1, 0, 1);
454 bound_constraint(M2
->p
[2]+1, T
->Dimension
, lower
+1, 1, -1, 0);
455 D
= AddConstraints(M2
->p_Init
, 3, T
, data
->MaxRays
);
456 POL_ENSURE_VERTICES(D
);
460 poly_l
= compute_poly_l(poly_l
, lower
, row
, dim
, faulhaber
, data
);
462 eadd(linear
, poly_l
);
463 data
->s
[data
->ns
].E
= poly_l
;
465 data
->s
[data
->ns
].D
= D
;
477 if (factor
!= data
->e
)
483 /* Looks for variable with integer bounds, i.e., with coefficients 0, 1 or -1.
484 * Returns 1 if such a variable is found and puts it in the first position,
485 * possibly changing *P_p and *E_p.
487 static int find_integer_bounds(Polyhedron
**P_p
, evalue
**E_p
, unsigned nvar
)
489 Polyhedron
*P
= *P_p
;
491 unsigned dim
= P
->Dimension
;
494 for (i
= 0; i
< nvar
; ++i
) {
495 for (j
= 0; j
< P
->NbConstraints
; ++j
) {
496 if (value_zero_p(P
->Constraint
[j
][1+i
]))
498 if (value_one_p(P
->Constraint
[j
][1+i
]))
500 if (value_mone_p(P
->Constraint
[j
][1+i
]))
504 if (j
== P
->NbConstraints
)
511 P
= Polyhedron_Copy(P
);
512 Polyhedron_ExchangeColumns(P
, 1, 1+i
);
515 if (value_zero_p(E
->d
)) {
517 subs
= ALLOCN(evalue
*, dim
);
518 for (j
= 0; j
< dim
; ++j
)
519 subs
[j
] = evalue_var(j
);
523 E
= evalue_dup(*E_p
);
524 evalue_substitute(E
, subs
);
525 for (j
= 0; j
< dim
; ++j
)
526 evalue_free(subs
[j
]);
534 static evalue
*sum_over_polytope_base(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
535 struct Bernoulli_data
*data
,
536 struct barvinok_options
*options
)
540 assert(P
->NbEq
== 0);
545 for_each_lower_upper_bound(P
, Bernoulli_init
, Bernoulli_cb
, data
);
547 res
= evalue_from_section_array(data
->s
, data
->ns
);
550 evalue
*tmp
= Bernoulli_sum_evalue(res
, nvar
-1, options
);
558 static evalue
*sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
559 struct Bernoulli_data
*data
,
560 struct barvinok_options
*options
);
562 static evalue
*sum_over_polytope_with_equalities(Polyhedron
*P
, evalue
*E
,
564 struct Bernoulli_data
*data
,
565 struct barvinok_options
*options
)
567 unsigned dim
= P
->Dimension
;
568 unsigned new_dim
, new_nparam
;
569 Matrix
*T
= NULL
, *CP
= NULL
;
576 remove_all_equalities(&P
, NULL
, &CP
, &T
, dim
-nvar
, options
->MaxRays
);
580 return evalue_zero();
583 new_nparam
= CP
? CP
->NbColumns
-1 : dim
- nvar
;
584 new_dim
= T
? T
->NbColumns
-1 : nvar
+ new_nparam
;
586 /* We can avoid these substitutions if E is a constant */
587 subs
= ALLOCN(evalue
*, dim
);
588 for (j
= 0; j
< nvar
; ++j
) {
590 subs
[j
] = affine2evalue(T
->p
[j
], T
->p
[nvar
+new_nparam
][new_dim
],
593 subs
[j
] = evalue_var(j
);
595 for (j
= 0; j
< dim
-nvar
; ++j
) {
597 subs
[nvar
+j
] = affine2evalue(CP
->p
[j
], CP
->p
[dim
-nvar
][new_nparam
],
600 subs
[nvar
+j
] = evalue_var(j
);
601 unshift(subs
[nvar
+j
], new_dim
-new_nparam
);
605 evalue_substitute(E
, subs
);
608 for (j
= 0; j
< dim
; ++j
)
609 evalue_free(subs
[j
]);
612 if (new_dim
-new_nparam
> 0) {
613 sum
= sum_over_polytope(P
, E
, new_dim
-new_nparam
, data
, options
);
619 sum
->x
.p
= new_enode(partition
, 2, new_dim
);
620 EVALUE_SET_DOMAIN(sum
->x
.p
->arr
[0], P
);
621 value_clear(sum
->x
.p
->arr
[1].d
);
622 sum
->x
.p
->arr
[1] = *E
;
627 evalue_backsubstitute(sum
, CP
, options
->MaxRays
);
637 static evalue
*sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
638 struct Bernoulli_data
*data
,
639 struct barvinok_options
*options
)
642 return sum_over_polytope_with_equalities(P
, E
, nvar
, data
, options
);
644 return sum_over_polytope_base(P
, E
, nvar
, data
, options
);
647 evalue
*Bernoulli_sum_evalue(evalue
*e
, unsigned nvar
,
648 struct barvinok_options
*options
)
650 struct Bernoulli_data data
;
652 evalue
*sum
= evalue_zero();
654 if (EVALUE_IS_ZERO(*e
))
662 assert(value_zero_p(e
->d
));
663 assert(e
->x
.p
->type
== partition
);
666 data
.s
= ALLOCN(struct evalue_section
, data
.size
);
667 data
.MaxRays
= options
->MaxRays
;
669 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
671 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
672 evalue
*E
= &e
->x
.p
->arr
[2*i
+1];
674 Polyhedron
*next
= D
->next
;
680 integer_bounds
= find_integer_bounds(&P
, &E
, nvar
);
681 if (options
->approximation_method
== BV_APPROX_NONE
&&
686 evalue
*tmp
= sum_over_polytope(P
, E
, nvar
, &data
, options
);
693 if (E
!= &e
->x
.p
->arr
[2*i
+1])
713 evalue
*Bernoulli_sum(Polyhedron
*P
, Polyhedron
*C
,
714 struct barvinok_options
*options
)
720 if (emptyQ(P
) || emptyQ(C
))
721 return evalue_zero();
723 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
724 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
729 return evalue_zero();
733 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
734 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
735 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
736 sum
= Bernoulli_sum_evalue(&e
, P
->Dimension
- C
->Dimension
, options
);
737 free_evalue_refs(&e
);