5 #include <barvinok/genfun.h>
6 #include <barvinok/barvinok.h>
7 #include "conversion.h"
8 #include "genfun_constructor.h"
17 bool short_rat_lex_smaller_denominator::operator()(const short_rat
* r1
,
18 const short_rat
* r2
) const
20 return lex_cmp(r1
->d
.power
, r2
->d
.power
) < 0;
23 static void lex_order_terms(struct short_rat
* rat
)
25 for (int i
= 0; i
< rat
->n
.power
.NumRows(); ++i
) {
27 for (int j
= i
+1; j
< rat
->n
.power
.NumRows(); ++j
)
28 if (lex_cmp(rat
->n
.power
[j
], rat
->n
.power
[m
]) < 0)
31 vec_ZZ tmp
= rat
->n
.power
[m
];
32 rat
->n
.power
[m
] = rat
->n
.power
[i
];
33 rat
->n
.power
[i
] = tmp
;
34 QQ tmp_coeff
= rat
->n
.coeff
[m
];
35 rat
->n
.coeff
[m
] = rat
->n
.coeff
[i
];
36 rat
->n
.coeff
[i
] = tmp_coeff
;
41 short_rat::short_rat(const short_rat
& r
)
48 short_rat::short_rat(Value c
)
51 value2zz(c
, n
.coeff
[0].n
);
53 n
.power
.SetDims(1, 0);
54 d
.power
.SetDims(0, 0);
57 short_rat::short_rat(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den
)
63 n
.power
.SetDims(1, num
.length());
69 short_rat::short_rat(const vec_QQ
& c
, const mat_ZZ
& num
, const mat_ZZ
& den
)
77 void short_rat::normalize()
79 /* Make all powers in denominator reverse-lexico-positive */
80 for (int i
= 0; i
< d
.power
.NumRows(); ++i
) {
82 for (j
= d
.power
.NumCols()-1; j
>= 0; --j
)
83 if (!IsZero(d
.power
[i
][j
]))
86 if (sign(d
.power
[i
][j
]) < 0) {
87 negate(d
.power
[i
], d
.power
[i
]);
88 for (int k
= 0; k
< n
.coeff
.length(); ++k
) {
89 negate(n
.coeff
[k
].n
, n
.coeff
[k
].n
);
90 n
.power
[k
] += d
.power
[i
];
95 /* Order powers in denominator */
96 lex_order_rows(d
.power
);
99 void short_rat::add(const short_rat
*r
)
101 for (int i
= 0; i
< r
->n
.power
.NumRows(); ++i
) {
102 int len
= n
.coeff
.length();
104 for (j
= 0; j
< len
; ++j
)
105 if (r
->n
.power
[i
] == n
.power
[j
])
108 n
.coeff
[j
] += r
->n
.coeff
[i
];
109 if (n
.coeff
[j
].n
== 0) {
111 n
.power
[j
] = n
.power
[len
-1];
112 n
.coeff
[j
] = n
.coeff
[len
-1];
114 int dim
= n
.power
.NumCols();
115 n
.coeff
.SetLength(len
-1);
116 n
.power
.SetDims(len
-1, dim
);
119 int dim
= n
.power
.NumCols();
120 n
.coeff
.SetLength(len
+1);
121 n
.power
.SetDims(len
+1, dim
);
122 n
.coeff
[len
] = r
->n
.coeff
[i
];
123 n
.power
[len
] = r
->n
.power
[i
];
128 QQ
short_rat::coefficient(Value
* params
, barvinok_options
*options
) const
130 unsigned nvar
= d
.power
.NumRows();
131 unsigned nparam
= d
.power
.NumCols();
132 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ 1);
138 for (int j
= 0; j
< n
.coeff
.length(); ++j
) {
139 C
->NbRows
= nparam
+nvar
;
140 for (int r
= 0; r
< nparam
; ++r
) {
141 value_set_si(C
->p
[r
][0], 0);
142 for (int c
= 0; c
< nvar
; ++c
) {
143 zz2value(d
.power
[c
][r
], C
->p
[r
][1+c
]);
145 zz2value(n
.power
[j
][r
], C
->p
[r
][1+nvar
]);
146 value_subtract(C
->p
[r
][1+nvar
], C
->p
[r
][1+nvar
], params
[r
]);
148 for (int r
= 0; r
< nvar
; ++r
) {
149 value_set_si(C
->p
[nparam
+r
][0], 1);
150 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ 1);
151 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
153 Polyhedron
*P
= Constraints2Polyhedron(C
, options
->MaxRays
);
158 barvinok_count_with_options(P
, &tmp
, options
);
160 if (value_zero_p(tmp
))
172 bool short_rat::reduced()
174 int dim
= n
.power
.NumCols();
175 lex_order_terms(this);
176 if (n
.power
.NumRows() % 2 == 0) {
177 if (n
.coeff
[0].n
== -n
.coeff
[1].n
&&
178 n
.coeff
[0].d
== n
.coeff
[1].d
) {
179 vec_ZZ step
= n
.power
[1] - n
.power
[0];
181 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
182 if (n
.coeff
[2*k
].n
!= -n
.coeff
[2*k
+1].n
||
183 n
.coeff
[2*k
].d
!= n
.coeff
[2*k
+1].d
)
185 if (step
!= n
.power
[2*k
+1] - n
.power
[2*k
])
188 if (k
== n
.power
.NumRows()/2) {
189 for (k
= 0; k
< d
.power
.NumRows(); ++k
)
190 if (d
.power
[k
] == step
)
192 if (k
< d
.power
.NumRows()) {
193 for (++k
; k
< d
.power
.NumRows(); ++k
)
194 d
.power
[k
-1] = d
.power
[k
];
195 d
.power
.SetDims(k
-1, dim
);
196 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
197 n
.coeff
[k
] = n
.coeff
[2*k
];
198 n
.power
[k
] = n
.power
[2*k
];
200 n
.coeff
.SetLength(k
);
201 n
.power
.SetDims(k
, dim
);
210 gen_fun::gen_fun(Value c
)
212 short_rat
*r
= new short_rat(c
);
213 context
= Universe_Polyhedron(0);
217 void gen_fun::add(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den
)
222 add(new short_rat(c
, num
, den
));
225 void gen_fun::add(short_rat
*r
)
227 short_rat_list::iterator i
= term
.find(r
);
228 while (i
!= term
.end()) {
230 if ((*i
)->n
.coeff
.length() == 0) {
233 } else if ((*i
)->reduced()) {
235 /* we've modified term[i], so remove it
236 * and add it back again
250 void gen_fun::add(const QQ
& c
, const gen_fun
*gf
)
253 for (short_rat_list::iterator i
= gf
->term
.begin(); i
!= gf
->term
.end(); ++i
) {
254 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
) {
256 p
*= (*i
)->n
.coeff
[j
];
257 add(p
, (*i
)->n
.power
[j
], (*i
)->d
.power
);
262 static void split_param_compression(Matrix
*CP
, mat_ZZ
& map
, vec_ZZ
& offset
)
264 Matrix
*T
= Transpose(CP
);
265 matrix2zz(T
, map
, T
->NbRows
-1, T
->NbColumns
-1);
266 values2zz(T
->p
[T
->NbRows
-1], offset
, T
->NbColumns
-1);
271 * Perform the substitution specified by CP
273 * CP is a homogeneous matrix that maps a set of "compressed parameters"
274 * to the original set of parameters.
276 * This function is applied to a gen_fun computed with the compressed parameters
277 * and adapts it to refer to the original parameters.
279 * That is, if y are the compressed parameters and x = A y + b are the original
280 * parameters, then we want the coefficient of the monomial t^y in the original
281 * generating function to be the coefficient of the monomial u^x in the resulting
282 * generating function.
283 * The original generating function has the form
285 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
287 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
289 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
291 * = a u^{A m + b}/(1-u^{A n})
293 * Therefore, we multiply the powers m and n in both numerator and denominator by A
294 * and add b to the power in the numerator.
295 * Since the above powers are stored as row vectors m^T and n^T,
296 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
298 * The pair (map, offset) contains the same information as CP.
299 * map is the transpose of the linear part of CP, while offset is the constant part.
301 void gen_fun::substitute(Matrix
*CP
)
305 split_param_compression(CP
, map
, offset
);
306 Polyhedron
*C
= Polyhedron_Image(context
, CP
, 0);
307 Polyhedron_Free(context
);
310 short_rat_list new_term
;
311 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
315 for (int j
= 0; j
< r
->n
.power
.NumRows(); ++j
)
316 r
->n
.power
[j
] += offset
;
323 struct parallel_cones
{
325 vector
<pair
<Vector
*, QQ
> > vertices
;
326 parallel_cones(int *pos
) : pos(pos
) {}
329 struct parallel_polytopes
{
336 vector
<parallel_cones
> cones
;
337 barvinok_options
*options
;
339 parallel_polytopes(int n
, Polyhedron
*context
, int nparam
,
340 barvinok_options
*options
) :
341 context(context
), dim(-1), nparam(nparam
),
348 bool add(const QQ
& c
, Polyhedron
*P
) {
351 for (i
= 0; i
< P
->NbEq
; ++i
)
352 if (First_Non_Zero(P
->Constraint
[i
]+1,
353 P
->Dimension
-nparam
) == -1)
358 Polyhedron
*Q
= remove_equalities_p(Polyhedron_Copy(P
), P
->Dimension
-nparam
,
359 NULL
, options
->MaxRays
);
360 POL_ENSURE_VERTICES(Q
);
370 M
= Matrix_Alloc(Q
->NbEq
, Q
->Dimension
+2);
371 Vector_Copy(Q
->Constraint
[0], M
->p
[0], Q
->NbEq
* (Q
->Dimension
+2));
372 CP
= compress_parms(M
, nparam
);
373 T
= align_matrix(CP
, Q
->Dimension
+1);
376 R
= Polyhedron_Preimage(Q
, T
, options
->MaxRays
);
378 Q
= remove_equalities_p(R
, R
->Dimension
-nparam
, NULL
,
381 assert(Q
->NbEq
== 0);
383 if (First_Non_Zero(Q
->Constraint
[Q
->NbConstraints
-1]+1, Q
->Dimension
) == -1)
388 red
= gf_base::create(Polyhedron_Copy(context
), dim
, nparam
, options
);
390 Constraints
= Matrix_Alloc(Q
->NbConstraints
, Q
->Dimension
);
391 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
392 Vector_Copy(Q
->Constraint
[i
]+1, Constraints
->p
[i
], Q
->Dimension
);
395 assert(Q
->Dimension
== dim
);
396 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
398 for (j
= 0; j
< Constraints
->NbRows
; ++j
)
399 if (Vector_Equal(Q
->Constraint
[i
]+1, Constraints
->p
[j
],
402 assert(j
< Constraints
->NbRows
);
406 for (int i
= 0; i
< Q
->NbRays
; ++i
) {
407 if (!value_pos_p(Q
->Ray
[i
][dim
+1]))
410 Polyhedron
*C
= supporting_cone(Q
, i
);
412 if (First_Non_Zero(C
->Constraint
[C
->NbConstraints
-1]+1,
416 int *pos
= new int[1+C
->NbConstraints
];
417 pos
[0] = C
->NbConstraints
;
419 for (int k
= 0; k
< Constraints
->NbRows
; ++k
) {
420 for (int j
= 0; j
< C
->NbConstraints
; ++j
) {
421 if (Vector_Equal(C
->Constraint
[j
]+1, Constraints
->p
[k
],
428 assert(l
== C
->NbConstraints
);
431 for (j
= 0; j
< cones
.size(); ++j
)
432 if (!memcmp(pos
, cones
[j
].pos
, (1+C
->NbConstraints
)*sizeof(int)))
434 if (j
== cones
.size())
435 cones
.push_back(parallel_cones(pos
));
442 for (k
= 0; k
< cones
[j
].vertices
.size(); ++k
)
443 if (Vector_Equal(Q
->Ray
[i
]+1, cones
[j
].vertices
[k
].first
->p
,
447 if (k
== cones
[j
].vertices
.size()) {
448 Vector
*vertex
= Vector_Alloc(Q
->Dimension
+1);
449 Vector_Copy(Q
->Ray
[i
]+1, vertex
->p
, Q
->Dimension
+1);
450 cones
[j
].vertices
.push_back(pair
<Vector
*,QQ
>(vertex
, c
));
452 cones
[j
].vertices
[k
].second
+= c
;
453 if (cones
[j
].vertices
[k
].second
.n
== 0) {
454 int size
= cones
[j
].vertices
.size();
455 Vector_Free(cones
[j
].vertices
[k
].first
);
457 cones
[j
].vertices
[k
] = cones
[j
].vertices
[size
-1];
458 cones
[j
].vertices
.pop_back();
469 for (int i
= 0; i
< cones
.size(); ++i
) {
470 Matrix
*M
= Matrix_Alloc(cones
[i
].pos
[0], 1+Constraints
->NbColumns
+1);
472 for (int j
= 0; j
<cones
[i
].pos
[0]; ++j
) {
473 value_set_si(M
->p
[j
][0], 1);
474 Vector_Copy(Constraints
->p
[cones
[i
].pos
[1+j
]], M
->p
[j
]+1,
475 Constraints
->NbColumns
);
477 Cone
= Constraints2Polyhedron(M
, options
->MaxRays
);
479 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
480 red
->base
->do_vertex_cone(cones
[i
].vertices
[j
].second
,
481 Polyhedron_Copy(Cone
),
482 cones
[i
].vertices
[j
].first
->p
, options
);
484 Polyhedron_Free(Cone
);
487 red
->gf
->substitute(CP
);
490 void print(std::ostream
& os
) const {
491 for (int i
= 0; i
< cones
.size(); ++i
) {
493 for (int j
= 0; j
< cones
[i
].pos
[0]; ++j
) {
496 os
<< cones
[i
].pos
[1+j
];
499 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
500 Vector_Print(stderr
, P_VALUE_FMT
, cones
[i
].vertices
[j
].first
);
501 os
<< cones
[i
].vertices
[j
].second
<< endl
;
505 ~parallel_polytopes() {
506 for (int i
= 0; i
< cones
.size(); ++i
) {
507 delete [] cones
[i
].pos
;
508 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
)
509 Vector_Free(cones
[i
].vertices
[j
].first
);
512 Matrix_Free(Constraints
);
521 gen_fun
*gen_fun::Hadamard_product(const gen_fun
*gf
, barvinok_options
*options
)
524 Polyhedron
*C
= DomainIntersection(context
, gf
->context
, options
->MaxRays
);
525 Polyhedron
*U
= Universe_Polyhedron(C
->Dimension
);
526 gen_fun
*sum
= new gen_fun(C
);
527 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
528 for (short_rat_list::iterator i2
= gf
->term
.begin(); i2
!= gf
->term
.end();
530 int d
= (*i
)->d
.power
.NumCols();
531 int k1
= (*i
)->d
.power
.NumRows();
532 int k2
= (*i2
)->d
.power
.NumRows();
533 assert((*i
)->d
.power
.NumCols() == (*i2
)->d
.power
.NumCols());
535 parallel_polytopes
pp((*i
)->n
.power
.NumRows() *
536 (*i2
)->n
.power
.NumRows(),
537 sum
->context
, d
, options
);
539 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
) {
540 for (int j2
= 0; j2
< (*i2
)->n
.power
.NumRows(); ++j2
) {
541 Matrix
*M
= Matrix_Alloc(k1
+k2
+d
+d
, 1+k1
+k2
+d
+1);
542 for (int k
= 0; k
< k1
+k2
; ++k
) {
543 value_set_si(M
->p
[k
][0], 1);
544 value_set_si(M
->p
[k
][1+k
], 1);
546 for (int k
= 0; k
< d
; ++k
) {
547 value_set_si(M
->p
[k1
+k2
+k
][1+k1
+k2
+k
], -1);
548 zz2value((*i
)->n
.power
[j
][k
], M
->p
[k1
+k2
+k
][1+k1
+k2
+d
]);
549 for (int l
= 0; l
< k1
; ++l
)
550 zz2value((*i
)->d
.power
[l
][k
], M
->p
[k1
+k2
+k
][1+l
]);
552 for (int k
= 0; k
< d
; ++k
) {
553 value_set_si(M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+k
], -1);
554 zz2value((*i2
)->n
.power
[j2
][k
],
555 M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+d
]);
556 for (int l
= 0; l
< k2
; ++l
)
557 zz2value((*i2
)->d
.power
[l
][k
],
558 M
->p
[k1
+k2
+d
+k
][1+k1
+l
]);
560 Polyhedron
*P
= Constraints2Polyhedron(M
, options
->MaxRays
);
563 QQ c
= (*i
)->n
.coeff
[j
];
564 c
*= (*i2
)->n
.coeff
[j2
];
566 gen_fun
*t
= barvinok_series_with_options(P
, U
, options
);
575 gen_fun
*t
= pp
.compute();
586 void gen_fun::add_union(gen_fun
*gf
, barvinok_options
*options
)
588 QQ
one(1, 1), mone(-1, 1);
590 gen_fun
*hp
= Hadamard_product(gf
, options
);
596 static void Polyhedron_Shift(Polyhedron
*P
, Vector
*offset
)
600 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
601 Inner_Product(P
->Constraint
[i
]+1, offset
->p
, P
->Dimension
, &tmp
);
602 value_subtract(P
->Constraint
[i
][1+P
->Dimension
],
603 P
->Constraint
[i
][1+P
->Dimension
], tmp
);
605 for (int i
= 0; i
< P
->NbRays
; ++i
) {
606 if (value_notone_p(P
->Ray
[i
][0]))
608 if (value_zero_p(P
->Ray
[i
][1+P
->Dimension
]))
610 Vector_Combine(P
->Ray
[i
]+1, offset
->p
, P
->Ray
[i
]+1,
611 P
->Ray
[i
][0], P
->Ray
[i
][1+P
->Dimension
], P
->Dimension
);
616 void gen_fun::shift(const vec_ZZ
& offset
)
618 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
619 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
)
620 (*i
)->n
.power
[j
] += offset
;
622 Vector
*v
= Vector_Alloc(offset
.length());
623 zz2values(offset
, v
->p
);
624 Polyhedron_Shift(context
, v
);
628 /* Divide the generating functin by 1/(1-z^power).
629 * The effect on the corresponding explicit function f(x) is
630 * f'(x) = \sum_{i=0}^\infty f(x - i * power)
632 void gen_fun::divide(const vec_ZZ
& power
)
634 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
635 int r
= (*i
)->d
.power
.NumRows();
636 int c
= (*i
)->d
.power
.NumCols();
637 (*i
)->d
.power
.SetDims(r
+1, c
);
638 (*i
)->d
.power
[r
] = power
;
641 Vector
*v
= Vector_Alloc(1+power
.length()+1);
642 value_set_si(v
->p
[0], 1);
643 zz2values(power
, v
->p
+1);
644 Polyhedron
*C
= AddRays(v
->p
, 1, context
, context
->NbConstraints
+1);
646 Polyhedron_Free(context
);
650 static void print_power(std::ostream
& os
, const QQ
& c
, const vec_ZZ
& p
,
651 unsigned int nparam
, char **param_name
)
655 for (int i
= 0; i
< p
.length(); ++i
) {
659 if (c
.n
== -1 && c
.d
== 1)
661 else if (c
.n
!= 1 || c
.d
!= 1) {
677 os
<< "^(" << p
[i
] << ")";
688 void short_rat::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
692 for (int j
= 0; j
< n
.coeff
.length(); ++j
) {
693 if (j
!= 0 && n
.coeff
[j
].n
> 0)
695 print_power(os
, n
.coeff
[j
], n
.power
[j
], nparam
, param_name
);
698 for (int j
= 0; j
< d
.power
.NumRows(); ++j
) {
702 print_power(os
, mone
, d
.power
[j
], nparam
, param_name
);
708 void gen_fun::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
710 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
711 if (i
!= term
.begin())
713 (*i
)->print(os
, nparam
, param_name
);
717 std::ostream
& operator<< (std::ostream
& os
, const short_rat
& r
)
719 os
<< r
.n
.coeff
<< endl
;
720 os
<< r
.n
.power
<< endl
;
721 os
<< r
.d
.power
<< endl
;
725 std::ostream
& operator<< (std::ostream
& os
, const Polyhedron
& P
)
728 void (*gmp_free
)(void *, size_t);
729 mp_get_memory_functions(NULL
, NULL
, &gmp_free
);
730 os
<< P
.NbConstraints
<< " " << P
.Dimension
+2 << endl
;
731 for (int i
= 0; i
< P
.NbConstraints
; ++i
) {
732 for (int j
= 0; j
< P
.Dimension
+2; ++j
) {
733 str
= mpz_get_str(0, 10, P
.Constraint
[i
][j
]);
734 os
<< std::setw(4) << str
<< " ";
735 (*gmp_free
)(str
, strlen(str
)+1);
742 std::ostream
& operator<< (std::ostream
& os
, const gen_fun
& gf
)
744 os
<< *gf
.context
<< endl
;
746 os
<< gf
.term
.size() << endl
;
747 for (short_rat_list::iterator i
= gf
.term
.begin(); i
!= gf
.term
.end(); ++i
)
752 static Matrix
*Matrix_Read(std::istream
& is
)
759 M
= Matrix_Alloc(r
, c
);
760 for (int i
= 0; i
< r
; ++i
)
761 for (int j
= 0; j
< c
; ++j
) {
763 zz2value(tmp
, M
->p
[i
][j
]);
768 gen_fun
*gen_fun::read(std::istream
& is
, barvinok_options
*options
)
770 Matrix
*M
= Matrix_Read(is
);
771 Polyhedron
*C
= Constraints2Polyhedron(M
, options
->MaxRays
);
774 gen_fun
*gf
= new gen_fun(C
);
782 for (int i
= 0; i
< n
; ++i
) {
783 is
>> c
>> num
>> den
;
784 gf
->add(new short_rat(c
, num
, den
));
790 gen_fun::operator evalue
*() const
794 value_init(factor
.d
);
795 value_init(factor
.x
.n
);
796 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
797 unsigned nvar
= (*i
)->d
.power
.NumRows();
798 unsigned nparam
= (*i
)->d
.power
.NumCols();
799 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ nparam
+ 1);
800 mat_ZZ
& d
= (*i
)->d
.power
;
801 Polyhedron
*U
= context
? context
: Universe_Polyhedron(nparam
);
803 for (int j
= 0; j
< (*i
)->n
.coeff
.length(); ++j
) {
804 for (int r
= 0; r
< nparam
; ++r
) {
805 value_set_si(C
->p
[r
][0], 0);
806 for (int c
= 0; c
< nvar
; ++c
) {
807 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
809 Vector_Set(&C
->p
[r
][1+nvar
], 0, nparam
);
810 value_set_si(C
->p
[r
][1+nvar
+r
], -1);
811 zz2value((*i
)->n
.power
[j
][r
], C
->p
[r
][1+nvar
+nparam
]);
813 for (int r
= 0; r
< nvar
; ++r
) {
814 value_set_si(C
->p
[nparam
+r
][0], 1);
815 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ nparam
+ 1);
816 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
818 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
819 evalue
*E
= barvinok_enumerate_ev(P
, U
, 0);
821 if (EVALUE_IS_ZERO(*E
)) {
825 zz2value((*i
)->n
.coeff
[j
].n
, factor
.x
.n
);
826 zz2value((*i
)->n
.coeff
[j
].d
, factor
.d
);
829 Matrix_Print(stdout, P_VALUE_FMT, C);
830 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
831 print_evalue(stdout, E, test);
844 value_clear(factor
.d
);
845 value_clear(factor
.x
.n
);
849 ZZ
gen_fun::coefficient(Value
* params
, barvinok_options
*options
) const
851 if (context
&& !in_domain(context
, params
))
856 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
857 sum
+= (*i
)->coefficient(params
, options
);
863 void gen_fun::coefficient(Value
* params
, Value
* c
) const
865 barvinok_options
*options
= barvinok_options_new_with_defaults();
867 ZZ coeff
= coefficient(params
, options
);
871 barvinok_options_free(options
);
874 gen_fun
*gen_fun::summate(int nvar
, barvinok_options
*options
) const
876 int dim
= context
->Dimension
;
877 int nparam
= dim
- nvar
;
881 if (options
->incremental_specialization
== 1) {
882 red
= new partial_ireducer(Polyhedron_Project(context
, nparam
), dim
, nparam
);
884 red
= new partial_reducer(Polyhedron_Project(context
, nparam
), dim
, nparam
);
888 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
889 red
->reduce((*i
)->n
.coeff
, (*i
)->n
.power
, (*i
)->d
.power
);
891 } catch (OrthogonalException
&e
) {
900 /* returns true if the set was finite and false otherwise */
901 bool gen_fun::summate(Value
*sum
) const
903 if (term
.size() == 0) {
904 value_set_si(*sum
, 0);
909 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
910 if ((*i
)->d
.power
.NumRows() > maxlen
)
911 maxlen
= (*i
)->d
.power
.NumRows();
913 infinite_icounter
cnt((*term
.begin())->d
.power
.NumCols(), maxlen
);
914 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
915 cnt
.reduce((*i
)->n
.coeff
, (*i
)->n
.power
, (*i
)->d
.power
);
917 for (int i
= 1; i
<= maxlen
; ++i
)
918 if (value_notzero_p(mpq_numref(cnt
.count
[i
]))) {
919 value_set_si(*sum
, -1);
923 assert(value_one_p(mpq_denref(cnt
.count
[0])));
924 value_assign(*sum
, mpq_numref(cnt
.count
[0]));