summate.c: sum_base: check equality constraints in Param_Polyhedron
[barvinok.git] / genfun.cc
blob3eac807e71f06e7047729bc933d7d49286400e0b
1 #include <iostream>
2 #include <iomanip>
3 #include <vector>
4 #include <assert.h>
5 #include <NTL/ZZ.h>
6 #include <NTL/vec_ZZ.h>
7 #include <NTL/mat_ZZ.h>
8 #include <barvinok/polylib.h>
9 #include <barvinok/genfun.h>
10 #include <barvinok/barvinok.h>
11 #include "conversion.h"
12 #include "counter.h"
13 #include "genfun_constructor.h"
14 #include "mat_util.h"
15 #include "matrix_read.h"
16 #include "remove_equalities.h"
18 using std::cout;
19 using std::cerr;
20 using std::endl;
21 using std::pair;
22 using std::vector;
24 bool short_rat_lex_smaller_denominator::operator()(const short_rat* r1,
25 const short_rat* r2) const
27 return lex_cmp(r1->d.power, r2->d.power) < 0;
30 static void lex_order_terms(struct short_rat* rat)
32 for (int i = 0; i < rat->n.power.NumRows(); ++i) {
33 int m = i;
34 for (int j = i+1; j < rat->n.power.NumRows(); ++j)
35 if (lex_cmp(rat->n.power[j], rat->n.power[m]) < 0)
36 m = j;
37 if (m != i) {
38 vec_ZZ tmp = rat->n.power[m];
39 rat->n.power[m] = rat->n.power[i];
40 rat->n.power[i] = tmp;
41 QQ tmp_coeff = rat->n.coeff[m];
42 rat->n.coeff[m] = rat->n.coeff[i];
43 rat->n.coeff[i] = tmp_coeff;
48 short_rat::short_rat(const short_rat& r)
50 n.coeff = r.n.coeff;
51 n.power = r.n.power;
52 d.power = r.d.power;
55 short_rat::short_rat(Value c)
57 n.coeff.SetLength(1);
58 value2zz(c, n.coeff[0].n);
59 n.coeff[0].d = 1;
60 n.power.SetDims(1, 0);
61 d.power.SetDims(0, 0);
64 short_rat::short_rat(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
66 n.coeff.SetLength(1);
67 ZZ g = GCD(c.n, c.d);
68 n.coeff[0].n = c.n/g;
69 n.coeff[0].d = c.d/g;
70 n.power.SetDims(1, num.length());
71 n.power[0] = num;
72 d.power = den;
73 normalize();
76 short_rat::short_rat(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den)
78 n.coeff = c;
79 n.power = num;
80 d.power = den;
81 normalize();
84 void short_rat::normalize()
86 /* Make all powers in denominator reverse-lexico-positive */
87 for (int i = 0; i < d.power.NumRows(); ++i) {
88 int j;
89 for (j = d.power.NumCols()-1; j >= 0; --j)
90 if (!IsZero(d.power[i][j]))
91 break;
92 assert(j >= 0);
93 if (sign(d.power[i][j]) < 0) {
94 negate(d.power[i], d.power[i]);
95 for (int k = 0; k < n.coeff.length(); ++k) {
96 negate(n.coeff[k].n, n.coeff[k].n);
97 n.power[k] += d.power[i];
102 /* Order powers in denominator */
103 lex_order_rows(d.power);
106 void short_rat::add(const short_rat *r)
108 for (int i = 0; i < r->n.power.NumRows(); ++i) {
109 int len = n.coeff.length();
110 int j;
111 for (j = 0; j < len; ++j)
112 if (r->n.power[i] == n.power[j])
113 break;
114 if (j < len) {
115 n.coeff[j] += r->n.coeff[i];
116 if (n.coeff[j].n == 0) {
117 if (j < len-1) {
118 n.power[j] = n.power[len-1];
119 n.coeff[j] = n.coeff[len-1];
121 int dim = n.power.NumCols();
122 n.coeff.SetLength(len-1);
123 n.power.SetDims(len-1, dim);
125 } else {
126 int dim = n.power.NumCols();
127 n.coeff.SetLength(len+1);
128 n.power.SetDims(len+1, dim);
129 n.coeff[len] = r->n.coeff[i];
130 n.power[len] = r->n.power[i];
135 QQ short_rat::coefficient(Value* params, barvinok_options *options) const
137 unsigned nvar = d.power.NumRows();
138 unsigned nparam = d.power.NumCols();
139 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
140 Value tmp;
141 value_init(tmp);
143 QQ c(0, 1);
145 for (int j = 0; j < n.coeff.length(); ++j) {
146 C->NbRows = nparam+nvar;
147 for (int r = 0; r < nparam; ++r) {
148 value_set_si(C->p[r][0], 0);
149 for (int c = 0; c < nvar; ++c) {
150 zz2value(d.power[c][r], C->p[r][1+c]);
152 zz2value(n.power[j][r], C->p[r][1+nvar]);
153 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
155 for (int r = 0; r < nvar; ++r) {
156 value_set_si(C->p[nparam+r][0], 1);
157 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
158 value_set_si(C->p[nparam+r][1+r], 1);
160 Polyhedron *P = Constraints2Polyhedron(C, options->MaxRays);
161 if (emptyQ2(P)) {
162 Polyhedron_Free(P);
163 continue;
165 barvinok_count_with_options(P, &tmp, options);
166 Polyhedron_Free(P);
167 if (value_zero_p(tmp))
168 continue;
169 QQ c2(0, 1);
170 value2zz(tmp, c2.n);
171 c2 *= n.coeff[j];
172 c += c2;
174 Matrix_Free(C);
175 value_clear(tmp);
176 return c;
179 bool short_rat::reduced()
181 int dim = n.power.NumCols();
182 lex_order_terms(this);
183 if (n.power.NumRows() % 2 == 0) {
184 if (n.coeff[0].n == -n.coeff[1].n &&
185 n.coeff[0].d == n.coeff[1].d) {
186 vec_ZZ step = n.power[1] - n.power[0];
187 int k;
188 for (k = 1; k < n.power.NumRows()/2; ++k) {
189 if (n.coeff[2*k].n != -n.coeff[2*k+1].n ||
190 n.coeff[2*k].d != n.coeff[2*k+1].d)
191 break;
192 if (step != n.power[2*k+1] - n.power[2*k])
193 break;
195 if (k == n.power.NumRows()/2) {
196 for (k = 0; k < d.power.NumRows(); ++k)
197 if (d.power[k] == step)
198 break;
199 if (k < d.power.NumRows()) {
200 for (++k; k < d.power.NumRows(); ++k)
201 d.power[k-1] = d.power[k];
202 d.power.SetDims(k-1, dim);
203 for (k = 1; k < n.power.NumRows()/2; ++k) {
204 n.coeff[k] = n.coeff[2*k];
205 n.power[k] = n.power[2*k];
207 n.coeff.SetLength(k);
208 n.power.SetDims(k, dim);
209 return true;
214 return false;
217 gen_fun::gen_fun(Value c)
219 short_rat *r = new short_rat(c);
220 context = Universe_Polyhedron(0);
221 term.insert(r);
224 void gen_fun::add(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
226 if (c.n == 0)
227 return;
229 add(new short_rat(c, num, den));
232 void gen_fun::add(short_rat *r)
234 short_rat_list::iterator i = term.find(r);
235 while (i != term.end()) {
236 (*i)->add(r);
237 if ((*i)->n.coeff.length() == 0) {
238 delete *i;
239 term.erase(i);
240 } else if ((*i)->reduced()) {
241 delete r;
242 /* we've modified term[i], so remove it
243 * and add it back again
245 r = *i;
246 term.erase(i);
247 i = term.find(r);
248 continue;
250 delete r;
251 return;
254 term.insert(r);
257 /* Extend the context of "this" to include the one of "gf".
259 void gen_fun::extend_context(const gen_fun *gf, barvinok_options *options)
261 Polyhedron *U = DomainUnion(context, gf->context, options->MaxRays);
262 Polyhedron *C = DomainConvex(U, options->MaxRays);
263 Domain_Free(U);
264 Domain_Free(context);
265 context = C;
268 /* Add the generating "gf" to "this" on the union of their domains.
270 void gen_fun::add(const QQ& c, const gen_fun *gf, barvinok_options *options)
272 extend_context(gf, options);
273 add(c, gf);
276 void gen_fun::add(const QQ& c, const gen_fun *gf)
278 QQ p;
279 for (short_rat_list::iterator i = gf->term.begin(); i != gf->term.end(); ++i) {
280 for (int j = 0; j < (*i)->n.power.NumRows(); ++j) {
281 p = c;
282 p *= (*i)->n.coeff[j];
283 add(p, (*i)->n.power[j], (*i)->d.power);
288 static void split_param_compression(Matrix *CP, mat_ZZ& map, vec_ZZ& offset)
290 Matrix *T = Transpose(CP);
291 matrix2zz(T, map, T->NbRows-1, T->NbColumns-1);
292 values2zz(T->p[T->NbRows-1], offset, T->NbColumns-1);
293 Matrix_Free(T);
297 * Perform the substitution specified by CP
299 * CP is a homogeneous matrix that maps a set of "compressed parameters"
300 * to the original set of parameters.
302 * This function is applied to a gen_fun computed with the compressed parameters
303 * and adapts it to refer to the original parameters.
305 * That is, if y are the compressed parameters and x = A y + b are the original
306 * parameters, then we want the coefficient of the monomial t^y in the original
307 * generating function to be the coefficient of the monomial u^x in the resulting
308 * generating function.
309 * The original generating function has the form
311 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
313 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
315 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
317 * = a u^{A m + b}/(1-u^{A n})
319 * Therefore, we multiply the powers m and n in both numerator and denominator by A
320 * and add b to the power in the numerator.
321 * Since the above powers are stored as row vectors m^T and n^T,
322 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
324 * The pair (map, offset) contains the same information as CP.
325 * map is the transpose of the linear part of CP, while offset is the constant part.
327 void gen_fun::substitute(Matrix *CP)
329 mat_ZZ map;
330 vec_ZZ offset;
331 split_param_compression(CP, map, offset);
332 Polyhedron *C = Polyhedron_Image(context, CP, 0);
333 Polyhedron_Free(context);
334 context = C;
336 short_rat_list new_term;
337 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
338 short_rat *r = (*i);
339 r->d.power *= map;
340 r->n.power *= map;
341 for (int j = 0; j < r->n.power.NumRows(); ++j)
342 r->n.power[j] += offset;
343 r->normalize();
344 new_term.insert(r);
346 term.swap(new_term);
349 static int Matrix_Equal(Matrix *M1, Matrix *M2)
351 int i, j;
353 if (M1->NbRows != M2->NbRows)
354 return 0;
355 if (M1->NbColumns != M2->NbColumns)
356 return 0;
357 for (i = 0; i < M1->NbRows; ++i)
358 for (j = 0; j < M1->NbColumns; ++j)
359 if (value_ne(M1->p[i][j], M2->p[i][j]))
360 return 0;
362 return 1;
365 struct parallel_cones {
366 int *pos;
367 vector<pair<Vector *, QQ> > vertices;
368 parallel_cones(int *pos) : pos(pos) {}
371 /* This structure helps in computing the generating functions
372 * of polytopes with pairwise parallel hyperplanes more efficiently.
373 * These occur when computing hadamard products of pairs of generating
374 * functions with the same denominators.
375 * If there are many such pairs then the same vertex cone
376 * may appear more than once. We therefore keep a list of all
377 * vertex cones and only compute the corresponding generating function
378 * once.
379 * However, even HPs of generating functions with the same denominators
380 * can result in polytopes of different "shapes", making them incomparable.
381 * In particular, they can have different equalities among the parameters
382 * and the variables. In such cases, only polytopes of the first "shape"
383 * that is encountered are kept in this way. The others are handled
384 * in the usual, non-optimized way.
386 struct parallel_polytopes {
387 gf_base *red;
388 Polyhedron *context;
389 Matrix *Constraints;
390 Matrix *CP, *T;
391 int dim;
392 int nparam;
393 unsigned reduced_nparam;
394 vector<parallel_cones> cones;
395 barvinok_options *options;
397 parallel_polytopes(int n, Polyhedron *context, int nparam,
398 barvinok_options *options) :
399 context(context), dim(-1), nparam(nparam),
400 options(options) {
401 red = NULL;
402 Constraints = NULL;
403 CP = NULL;
404 T = NULL;
406 bool add(const QQ& c, Polyhedron *P) {
407 int i;
409 for (i = 0; i < P->NbEq; ++i)
410 if (First_Non_Zero(P->Constraint[i]+1,
411 P->Dimension-nparam) == -1)
412 break;
413 if (i < P->NbEq)
414 return false;
416 Polyhedron *Q = remove_equalities_p(Polyhedron_Copy(P), P->Dimension-nparam,
417 NULL, options->MaxRays);
418 POL_ENSURE_VERTICES(Q);
419 if (emptyQ(Q)) {
420 Polyhedron_Free(Q);
421 return true;
424 if (Q->Dimension == 0) {
425 Polyhedron_Free(Q);
426 return false;
429 if (Q->NbEq != 0) {
430 Matrix *Q_CP;
431 Polyhedron *R = Q;
433 remove_all_equalities(&Q, NULL, &Q_CP, NULL, nparam,
434 options->MaxRays);
436 POL_ENSURE_VERTICES(Q);
437 if (emptyQ(Q) || Q->NbEq > 0 || Q->Dimension == 0) {
438 if (Q_CP)
439 Matrix_Free(Q_CP);
440 Polyhedron_Free(R);
441 Polyhedron_Free(Q);
442 return emptyQ(Q);
445 if (red) {
446 if ((!CP ^ !Q_CP) || (CP && !Matrix_Equal(CP, Q_CP))) {
447 Matrix_Free(Q_CP);
448 Polyhedron_Free(R);
449 Polyhedron_Free(Q);
450 return false;
452 Matrix_Free(Q_CP);
453 } else {
454 CP = Q_CP;
455 T = align_matrix(CP, R->Dimension+1);
458 reduced_nparam = CP->NbColumns-1;
459 Polyhedron_Free(R);
460 } else {
461 if (red && CP) {
462 Polyhedron_Free(Q);
463 return false;
465 reduced_nparam = nparam;
468 if (First_Non_Zero(Q->Constraint[Q->NbConstraints-1]+1, Q->Dimension) == -1)
469 Q->NbConstraints--;
471 if (!Constraints) {
472 Polyhedron *reduced_context;
473 dim = Q->Dimension;
474 if (CP)
475 reduced_context = Polyhedron_Preimage(context, CP, options->MaxRays);
476 else
477 reduced_context = Polyhedron_Copy(context);
478 red = gf_base::create(reduced_context, dim, reduced_nparam, options);
479 red->base->init(Q, 0);
480 Constraints = Matrix_Alloc(Q->NbConstraints, Q->Dimension);
481 for (int i = 0; i < Q->NbConstraints; ++i) {
482 Vector_Copy(Q->Constraint[i]+1, Constraints->p[i], Q->Dimension);
484 } else {
485 if (Q->Dimension != dim) {
486 Polyhedron_Free(Q);
487 return false;
489 assert(Q->Dimension == dim);
490 for (int i = 0; i < Q->NbConstraints; ++i) {
491 int j;
492 for (j = 0; j < Constraints->NbRows; ++j)
493 if (Vector_Equal(Q->Constraint[i]+1, Constraints->p[j],
494 Q->Dimension))
495 break;
496 if (j >= Constraints->NbRows) {
497 Matrix_Extend(Constraints, Constraints->NbRows+1);
498 Vector_Copy(Q->Constraint[i]+1,
499 Constraints->p[Constraints->NbRows-1],
500 Q->Dimension);
505 for (int i = 0; i < Q->NbRays; ++i) {
506 if (!value_pos_p(Q->Ray[i][dim+1]))
507 continue;
509 Polyhedron *C = supporting_cone(Q, i);
511 if (First_Non_Zero(C->Constraint[C->NbConstraints-1]+1,
512 C->Dimension) == -1)
513 C->NbConstraints--;
515 int *pos = new int[1+C->NbConstraints];
516 int l = 0;
517 for (int k = 0; k < Constraints->NbRows; ++k) {
518 for (int j = 0; j < C->NbConstraints; ++j) {
519 if (Vector_Equal(C->Constraint[j]+1, Constraints->p[k],
520 C->Dimension)) {
521 pos[1+l++] = k;
522 break;
526 pos[0] = l;
528 int j;
529 for (j = 0; j < cones.size(); ++j)
530 if (!memcmp(pos, cones[j].pos, (1+C->NbConstraints)*sizeof(int)))
531 break;
532 if (j == cones.size())
533 cones.push_back(parallel_cones(pos));
534 else
535 delete [] pos;
537 Polyhedron_Free(C);
539 int k;
540 for (k = 0; k < cones[j].vertices.size(); ++k)
541 if (Vector_Equal(Q->Ray[i]+1, cones[j].vertices[k].first->p,
542 Q->Dimension+1))
543 break;
545 if (k == cones[j].vertices.size()) {
546 Vector *vertex = Vector_Alloc(Q->Dimension+1);
547 Vector_Copy(Q->Ray[i]+1, vertex->p, Q->Dimension+1);
548 cones[j].vertices.push_back(pair<Vector*,QQ>(vertex, c));
549 } else {
550 cones[j].vertices[k].second += c;
551 if (cones[j].vertices[k].second.n == 0) {
552 int size = cones[j].vertices.size();
553 Vector_Free(cones[j].vertices[k].first);
554 if (k < size-1)
555 cones[j].vertices[k] = cones[j].vertices[size-1];
556 cones[j].vertices.pop_back();
561 Polyhedron_Free(Q);
562 return true;
564 gen_fun *compute() {
565 if (!red)
566 return NULL;
567 for (int i = 0; i < cones.size(); ++i) {
568 Matrix *M = Matrix_Alloc(cones[i].pos[0], 1+Constraints->NbColumns+1);
569 Polyhedron *Cone;
570 for (int j = 0; j <cones[i].pos[0]; ++j) {
571 value_set_si(M->p[j][0], 1);
572 Vector_Copy(Constraints->p[cones[i].pos[1+j]], M->p[j]+1,
573 Constraints->NbColumns);
575 Cone = Constraints2Polyhedron(M, options->MaxRays);
576 Matrix_Free(M);
577 for (int j = 0; j < cones[i].vertices.size(); ++j) {
578 red->base->do_vertex_cone(cones[i].vertices[j].second,
579 Polyhedron_Copy(Cone),
580 cones[i].vertices[j].first->p, options);
582 Polyhedron_Free(Cone);
584 if (CP)
585 red->gf->substitute(CP);
586 return red->gf;
588 void print(std::ostream& os) const {
589 for (int i = 0; i < cones.size(); ++i) {
590 os << "[";
591 for (int j = 0; j < cones[i].pos[0]; ++j) {
592 if (j)
593 os << ", ";
594 os << cones[i].pos[1+j];
596 os << "]" << endl;
597 for (int j = 0; j < cones[i].vertices.size(); ++j) {
598 Vector_Print(stderr, P_VALUE_FMT, cones[i].vertices[j].first);
599 os << cones[i].vertices[j].second << endl;
603 ~parallel_polytopes() {
604 for (int i = 0; i < cones.size(); ++i) {
605 delete [] cones[i].pos;
606 for (int j = 0; j < cones[i].vertices.size(); ++j)
607 Vector_Free(cones[i].vertices[j].first);
609 if (Constraints)
610 Matrix_Free(Constraints);
611 if (CP)
612 Matrix_Free(CP);
613 if (T)
614 Matrix_Free(T);
615 delete red;
619 gen_fun *gen_fun::Hadamard_product(const gen_fun *gf, barvinok_options *options)
621 QQ one(1, 1);
622 Polyhedron *C = DomainIntersection(context, gf->context, options->MaxRays);
623 gen_fun *sum = new gen_fun(C);
625 int j = 0;
626 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i, j++) {
627 int k = 0;
628 for (short_rat_list::iterator i2 = gf->term.begin();
629 i2 != gf->term.end();
630 ++i2, k++) {
631 int d = (*i)->d.power.NumCols();
632 int k1 = (*i)->d.power.NumRows();
633 int k2 = (*i2)->d.power.NumRows();
634 assert((*i)->d.power.NumCols() == (*i2)->d.power.NumCols());
636 if (options->verbose)
637 fprintf(stderr, "HP: %d/%zd %d/%zd \r",
638 j, term.size(), k, gf->term.size());
640 parallel_polytopes pp((*i)->n.power.NumRows() *
641 (*i2)->n.power.NumRows(),
642 sum->context, d, options);
644 for (int j = 0; j < (*i)->n.power.NumRows(); ++j) {
645 for (int j2 = 0; j2 < (*i2)->n.power.NumRows(); ++j2) {
646 Matrix *M = Matrix_Alloc(k1+k2+d+d, 1+k1+k2+d+1);
647 for (int k = 0; k < k1+k2; ++k) {
648 value_set_si(M->p[k][0], 1);
649 value_set_si(M->p[k][1+k], 1);
651 for (int k = 0; k < d; ++k) {
652 value_set_si(M->p[k1+k2+k][1+k1+k2+k], -1);
653 zz2value((*i)->n.power[j][k], M->p[k1+k2+k][1+k1+k2+d]);
654 for (int l = 0; l < k1; ++l)
655 zz2value((*i)->d.power[l][k], M->p[k1+k2+k][1+l]);
657 for (int k = 0; k < d; ++k) {
658 value_set_si(M->p[k1+k2+d+k][1+k1+k2+k], -1);
659 zz2value((*i2)->n.power[j2][k],
660 M->p[k1+k2+d+k][1+k1+k2+d]);
661 for (int l = 0; l < k2; ++l)
662 zz2value((*i2)->d.power[l][k],
663 M->p[k1+k2+d+k][1+k1+l]);
665 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
666 Matrix_Free(M);
668 QQ c = (*i)->n.coeff[j];
669 c *= (*i2)->n.coeff[j2];
670 if (!pp.add(c, P)) {
671 gen_fun *t = barvinok_enumerate_series(P, C->Dimension, options);
672 sum->add(c, t);
673 delete t;
676 Polyhedron_Free(P);
680 gen_fun *t = pp.compute();
681 if (t) {
682 sum->add(one, t);
683 delete t;
687 return sum;
690 /* Assuming "this" and "gf" are generating functions for sets,
691 * replace "this" by the generating function for the union of these sets.
693 * In particular, compute
695 * this + gf - gen_fun(intersection of sets)
697 void gen_fun::add_union(gen_fun *gf, barvinok_options *options)
699 QQ one(1, 1), mone(-1, 1);
701 gen_fun *hp = Hadamard_product(gf, options);
702 extend_context(gf, options);
703 add(one, gf);
704 add(mone, hp);
705 delete hp;
708 static void Polyhedron_Shift(Polyhedron *P, Vector *offset)
710 Value tmp;
711 value_init(tmp);
712 for (int i = 0; i < P->NbConstraints; ++i) {
713 Inner_Product(P->Constraint[i]+1, offset->p, P->Dimension, &tmp);
714 value_subtract(P->Constraint[i][1+P->Dimension],
715 P->Constraint[i][1+P->Dimension], tmp);
717 for (int i = 0; i < P->NbRays; ++i) {
718 if (value_notone_p(P->Ray[i][0]))
719 continue;
720 if (value_zero_p(P->Ray[i][1+P->Dimension]))
721 continue;
722 Vector_Combine(P->Ray[i]+1, offset->p, P->Ray[i]+1,
723 P->Ray[i][0], P->Ray[i][1+P->Dimension], P->Dimension);
725 value_clear(tmp);
728 void gen_fun::shift(const vec_ZZ& offset)
730 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
731 for (int j = 0; j < (*i)->n.power.NumRows(); ++j)
732 (*i)->n.power[j] += offset;
734 Vector *v = Vector_Alloc(offset.length());
735 zz2values(offset, v->p);
736 Polyhedron_Shift(context, v);
737 Vector_Free(v);
740 /* Divide the generating functin by 1/(1-z^power).
741 * The effect on the corresponding explicit function f(x) is
742 * f'(x) = \sum_{i=0}^\infty f(x - i * power)
744 void gen_fun::divide(const vec_ZZ& power)
746 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
747 int r = (*i)->d.power.NumRows();
748 int c = (*i)->d.power.NumCols();
749 (*i)->d.power.SetDims(r+1, c);
750 (*i)->d.power[r] = power;
753 Vector *v = Vector_Alloc(1+power.length()+1);
754 value_set_si(v->p[0], 1);
755 zz2values(power, v->p+1);
756 Polyhedron *C = AddRays(v->p, 1, context, context->NbConstraints+1);
757 Vector_Free(v);
758 Polyhedron_Free(context);
759 context = C;
762 static void print_power(std::ostream& os, const QQ& c, const vec_ZZ& p,
763 unsigned int nparam, const char **param_name)
765 bool first = true;
767 for (int i = 0; i < p.length(); ++i) {
768 if (p[i] == 0)
769 continue;
770 if (first) {
771 if (c.n == -1 && c.d == 1)
772 os << "-";
773 else if (c.n != 1 || c.d != 1) {
774 os << c.n;
775 if (c.d != 1)
776 os << "/" << c.d;
777 os << "*";
779 first = false;
780 } else
781 os << "*";
782 if (i < nparam)
783 os << param_name[i];
784 else
785 os << "x" << i;
786 if (p[i] == 1)
787 continue;
788 if (p[i] < 0)
789 os << "^(" << p[i] << ")";
790 else
791 os << "^" << p[i];
793 if (first) {
794 os << c.n;
795 if (c.d != 1)
796 os << "/" << c.d;
800 void short_rat::print(std::ostream& os, unsigned int nparam,
801 const char **param_name) const
803 QQ mone(-1, 1);
804 os << "(";
805 for (int j = 0; j < n.coeff.length(); ++j) {
806 if (j != 0 && n.coeff[j].n >= 0)
807 os << "+";
808 print_power(os, n.coeff[j], n.power[j], nparam, param_name);
810 os << ")";
811 if (d.power.NumRows() == 0)
812 return;
813 os << "/(";
814 for (int j = 0; j < d.power.NumRows(); ++j) {
815 if (j != 0)
816 os << " * ";
817 os << "(1";
818 print_power(os, mone, d.power[j], nparam, param_name);
819 os << ")";
821 os << ")";
824 void gen_fun::print(std::ostream& os, unsigned int nparam,
825 const char **param_name) const
827 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
828 if (i != term.begin())
829 os << " + ";
830 (*i)->print(os, nparam, param_name);
834 std::ostream & operator<< (std::ostream & os, const short_rat& r)
836 os << r.n.coeff << endl;
837 os << r.n.power << endl;
838 os << r.d.power << endl;
839 return os;
842 extern "C" { typedef void (*gmp_free_t)(void *, size_t); }
844 std::ostream & operator<< (std::ostream & os, const Polyhedron& P)
846 char *str;
847 gmp_free_t gmp_free;
848 mp_get_memory_functions(NULL, NULL, &gmp_free);
849 os << P.NbConstraints << " " << P.Dimension+2 << endl;
850 for (int i = 0; i < P.NbConstraints; ++i) {
851 for (int j = 0; j < P.Dimension+2; ++j) {
852 str = mpz_get_str(0, 10, P.Constraint[i][j]);
853 os << std::setw(4) << str << " ";
854 (*gmp_free)(str, strlen(str)+1);
856 os << endl;
858 return os;
861 std::ostream & operator<< (std::ostream & os, const gen_fun& gf)
863 os << *gf.context << endl;
864 os << endl;
865 os << gf.term.size() << endl;
866 for (short_rat_list::iterator i = gf.term.begin(); i != gf.term.end(); ++i)
867 os << **i;
868 return os;
871 gen_fun *gen_fun::read(std::istream& is, barvinok_options *options)
873 Matrix *M = Matrix_Read(is);
874 Polyhedron *C = Constraints2Polyhedron(M, options->MaxRays);
875 Matrix_Free(M);
877 gen_fun *gf = new gen_fun(C);
879 int n;
880 is >> n;
882 vec_QQ c;
883 mat_ZZ num;
884 mat_ZZ den;
885 for (int i = 0; i < n; ++i) {
886 is >> c >> num >> den;
887 gf->add(new short_rat(c, num, den));
890 return gf;
893 gen_fun::operator evalue *() const
895 evalue *EP = NULL;
896 evalue factor;
897 value_init(factor.d);
898 value_init(factor.x.n);
899 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
900 unsigned nvar = (*i)->d.power.NumRows();
901 unsigned nparam = (*i)->d.power.NumCols();
902 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
903 mat_ZZ& d = (*i)->d.power;
904 Polyhedron *U = context;
906 for (int j = 0; j < (*i)->n.coeff.length(); ++j) {
907 for (int r = 0; r < nparam; ++r) {
908 value_set_si(C->p[r][0], 0);
909 for (int c = 0; c < nvar; ++c) {
910 zz2value(d[c][r], C->p[r][1+c]);
912 Vector_Set(&C->p[r][1+nvar], 0, nparam);
913 value_set_si(C->p[r][1+nvar+r], -1);
914 zz2value((*i)->n.power[j][r], C->p[r][1+nvar+nparam]);
916 for (int r = 0; r < nvar; ++r) {
917 value_set_si(C->p[nparam+r][0], 1);
918 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
919 value_set_si(C->p[nparam+r][1+r], 1);
921 Polyhedron *P = Constraints2Polyhedron(C, 0);
922 evalue *E = barvinok_enumerate_ev(P, U, 0);
923 Polyhedron_Free(P);
924 if (EVALUE_IS_ZERO(*E)) {
925 evalue_free(E);
926 continue;
928 zz2value((*i)->n.coeff[j].n, factor.x.n);
929 zz2value((*i)->n.coeff[j].d, factor.d);
930 emul(&factor, E);
931 if (!EP)
932 EP = E;
933 else {
934 eadd(E, EP);
935 evalue_free(E);
938 Matrix_Free(C);
940 value_clear(factor.d);
941 value_clear(factor.x.n);
942 return EP ? EP : evalue_zero();
945 ZZ gen_fun::coefficient(Value* params, barvinok_options *options) const
947 if (!in_domain(context, params))
948 return ZZ::zero();
950 QQ sum(0, 1);
952 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
953 sum += (*i)->coefficient(params, options);
955 assert(sum.d == 1);
956 return sum.n;
959 void gen_fun::coefficient(Value* params, Value* c) const
961 barvinok_options *options = barvinok_options_new_with_defaults();
963 ZZ coeff = coefficient(params, options);
965 zz2value(coeff, *c);
967 barvinok_options_free(options);
970 gen_fun *gen_fun::summate(int nvar, barvinok_options *options) const
972 int n_try;
973 int dim = context->Dimension;
974 int nparam = dim - nvar;
975 reducer *red;
976 gen_fun *gf;
978 if (nparam == 0) {
979 bool finite;
980 Value c;
981 value_init(c);
982 finite = summate(&c);
983 assert(finite);
984 gf = new gen_fun(c);
985 value_clear(c);
986 return gf;
989 if (options->incremental_specialization == 1) {
990 red = new partial_ireducer(Polyhedron_Project(context, nparam), dim, nparam);
991 } else
992 red = new partial_reducer(Polyhedron_Project(context, nparam), dim, nparam);
993 n_try = 0;
994 for (;;) {
995 try {
996 red->init(context, n_try);
997 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
998 red->reduce((*i)->n.coeff, (*i)->n.power, (*i)->d.power);
999 break;
1000 } catch (OrthogonalException &e) {
1001 red->reset();
1002 n_try++;
1005 gf = red->get_gf();
1006 delete red;
1007 return gf;
1010 /* returns true if the set was finite and false otherwise */
1011 bool gen_fun::summate(Value *sum) const
1013 if (term.size() == 0) {
1014 value_set_si(*sum, 0);
1015 return true;
1018 int maxlen = 0;
1019 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
1020 if ((*i)->d.power.NumRows() > maxlen)
1021 maxlen = (*i)->d.power.NumRows();
1023 infinite_counter cnt((*term.begin())->d.power.NumCols(), maxlen);
1024 cnt.init(context, 0);
1025 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
1026 cnt.reduce((*i)->n.coeff, (*i)->n.power, (*i)->d.power);
1028 for (int i = 1; i <= maxlen; ++i)
1029 if (value_notzero_p(mpq_numref(cnt.count[i]))) {
1030 value_set_si(*sum, -1);
1031 return false;
1034 assert(value_one_p(mpq_denref(cnt.count[0])));
1035 value_assign(*sum, mpq_numref(cnt.count[0]));
1036 return true;
1039 bool gen_fun::is_zero() const
1041 bool empty;
1042 Value c;
1044 value_init(c);
1046 empty = summate(&c) && value_zero_p(c);
1048 value_clear(c);
1050 return empty;