volume.c: fix typo in comment
[barvinok.git] / bernoulli.c
blobc67e87832114ac2d200975ac41ef685da6968e2a
1 #include <barvinok/options.h>
2 #include <barvinok/util.h>
3 #include "bernoulli.h"
5 #define ALLOC(type) (type*)malloc(sizeof(type))
6 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
8 static struct bernoulli bernoulli;
10 struct bernoulli *bernoulli_compute(int n)
12 int i, j;
13 Value lcm, factor, tmp;
15 if (n < bernoulli.n)
16 return &bernoulli;
18 if (n >= bernoulli.size) {
19 int size = n + 5;
20 Vector *b;
21 Vector **sum;
23 b = Vector_Alloc(size);
24 if (bernoulli.n) {
25 Vector_Copy(bernoulli.b_num->p, b->p, bernoulli.n);
26 Vector_Free(bernoulli.b_num);
28 bernoulli.b_num = b;
29 b = Vector_Alloc(size);
30 if (bernoulli.n) {
31 Vector_Copy(bernoulli.b_den->p, b->p, bernoulli.n);
32 Vector_Free(bernoulli.b_den);
34 bernoulli.b_den = b;
35 sum = ALLOCN(Vector *, size);
36 for (i = 0; i < bernoulli.n; ++i)
37 sum[i] = bernoulli.sum[i];
38 free(bernoulli.sum);
39 bernoulli.sum = sum;
41 bernoulli.size = size;
43 value_init(lcm);
44 value_init(factor);
45 value_init(tmp);
46 value_set_si(lcm, 1);
47 for (i = 0; i < bernoulli.n; ++i)
48 value_lcm(lcm, bernoulli.b_den->p[i], &lcm);
49 for (i = bernoulli.n; i <= n; ++i) {
50 if (i == 0) {
51 value_set_si(bernoulli.b_num->p[0], 1);
52 value_set_si(bernoulli.b_den->p[0], 1);
53 continue;
55 value_set_si(bernoulli.b_num->p[i], 0);
56 value_set_si(factor, -(i+1));
57 for (j = i-1; j >= 0; --j) {
58 mpz_mul_ui(factor, factor, j+1);
59 mpz_divexact_ui(factor, factor, i+1-j);
60 value_division(tmp, lcm, bernoulli.b_den->p[j]);
61 value_multiply(tmp, tmp, bernoulli.b_num->p[j]);
62 value_multiply(tmp, tmp, factor);
63 value_addto(bernoulli.b_num->p[i], bernoulli.b_num->p[i], tmp);
65 mpz_mul_ui(bernoulli.b_den->p[i], lcm, i+1);
66 Gcd(bernoulli.b_num->p[i], bernoulli.b_den->p[i], &tmp);
67 if (value_notone_p(tmp)) {
68 value_division(bernoulli.b_num->p[i], bernoulli.b_num->p[i], tmp);
69 value_division(bernoulli.b_den->p[i], bernoulli.b_den->p[i], tmp);
71 value_lcm(lcm, bernoulli.b_den->p[i], &lcm);
73 for (i = bernoulli.n; i <= n; ++i) {
74 bernoulli.sum[i] = Vector_Alloc(i+3);
75 value_assign(bernoulli.sum[i]->p[i+1], lcm);
76 mpz_mul_ui(bernoulli.sum[i]->p[i+2], lcm, i+1);
77 value_set_si(factor, 1);
78 for (j = 1; j <= i; ++j) {
79 mpz_mul_ui(factor, factor, i+2-j);
80 mpz_divexact_ui(factor, factor, j);
81 value_division(bernoulli.sum[i]->p[i+1-j], lcm, bernoulli.b_den->p[j]);
82 value_multiply(bernoulli.sum[i]->p[i+1-j],
83 bernoulli.sum[i]->p[i+1-j], bernoulli.b_num->p[j]);
84 value_multiply(bernoulli.sum[i]->p[i+1-j],
85 bernoulli.sum[i]->p[i+1-j], factor);
87 Vector_Normalize(bernoulli.sum[i]->p, i+3);
89 bernoulli.n = n+1;
90 value_clear(lcm);
91 value_clear(factor);
92 value_clear(tmp);
94 return &bernoulli;
97 /* shift variables in polynomial one down */
98 static void shift(evalue *e)
100 int i;
101 if (value_notzero_p(e->d))
102 return;
103 assert(e->x.p->type == polynomial);
104 assert(e->x.p->pos > 1);
105 e->x.p->pos--;
106 for (i = 0; i < e->x.p->size; ++i)
107 shift(&e->x.p->arr[i]);
110 static evalue *shifted_copy(evalue *src)
112 evalue *e = ALLOC(evalue);
113 value_init(e->d);
114 evalue_copy(e, src);
115 shift(e);
116 return e;
119 static evalue *power_sums(struct bernoulli *bernoulli, evalue *poly,
120 Vector *arg, Value denom, int neg, int alt_neg)
122 int i;
123 evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
124 evalue *sum = evalue_zero();
126 for (i = 1; i < poly->x.p->size; ++i) {
127 evalue *term = evalue_polynomial(bernoulli->sum[i], base);
128 evalue *factor = shifted_copy(&poly->x.p->arr[i]);
129 emul(factor, term);
130 if (alt_neg && (i % 2))
131 evalue_negate(term);
132 eadd(term, sum);
133 free_evalue_refs(factor);
134 free_evalue_refs(term);
135 free(factor);
136 free(term);
138 if (neg)
139 evalue_negate(sum);
140 free_evalue_refs(base);
141 free(base);
143 return sum;
146 struct section { Polyhedron *D; evalue *E; };
148 struct Bernoulli_data {
149 unsigned MaxRays;
150 struct section *s;
151 int size;
152 int ns;
153 evalue *e;
156 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
158 struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
159 Matrix *M2;
160 Polyhedron *T;
161 evalue *factor = NULL;
162 evalue *linear = NULL;
163 int constant = 0;
164 Value tmp;
165 unsigned dim = M->NbColumns-2;
166 Vector *row;
168 assert(data->ns < data->size);
170 M2 = Matrix_Copy(M);
171 T = Constraints2Polyhedron(M2, data->MaxRays);
172 Matrix_Free(M2);
174 POL_ENSURE_VERTICES(T);
175 if (emptyQ(T)) {
176 Polyhedron_Free(T);
177 return;
180 assert(lower != upper);
182 row = Vector_Alloc(dim+1);
183 value_init(tmp);
184 if (value_notzero_p(data->e->d)) {
185 factor = data->e;
186 constant = 1;
187 } else {
188 assert(data->e->x.p->type == polynomial);
189 if (data->e->x.p->pos > 1) {
190 factor = shifted_copy(data->e);
191 constant = 1;
192 } else
193 factor = shifted_copy(&data->e->x.p->arr[0]);
195 if (!EVALUE_IS_ZERO(*factor)) {
196 value_absolute(tmp, upper[1]);
197 /* upper - lower */
198 Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
199 value_multiply(tmp, tmp, lower[1]);
200 /* upper - lower + 1 */
201 value_addto(row->p[dim], row->p[dim], tmp);
203 linear = affine2evalue(row->p, tmp, dim);
204 emul(factor, linear);
205 } else
206 linear = evalue_zero();
208 if (constant) {
209 data->s[data->ns].E = linear;
210 data->s[data->ns].D = T;
211 ++data->ns;
212 } else {
213 evalue *poly_u = NULL, *poly_l = NULL;
214 Polyhedron *D;
215 struct bernoulli *bernoulli;
216 assert(data->e->x.p->type == polynomial);
217 assert(data->e->x.p->pos == 1);
218 bernoulli = bernoulli_compute(data->e->x.p->size-1);
219 /* lower is the constraint
220 * b i - l' >= 0 i >= l'/b = l
221 * upper is the constraint
222 * -a i + u' >= 0 i <= u'/a = u
224 M2 = Matrix_Alloc(3, 2+T->Dimension);
225 value_set_si(M2->p[0][0], 1);
226 value_set_si(M2->p[1][0], 1);
227 value_set_si(M2->p[2][0], 1);
228 /* Case 1:
229 * 1 <= l l' - b >= 0
231 Vector_Oppose(lower+2, M2->p[0]+1, T->Dimension+1);
232 value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
233 lower[1]);
234 D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
235 if (emptyQ2(D))
236 Polyhedron_Free(D);
237 else {
238 evalue *extra;
239 if (!poly_u) {
240 Vector_Copy(upper+2, row->p, dim+1);
241 value_oppose(tmp, upper[1]);
242 value_addto(row->p[dim], row->p[dim], tmp);
243 poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0);
245 Vector_Oppose(lower+2, row->p, dim+1);
246 extra = power_sums(bernoulli, data->e, row, lower[1], 1, 0);
247 eadd(poly_u, extra);
248 eadd(linear, extra);
250 data->s[data->ns].E = extra;
251 data->s[data->ns].D = D;
252 ++data->ns;
255 /* Case 2:
256 * 1 <= -u -u' - a >= 0
258 Vector_Oppose(upper+2, M2->p[0]+1, T->Dimension+1);
259 value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
260 upper[1]);
261 D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
262 if (emptyQ2(D))
263 Polyhedron_Free(D);
264 else {
265 evalue *extra;
266 if (!poly_l) {
267 Vector_Copy(lower+2, row->p, dim+1);
268 value_addto(row->p[dim], row->p[dim], lower[1]);
269 poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1);
271 Vector_Oppose(upper+2, row->p, dim+1);
272 value_oppose(tmp, upper[1]);
273 extra = power_sums(bernoulli, data->e, row, tmp, 1, 1);
274 eadd(poly_l, extra);
275 eadd(linear, extra);
277 data->s[data->ns].E = extra;
278 data->s[data->ns].D = D;
279 ++data->ns;
282 /* Case 3:
283 * u >= 0 u' >= 0
284 * -l >= 0 -l' >= 0
286 Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
287 Vector_Copy(lower+2, M2->p[1]+1, T->Dimension+1);
288 D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
289 if (emptyQ2(D))
290 Polyhedron_Free(D);
291 else {
292 if (!poly_l) {
293 Vector_Copy(lower+2, row->p, dim+1);
294 value_addto(row->p[dim], row->p[dim], lower[1]);
295 poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1);
297 if (!poly_u) {
298 Vector_Copy(upper+2, row->p, dim+1);
299 value_oppose(tmp, upper[1]);
300 value_addto(row->p[dim], row->p[dim], tmp);
301 poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0);
304 data->s[data->ns].E = ALLOC(evalue);
305 value_init(data->s[data->ns].E->d);
306 evalue_copy(data->s[data->ns].E, poly_u);
307 eadd(poly_l, data->s[data->ns].E);
308 eadd(linear, data->s[data->ns].E);
309 data->s[data->ns].D = D;
310 ++data->ns;
313 /* Case 4:
314 * l < 1 -l' + b - 1 >= 0
315 * 0 < l l' - 1 >= 0
317 Vector_Copy(lower+2, M2->p[0]+1, T->Dimension+1);
318 value_addto(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension], lower[1]);
319 value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
320 Vector_Oppose(lower+2, M2->p[1]+1, T->Dimension+1);
321 value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
322 D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
323 if (emptyQ2(D))
324 Polyhedron_Free(D);
325 else {
326 if (!poly_u) {
327 Vector_Copy(upper+2, row->p, dim+1);
328 value_oppose(tmp, upper[1]);
329 value_addto(row->p[dim], row->p[dim], tmp);
330 poly_u = power_sums(bernoulli, data->e, row, tmp, 0, 0);
333 eadd(linear, poly_u);
334 data->s[data->ns].E = poly_u;
335 poly_u = NULL;
336 data->s[data->ns].D = D;
337 ++data->ns;
340 /* Case 5:
341 * -u < 1 u' + a - 1 >= 0
342 * 0 < -u -u' - 1 >= 0
343 * l <= 0 -l' >= 0
345 Vector_Copy(upper+2, M2->p[0]+1, T->Dimension+1);
346 value_subtract(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension],
347 upper[1]);
348 value_decrement(M2->p[0][1+T->Dimension], M2->p[0][1+T->Dimension]);
349 Vector_Oppose(upper+2, M2->p[1]+1, T->Dimension+1);
350 value_decrement(M2->p[1][1+T->Dimension], M2->p[1][1+T->Dimension]);
351 Vector_Copy(lower+2, M2->p[2]+1, T->Dimension+1);
352 D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
353 if (emptyQ2(D))
354 Polyhedron_Free(D);
355 else {
356 if (!poly_l) {
357 Vector_Copy(lower+2, row->p, dim+1);
358 value_addto(row->p[dim], row->p[dim], lower[1]);
359 poly_l = power_sums(bernoulli, data->e, row, lower[1], 0, 1);
362 eadd(linear, poly_l);
363 data->s[data->ns].E = poly_l;
364 poly_l = NULL;
365 data->s[data->ns].D = D;
366 ++data->ns;
369 Matrix_Free(M2);
370 Polyhedron_Free(T);
371 if (poly_l) {
372 free_evalue_refs(poly_l);
373 free(poly_l);
375 if (poly_u) {
376 free_evalue_refs(poly_u);
377 free(poly_u);
379 free_evalue_refs(linear);
380 free(linear);
382 if (factor != data->e) {
383 free_evalue_refs(factor);
384 free(factor);
386 value_clear(tmp);
387 Vector_Free(row);
390 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
391 struct barvinok_options *options)
393 struct Bernoulli_data data;
394 int i, j;
395 evalue *sum = evalue_zero();
397 if (EVALUE_IS_ZERO(*e))
398 return sum;
400 if (nvar == 0) {
401 eadd(e, sum);
402 return sum;
405 assert(value_zero_p(e->d));
406 assert(e->x.p->type == partition);
408 data.size = e->x.p->size * 2 + 16;
409 data.s = ALLOCN(struct section, data.size);
410 data.MaxRays = options->MaxRays;
412 for (i = 0; i < e->x.p->size/2; ++i) {
413 Polyhedron *D;
414 for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
415 unsigned dim = D->Dimension - 1;
416 Polyhedron *next = D->next;
417 evalue tmp;
418 D->next = NULL;
420 value_init(tmp.d);
421 value_set_si(tmp.d, 0);
423 if (value_zero_p(D->Constraint[0][0]) &&
424 value_notzero_p(D->Constraint[0][1])) {
425 tmp.x.p = new_enode(partition, 2, dim);
426 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Project(D, dim));
427 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
428 reduce_evalue_in_domain(&tmp.x.p->arr[1], D);
429 shift(&tmp.x.p->arr[1]);
430 } else {
431 data.ns = 0;
432 data.e = &e->x.p->arr[2*i+1];
434 for_each_lower_upper_bound(D, Bernoulli_cb, &data);
436 if (data.ns == 0)
437 evalue_set_si(&tmp, 0, 1);
438 else {
439 tmp.x.p = new_enode(partition, 2*data.ns, dim);
440 for (j = 0; j < data.ns; ++j) {
441 EVALUE_SET_DOMAIN(tmp.x.p->arr[2*j], data.s[j].D);
442 value_clear(tmp.x.p->arr[2*j+1].d);
443 tmp.x.p->arr[2*j+1] = *data.s[j].E;
444 free(data.s[j].E);
449 if (nvar > 1) {
450 evalue *res = Bernoulli_sum_evalue(&tmp, nvar-1, options);
451 eadd(res, sum);
452 free_evalue_refs(res);
453 free(res);
454 } else
455 eadd(&tmp, sum);
457 free_evalue_refs(&tmp);
458 D->next = next;;
462 free(data.s);
464 reduce_evalue(sum);
465 return sum;
468 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
469 struct barvinok_options *options)
471 evalue e;
472 evalue *sum;
474 value_init(e.d);
475 e.x.p = new_enode(partition, 2, P->Dimension);
476 EVALUE_SET_DOMAIN(e.x.p->arr[0], Polyhedron_Copy(P));
477 evalue_set_si(&e.x.p->arr[1], 1, 1);
478 sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
479 free_evalue_refs(&e);
480 return sum;