series: leave freeing of P argument to calling function
[barvinok.git] / genfun.cc
blob7f4263b3db05c81209af9a59287c670220f20e55
1 #include <iostream>
2 #include <iomanip>
3 #include <vector>
4 #include <assert.h>
5 #include <barvinok/genfun.h>
6 #include <barvinok/barvinok.h>
7 #include "conversion.h"
8 #include "counter.h"
9 #include "genfun_constructor.h"
10 #include "mat_util.h"
11 #include "matrix_read.h"
12 #include "remove_equalities.h"
14 using std::cout;
15 using std::cerr;
16 using std::endl;
17 using std::pair;
18 using std::vector;
20 bool short_rat_lex_smaller_denominator::operator()(const short_rat* r1,
21 const short_rat* r2) const
23 return lex_cmp(r1->d.power, r2->d.power) < 0;
26 static void lex_order_terms(struct short_rat* rat)
28 for (int i = 0; i < rat->n.power.NumRows(); ++i) {
29 int m = i;
30 for (int j = i+1; j < rat->n.power.NumRows(); ++j)
31 if (lex_cmp(rat->n.power[j], rat->n.power[m]) < 0)
32 m = j;
33 if (m != i) {
34 vec_ZZ tmp = rat->n.power[m];
35 rat->n.power[m] = rat->n.power[i];
36 rat->n.power[i] = tmp;
37 QQ tmp_coeff = rat->n.coeff[m];
38 rat->n.coeff[m] = rat->n.coeff[i];
39 rat->n.coeff[i] = tmp_coeff;
44 short_rat::short_rat(const short_rat& r)
46 n.coeff = r.n.coeff;
47 n.power = r.n.power;
48 d.power = r.d.power;
51 short_rat::short_rat(Value c)
53 n.coeff.SetLength(1);
54 value2zz(c, n.coeff[0].n);
55 n.coeff[0].d = 1;
56 n.power.SetDims(1, 0);
57 d.power.SetDims(0, 0);
60 short_rat::short_rat(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
62 n.coeff.SetLength(1);
63 ZZ g = GCD(c.n, c.d);
64 n.coeff[0].n = c.n/g;
65 n.coeff[0].d = c.d/g;
66 n.power.SetDims(1, num.length());
67 n.power[0] = num;
68 d.power = den;
69 normalize();
72 short_rat::short_rat(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den)
74 n.coeff = c;
75 n.power = num;
76 d.power = den;
77 normalize();
80 void short_rat::normalize()
82 /* Make all powers in denominator reverse-lexico-positive */
83 for (int i = 0; i < d.power.NumRows(); ++i) {
84 int j;
85 for (j = d.power.NumCols()-1; j >= 0; --j)
86 if (!IsZero(d.power[i][j]))
87 break;
88 assert(j >= 0);
89 if (sign(d.power[i][j]) < 0) {
90 negate(d.power[i], d.power[i]);
91 for (int k = 0; k < n.coeff.length(); ++k) {
92 negate(n.coeff[k].n, n.coeff[k].n);
93 n.power[k] += d.power[i];
98 /* Order powers in denominator */
99 lex_order_rows(d.power);
102 void short_rat::add(const short_rat *r)
104 for (int i = 0; i < r->n.power.NumRows(); ++i) {
105 int len = n.coeff.length();
106 int j;
107 for (j = 0; j < len; ++j)
108 if (r->n.power[i] == n.power[j])
109 break;
110 if (j < len) {
111 n.coeff[j] += r->n.coeff[i];
112 if (n.coeff[j].n == 0) {
113 if (j < len-1) {
114 n.power[j] = n.power[len-1];
115 n.coeff[j] = n.coeff[len-1];
117 int dim = n.power.NumCols();
118 n.coeff.SetLength(len-1);
119 n.power.SetDims(len-1, dim);
121 } else {
122 int dim = n.power.NumCols();
123 n.coeff.SetLength(len+1);
124 n.power.SetDims(len+1, dim);
125 n.coeff[len] = r->n.coeff[i];
126 n.power[len] = r->n.power[i];
131 QQ short_rat::coefficient(Value* params, barvinok_options *options) const
133 unsigned nvar = d.power.NumRows();
134 unsigned nparam = d.power.NumCols();
135 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
136 Value tmp;
137 value_init(tmp);
139 QQ c(0, 1);
141 for (int j = 0; j < n.coeff.length(); ++j) {
142 C->NbRows = nparam+nvar;
143 for (int r = 0; r < nparam; ++r) {
144 value_set_si(C->p[r][0], 0);
145 for (int c = 0; c < nvar; ++c) {
146 zz2value(d.power[c][r], C->p[r][1+c]);
148 zz2value(n.power[j][r], C->p[r][1+nvar]);
149 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
151 for (int r = 0; r < nvar; ++r) {
152 value_set_si(C->p[nparam+r][0], 1);
153 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
154 value_set_si(C->p[nparam+r][1+r], 1);
156 Polyhedron *P = Constraints2Polyhedron(C, options->MaxRays);
157 if (emptyQ2(P)) {
158 Polyhedron_Free(P);
159 continue;
161 barvinok_count_with_options(P, &tmp, options);
162 Polyhedron_Free(P);
163 if (value_zero_p(tmp))
164 continue;
165 QQ c2(0, 1);
166 value2zz(tmp, c2.n);
167 c2 *= n.coeff[j];
168 c += c2;
170 Matrix_Free(C);
171 value_clear(tmp);
172 return c;
175 bool short_rat::reduced()
177 int dim = n.power.NumCols();
178 lex_order_terms(this);
179 if (n.power.NumRows() % 2 == 0) {
180 if (n.coeff[0].n == -n.coeff[1].n &&
181 n.coeff[0].d == n.coeff[1].d) {
182 vec_ZZ step = n.power[1] - n.power[0];
183 int k;
184 for (k = 1; k < n.power.NumRows()/2; ++k) {
185 if (n.coeff[2*k].n != -n.coeff[2*k+1].n ||
186 n.coeff[2*k].d != n.coeff[2*k+1].d)
187 break;
188 if (step != n.power[2*k+1] - n.power[2*k])
189 break;
191 if (k == n.power.NumRows()/2) {
192 for (k = 0; k < d.power.NumRows(); ++k)
193 if (d.power[k] == step)
194 break;
195 if (k < d.power.NumRows()) {
196 for (++k; k < d.power.NumRows(); ++k)
197 d.power[k-1] = d.power[k];
198 d.power.SetDims(k-1, dim);
199 for (k = 1; k < n.power.NumRows()/2; ++k) {
200 n.coeff[k] = n.coeff[2*k];
201 n.power[k] = n.power[2*k];
203 n.coeff.SetLength(k);
204 n.power.SetDims(k, dim);
205 return true;
210 return false;
213 gen_fun::gen_fun(Value c)
215 short_rat *r = new short_rat(c);
216 context = Universe_Polyhedron(0);
217 term.insert(r);
220 void gen_fun::add(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
222 if (c.n == 0)
223 return;
225 add(new short_rat(c, num, den));
228 void gen_fun::add(short_rat *r)
230 short_rat_list::iterator i = term.find(r);
231 while (i != term.end()) {
232 (*i)->add(r);
233 if ((*i)->n.coeff.length() == 0) {
234 delete *i;
235 term.erase(i);
236 } else if ((*i)->reduced()) {
237 delete r;
238 /* we've modified term[i], so remove it
239 * and add it back again
241 r = *i;
242 term.erase(i);
243 i = term.find(r);
244 continue;
246 delete r;
247 return;
250 term.insert(r);
253 void gen_fun::add(const QQ& c, const gen_fun *gf, barvinok_options *options)
255 Polyhedron *U = DomainUnion(context, gf->context, options->MaxRays);
256 Polyhedron *C = DomainConvex(U, options->MaxRays);
257 Domain_Free(U);
258 Domain_Free(context);
259 context = C;
261 add(c, gf);
264 void gen_fun::add(const QQ& c, const gen_fun *gf)
266 QQ p;
267 for (short_rat_list::iterator i = gf->term.begin(); i != gf->term.end(); ++i) {
268 for (int j = 0; j < (*i)->n.power.NumRows(); ++j) {
269 p = c;
270 p *= (*i)->n.coeff[j];
271 add(p, (*i)->n.power[j], (*i)->d.power);
276 static void split_param_compression(Matrix *CP, mat_ZZ& map, vec_ZZ& offset)
278 Matrix *T = Transpose(CP);
279 matrix2zz(T, map, T->NbRows-1, T->NbColumns-1);
280 values2zz(T->p[T->NbRows-1], offset, T->NbColumns-1);
281 Matrix_Free(T);
285 * Perform the substitution specified by CP
287 * CP is a homogeneous matrix that maps a set of "compressed parameters"
288 * to the original set of parameters.
290 * This function is applied to a gen_fun computed with the compressed parameters
291 * and adapts it to refer to the original parameters.
293 * That is, if y are the compressed parameters and x = A y + b are the original
294 * parameters, then we want the coefficient of the monomial t^y in the original
295 * generating function to be the coefficient of the monomial u^x in the resulting
296 * generating function.
297 * The original generating function has the form
299 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
301 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
303 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
305 * = a u^{A m + b}/(1-u^{A n})
307 * Therefore, we multiply the powers m and n in both numerator and denominator by A
308 * and add b to the power in the numerator.
309 * Since the above powers are stored as row vectors m^T and n^T,
310 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
312 * The pair (map, offset) contains the same information as CP.
313 * map is the transpose of the linear part of CP, while offset is the constant part.
315 void gen_fun::substitute(Matrix *CP)
317 mat_ZZ map;
318 vec_ZZ offset;
319 split_param_compression(CP, map, offset);
320 Polyhedron *C = Polyhedron_Image(context, CP, 0);
321 Polyhedron_Free(context);
322 context = C;
324 short_rat_list new_term;
325 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
326 short_rat *r = (*i);
327 r->d.power *= map;
328 r->n.power *= map;
329 for (int j = 0; j < r->n.power.NumRows(); ++j)
330 r->n.power[j] += offset;
331 r->normalize();
332 new_term.insert(r);
334 term.swap(new_term);
337 static int Matrix_Equal(Matrix *M1, Matrix *M2)
339 int i, j;
341 if (M1->NbRows != M2->NbRows)
342 return 0;
343 if (M1->NbColumns != M2->NbColumns)
344 return 0;
345 for (i = 0; i < M1->NbRows; ++i)
346 for (j = 0; j < M1->NbColumns; ++j)
347 if (value_ne(M1->p[i][j], M2->p[i][j]))
348 return 0;
350 return 1;
353 struct parallel_cones {
354 int *pos;
355 vector<pair<Vector *, QQ> > vertices;
356 parallel_cones(int *pos) : pos(pos) {}
359 /* This structure helps in computing the generating functions
360 * of polytopes with pairwise parallel hyperplanes more efficiently.
361 * These occur when computing hadamard products of pairs of generating
362 * functions with the same denominators.
363 * If there are many such pairs then the same vertex cone
364 * may appear more than once. We therefore keep a list of all
365 * vertex cones and only compute the corresponding generating function
366 * once.
367 * However, even HPs of generating functions with the same denominators
368 * can result in polytopes of different "shapes", making them incomparable.
369 * In particular, they can have different equalities among the parameters
370 * and the variables. In such cases, only polytopes of the first "shape"
371 * that is encountered are kept in this way. The others are handled
372 * in the usual, non-optimized way.
374 struct parallel_polytopes {
375 gf_base *red;
376 Polyhedron *context;
377 Matrix *Constraints;
378 Matrix *CP, *T;
379 int dim;
380 int nparam;
381 unsigned reduced_nparam;
382 vector<parallel_cones> cones;
383 barvinok_options *options;
385 parallel_polytopes(int n, Polyhedron *context, int nparam,
386 barvinok_options *options) :
387 context(context), dim(-1), nparam(nparam),
388 options(options) {
389 red = NULL;
390 Constraints = NULL;
391 CP = NULL;
392 T = NULL;
394 bool add(const QQ& c, Polyhedron *P) {
395 int i;
397 for (i = 0; i < P->NbEq; ++i)
398 if (First_Non_Zero(P->Constraint[i]+1,
399 P->Dimension-nparam) == -1)
400 break;
401 if (i < P->NbEq)
402 return false;
404 Polyhedron *Q = remove_equalities_p(Polyhedron_Copy(P), P->Dimension-nparam,
405 NULL, options->MaxRays);
406 POL_ENSURE_VERTICES(Q);
407 if (emptyQ(Q)) {
408 Polyhedron_Free(Q);
409 return true;
412 if (Q->NbEq != 0) {
413 Matrix *Q_CP;
414 Polyhedron *R = Q;
416 remove_all_equalities(&Q, NULL, &Q_CP, NULL, nparam,
417 options->MaxRays);
419 POL_ENSURE_VERTICES(Q);
420 if (emptyQ(Q) || Q->NbEq > 0) {
421 if (Q_CP)
422 Matrix_Free(Q_CP);
423 Polyhedron_Free(R);
424 Polyhedron_Free(Q);
425 return emptyQ(Q);
428 if (red) {
429 if ((!CP ^ !Q_CP) || (CP && !Matrix_Equal(CP, Q_CP))) {
430 Matrix_Free(Q_CP);
431 Polyhedron_Free(R);
432 Polyhedron_Free(Q);
433 return false;
435 Matrix_Free(Q_CP);
436 } else {
437 CP = Q_CP;
438 T = align_matrix(CP, R->Dimension+1);
441 reduced_nparam = CP->NbColumns-1;
442 Polyhedron_Free(R);
443 } else {
444 if (red && CP) {
445 Polyhedron_Free(Q);
446 return false;
448 reduced_nparam = nparam;
451 if (First_Non_Zero(Q->Constraint[Q->NbConstraints-1]+1, Q->Dimension) == -1)
452 Q->NbConstraints--;
454 if (!Constraints) {
455 Polyhedron *reduced_context;
456 dim = Q->Dimension;
457 if (CP)
458 reduced_context = Polyhedron_Preimage(context, CP, options->MaxRays);
459 else
460 reduced_context = Polyhedron_Copy(context);
461 red = gf_base::create(reduced_context, dim, reduced_nparam, options);
462 red->base->init(Q);
463 Constraints = Matrix_Alloc(Q->NbConstraints, Q->Dimension);
464 for (int i = 0; i < Q->NbConstraints; ++i) {
465 Vector_Copy(Q->Constraint[i]+1, Constraints->p[i], Q->Dimension);
467 } else {
468 if (Q->Dimension != dim) {
469 Polyhedron_Free(Q);
470 return false;
472 assert(Q->Dimension == dim);
473 for (int i = 0; i < Q->NbConstraints; ++i) {
474 int j;
475 for (j = 0; j < Constraints->NbRows; ++j)
476 if (Vector_Equal(Q->Constraint[i]+1, Constraints->p[j],
477 Q->Dimension))
478 break;
479 if (j >= Constraints->NbRows) {
480 Matrix_Extend(Constraints, Constraints->NbRows+1);
481 Vector_Copy(Q->Constraint[i]+1,
482 Constraints->p[Constraints->NbRows-1],
483 Q->Dimension);
488 for (int i = 0; i < Q->NbRays; ++i) {
489 if (!value_pos_p(Q->Ray[i][dim+1]))
490 continue;
492 Polyhedron *C = supporting_cone(Q, i);
494 if (First_Non_Zero(C->Constraint[C->NbConstraints-1]+1,
495 C->Dimension) == -1)
496 C->NbConstraints--;
498 int *pos = new int[1+C->NbConstraints];
499 int l = 0;
500 for (int k = 0; k < Constraints->NbRows; ++k) {
501 for (int j = 0; j < C->NbConstraints; ++j) {
502 if (Vector_Equal(C->Constraint[j]+1, Constraints->p[k],
503 C->Dimension)) {
504 pos[1+l++] = k;
505 break;
509 pos[0] = l;
511 int j;
512 for (j = 0; j < cones.size(); ++j)
513 if (!memcmp(pos, cones[j].pos, (1+C->NbConstraints)*sizeof(int)))
514 break;
515 if (j == cones.size())
516 cones.push_back(parallel_cones(pos));
517 else
518 delete [] pos;
520 Polyhedron_Free(C);
522 int k;
523 for (k = 0; k < cones[j].vertices.size(); ++k)
524 if (Vector_Equal(Q->Ray[i]+1, cones[j].vertices[k].first->p,
525 Q->Dimension+1))
526 break;
528 if (k == cones[j].vertices.size()) {
529 Vector *vertex = Vector_Alloc(Q->Dimension+1);
530 Vector_Copy(Q->Ray[i]+1, vertex->p, Q->Dimension+1);
531 cones[j].vertices.push_back(pair<Vector*,QQ>(vertex, c));
532 } else {
533 cones[j].vertices[k].second += c;
534 if (cones[j].vertices[k].second.n == 0) {
535 int size = cones[j].vertices.size();
536 Vector_Free(cones[j].vertices[k].first);
537 if (k < size-1)
538 cones[j].vertices[k] = cones[j].vertices[size-1];
539 cones[j].vertices.pop_back();
544 Polyhedron_Free(Q);
545 return true;
547 gen_fun *compute() {
548 if (!red)
549 return NULL;
550 for (int i = 0; i < cones.size(); ++i) {
551 Matrix *M = Matrix_Alloc(cones[i].pos[0], 1+Constraints->NbColumns+1);
552 Polyhedron *Cone;
553 for (int j = 0; j <cones[i].pos[0]; ++j) {
554 value_set_si(M->p[j][0], 1);
555 Vector_Copy(Constraints->p[cones[i].pos[1+j]], M->p[j]+1,
556 Constraints->NbColumns);
558 Cone = Constraints2Polyhedron(M, options->MaxRays);
559 Matrix_Free(M);
560 for (int j = 0; j < cones[i].vertices.size(); ++j) {
561 red->base->do_vertex_cone(cones[i].vertices[j].second,
562 Polyhedron_Copy(Cone),
563 cones[i].vertices[j].first->p, options);
565 Polyhedron_Free(Cone);
567 if (CP)
568 red->gf->substitute(CP);
569 return red->gf;
571 void print(std::ostream& os) const {
572 for (int i = 0; i < cones.size(); ++i) {
573 os << "[";
574 for (int j = 0; j < cones[i].pos[0]; ++j) {
575 if (j)
576 os << ", ";
577 os << cones[i].pos[1+j];
579 os << "]" << endl;
580 for (int j = 0; j < cones[i].vertices.size(); ++j) {
581 Vector_Print(stderr, P_VALUE_FMT, cones[i].vertices[j].first);
582 os << cones[i].vertices[j].second << endl;
586 ~parallel_polytopes() {
587 for (int i = 0; i < cones.size(); ++i) {
588 delete [] cones[i].pos;
589 for (int j = 0; j < cones[i].vertices.size(); ++j)
590 Vector_Free(cones[i].vertices[j].first);
592 if (Constraints)
593 Matrix_Free(Constraints);
594 if (CP)
595 Matrix_Free(CP);
596 if (T)
597 Matrix_Free(T);
598 delete red;
602 gen_fun *gen_fun::Hadamard_product(const gen_fun *gf, barvinok_options *options)
604 QQ one(1, 1);
605 Polyhedron *C = DomainIntersection(context, gf->context, options->MaxRays);
606 Polyhedron *U = Universe_Polyhedron(C->Dimension);
607 gen_fun *sum = new gen_fun(C);
609 int j = 0;
610 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i, j++) {
611 int k = 0;
612 for (short_rat_list::iterator i2 = gf->term.begin();
613 i2 != gf->term.end();
614 ++i2, k++) {
615 int d = (*i)->d.power.NumCols();
616 int k1 = (*i)->d.power.NumRows();
617 int k2 = (*i2)->d.power.NumRows();
618 assert((*i)->d.power.NumCols() == (*i2)->d.power.NumCols());
620 if (options->verbose)
621 fprintf(stderr, "HP: %d/%d %d/%d \r",
622 j, term.size(), k, gf->term.size());
624 parallel_polytopes pp((*i)->n.power.NumRows() *
625 (*i2)->n.power.NumRows(),
626 sum->context, d, options);
628 for (int j = 0; j < (*i)->n.power.NumRows(); ++j) {
629 for (int j2 = 0; j2 < (*i2)->n.power.NumRows(); ++j2) {
630 Matrix *M = Matrix_Alloc(k1+k2+d+d, 1+k1+k2+d+1);
631 for (int k = 0; k < k1+k2; ++k) {
632 value_set_si(M->p[k][0], 1);
633 value_set_si(M->p[k][1+k], 1);
635 for (int k = 0; k < d; ++k) {
636 value_set_si(M->p[k1+k2+k][1+k1+k2+k], -1);
637 zz2value((*i)->n.power[j][k], M->p[k1+k2+k][1+k1+k2+d]);
638 for (int l = 0; l < k1; ++l)
639 zz2value((*i)->d.power[l][k], M->p[k1+k2+k][1+l]);
641 for (int k = 0; k < d; ++k) {
642 value_set_si(M->p[k1+k2+d+k][1+k1+k2+k], -1);
643 zz2value((*i2)->n.power[j2][k],
644 M->p[k1+k2+d+k][1+k1+k2+d]);
645 for (int l = 0; l < k2; ++l)
646 zz2value((*i2)->d.power[l][k],
647 M->p[k1+k2+d+k][1+k1+l]);
649 Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
650 Matrix_Free(M);
652 QQ c = (*i)->n.coeff[j];
653 c *= (*i2)->n.coeff[j2];
654 if (!pp.add(c, P)) {
655 gen_fun *t = barvinok_series_with_options(P, U, options);
656 sum->add(c, t);
657 delete t;
660 Polyhedron_Free(P);
664 gen_fun *t = pp.compute();
665 if (t) {
666 sum->add(one, t);
667 delete t;
671 Polyhedron_Free(U);
672 return sum;
675 void gen_fun::add_union(gen_fun *gf, barvinok_options *options)
677 QQ one(1, 1), mone(-1, 1);
679 gen_fun *hp = Hadamard_product(gf, options);
680 add(one, gf);
681 add(mone, hp);
682 delete hp;
685 static void Polyhedron_Shift(Polyhedron *P, Vector *offset)
687 Value tmp;
688 value_init(tmp);
689 for (int i = 0; i < P->NbConstraints; ++i) {
690 Inner_Product(P->Constraint[i]+1, offset->p, P->Dimension, &tmp);
691 value_subtract(P->Constraint[i][1+P->Dimension],
692 P->Constraint[i][1+P->Dimension], tmp);
694 for (int i = 0; i < P->NbRays; ++i) {
695 if (value_notone_p(P->Ray[i][0]))
696 continue;
697 if (value_zero_p(P->Ray[i][1+P->Dimension]))
698 continue;
699 Vector_Combine(P->Ray[i]+1, offset->p, P->Ray[i]+1,
700 P->Ray[i][0], P->Ray[i][1+P->Dimension], P->Dimension);
702 value_clear(tmp);
705 void gen_fun::shift(const vec_ZZ& offset)
707 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
708 for (int j = 0; j < (*i)->n.power.NumRows(); ++j)
709 (*i)->n.power[j] += offset;
711 Vector *v = Vector_Alloc(offset.length());
712 zz2values(offset, v->p);
713 Polyhedron_Shift(context, v);
714 Vector_Free(v);
717 /* Divide the generating functin by 1/(1-z^power).
718 * The effect on the corresponding explicit function f(x) is
719 * f'(x) = \sum_{i=0}^\infty f(x - i * power)
721 void gen_fun::divide(const vec_ZZ& power)
723 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
724 int r = (*i)->d.power.NumRows();
725 int c = (*i)->d.power.NumCols();
726 (*i)->d.power.SetDims(r+1, c);
727 (*i)->d.power[r] = power;
730 Vector *v = Vector_Alloc(1+power.length()+1);
731 value_set_si(v->p[0], 1);
732 zz2values(power, v->p+1);
733 Polyhedron *C = AddRays(v->p, 1, context, context->NbConstraints+1);
734 Vector_Free(v);
735 Polyhedron_Free(context);
736 context = C;
739 static void print_power(std::ostream& os, const QQ& c, const vec_ZZ& p,
740 unsigned int nparam, char **param_name)
742 bool first = true;
744 for (int i = 0; i < p.length(); ++i) {
745 if (p[i] == 0)
746 continue;
747 if (first) {
748 if (c.n == -1 && c.d == 1)
749 os << "-";
750 else if (c.n != 1 || c.d != 1) {
751 os << c.n;
752 if (c.d != 1)
753 os << "/" << c.d;
754 os << "*";
756 first = false;
757 } else
758 os << "*";
759 if (i < nparam)
760 os << param_name[i];
761 else
762 os << "x" << i;
763 if (p[i] == 1)
764 continue;
765 if (p[i] < 0)
766 os << "^(" << p[i] << ")";
767 else
768 os << "^" << p[i];
770 if (first) {
771 os << c.n;
772 if (c.d != 1)
773 os << "/" << c.d;
777 void short_rat::print(std::ostream& os, unsigned int nparam, char **param_name) const
779 QQ mone(-1, 1);
780 os << "(";
781 for (int j = 0; j < n.coeff.length(); ++j) {
782 if (j != 0 && n.coeff[j].n >= 0)
783 os << "+";
784 print_power(os, n.coeff[j], n.power[j], nparam, param_name);
786 os << ")/(";
787 for (int j = 0; j < d.power.NumRows(); ++j) {
788 if (j != 0)
789 os << " * ";
790 os << "(1";
791 print_power(os, mone, d.power[j], nparam, param_name);
792 os << ")";
794 os << ")";
797 void gen_fun::print(std::ostream& os, unsigned int nparam, char **param_name) const
799 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
800 if (i != term.begin())
801 os << " + ";
802 (*i)->print(os, nparam, param_name);
806 std::ostream & operator<< (std::ostream & os, const short_rat& r)
808 os << r.n.coeff << endl;
809 os << r.n.power << endl;
810 os << r.d.power << endl;
811 return os;
814 std::ostream & operator<< (std::ostream & os, const Polyhedron& P)
816 char *str;
817 void (*gmp_free)(void *, size_t);
818 mp_get_memory_functions(NULL, NULL, &gmp_free);
819 os << P.NbConstraints << " " << P.Dimension+2 << endl;
820 for (int i = 0; i < P.NbConstraints; ++i) {
821 for (int j = 0; j < P.Dimension+2; ++j) {
822 str = mpz_get_str(0, 10, P.Constraint[i][j]);
823 os << std::setw(4) << str << " ";
824 (*gmp_free)(str, strlen(str)+1);
826 os << endl;
828 return os;
831 std::ostream & operator<< (std::ostream & os, const gen_fun& gf)
833 os << *gf.context << endl;
834 os << endl;
835 os << gf.term.size() << endl;
836 for (short_rat_list::iterator i = gf.term.begin(); i != gf.term.end(); ++i)
837 os << **i;
838 return os;
841 gen_fun *gen_fun::read(std::istream& is, barvinok_options *options)
843 Matrix *M = Matrix_Read(is);
844 Polyhedron *C = Constraints2Polyhedron(M, options->MaxRays);
845 Matrix_Free(M);
847 gen_fun *gf = new gen_fun(C);
849 int n;
850 is >> n;
852 vec_QQ c;
853 mat_ZZ num;
854 mat_ZZ den;
855 for (int i = 0; i < n; ++i) {
856 is >> c >> num >> den;
857 gf->add(new short_rat(c, num, den));
860 return gf;
863 gen_fun::operator evalue *() const
865 evalue *EP = NULL;
866 evalue factor;
867 value_init(factor.d);
868 value_init(factor.x.n);
869 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
870 unsigned nvar = (*i)->d.power.NumRows();
871 unsigned nparam = (*i)->d.power.NumCols();
872 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
873 mat_ZZ& d = (*i)->d.power;
874 Polyhedron *U = context;
876 for (int j = 0; j < (*i)->n.coeff.length(); ++j) {
877 for (int r = 0; r < nparam; ++r) {
878 value_set_si(C->p[r][0], 0);
879 for (int c = 0; c < nvar; ++c) {
880 zz2value(d[c][r], C->p[r][1+c]);
882 Vector_Set(&C->p[r][1+nvar], 0, nparam);
883 value_set_si(C->p[r][1+nvar+r], -1);
884 zz2value((*i)->n.power[j][r], C->p[r][1+nvar+nparam]);
886 for (int r = 0; r < nvar; ++r) {
887 value_set_si(C->p[nparam+r][0], 1);
888 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
889 value_set_si(C->p[nparam+r][1+r], 1);
891 Polyhedron *P = Constraints2Polyhedron(C, 0);
892 evalue *E = barvinok_enumerate_ev(P, U, 0);
893 Polyhedron_Free(P);
894 if (EVALUE_IS_ZERO(*E)) {
895 evalue_free(E);
896 continue;
898 zz2value((*i)->n.coeff[j].n, factor.x.n);
899 zz2value((*i)->n.coeff[j].d, factor.d);
900 emul(&factor, E);
901 if (!EP)
902 EP = E;
903 else {
904 eadd(E, EP);
905 evalue_free(E);
908 Matrix_Free(C);
910 value_clear(factor.d);
911 value_clear(factor.x.n);
912 return EP ? EP : evalue_zero();
915 ZZ gen_fun::coefficient(Value* params, barvinok_options *options) const
917 if (!in_domain(context, params))
918 return ZZ::zero();
920 QQ sum(0, 1);
922 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
923 sum += (*i)->coefficient(params, options);
925 assert(sum.d == 1);
926 return sum.n;
929 void gen_fun::coefficient(Value* params, Value* c) const
931 barvinok_options *options = barvinok_options_new_with_defaults();
933 ZZ coeff = coefficient(params, options);
935 zz2value(coeff, *c);
937 barvinok_options_free(options);
940 gen_fun *gen_fun::summate(int nvar, barvinok_options *options) const
942 int dim = context->Dimension;
943 int nparam = dim - nvar;
944 reducer *red;
945 gen_fun *gf;
947 if (nparam == 0) {
948 bool finite;
949 Value c;
950 value_init(c);
951 finite = summate(&c);
952 assert(finite);
953 gf = new gen_fun(c);
954 value_clear(c);
955 return gf;
958 if (options->incremental_specialization == 1) {
959 red = new partial_ireducer(Polyhedron_Project(context, nparam), dim, nparam);
960 } else
961 red = new partial_reducer(Polyhedron_Project(context, nparam), dim, nparam);
962 for (;;) {
963 try {
964 red->init(context);
965 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
966 red->reduce((*i)->n.coeff, (*i)->n.power, (*i)->d.power);
967 break;
968 } catch (OrthogonalException &e) {
969 red->reset();
972 gf = red->get_gf();
973 delete red;
974 return gf;
977 /* returns true if the set was finite and false otherwise */
978 bool gen_fun::summate(Value *sum) const
980 if (term.size() == 0) {
981 value_set_si(*sum, 0);
982 return true;
985 int maxlen = 0;
986 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
987 if ((*i)->d.power.NumRows() > maxlen)
988 maxlen = (*i)->d.power.NumRows();
990 infinite_counter cnt((*term.begin())->d.power.NumCols(), maxlen);
991 cnt.init(context);
992 for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
993 cnt.reduce((*i)->n.coeff, (*i)->n.power, (*i)->d.power);
995 for (int i = 1; i <= maxlen; ++i)
996 if (value_notzero_p(mpq_numref(cnt.count[i]))) {
997 value_set_si(*sum, -1);
998 return false;
1001 assert(value_one_p(mpq_denref(cnt.count[0])));
1002 value_assign(*sum, mpq_numref(cnt.count[0]));
1003 return true;