5 #include <barvinok/genfun.h>
6 #include <barvinok/barvinok.h>
7 #include "conversion.h"
8 #include "genfun_constructor.h"
10 #include "matrix_read.h"
18 bool short_rat_lex_smaller_denominator::operator()(const short_rat
* r1
,
19 const short_rat
* r2
) const
21 return lex_cmp(r1
->d
.power
, r2
->d
.power
) < 0;
24 static void lex_order_terms(struct short_rat
* rat
)
26 for (int i
= 0; i
< rat
->n
.power
.NumRows(); ++i
) {
28 for (int j
= i
+1; j
< rat
->n
.power
.NumRows(); ++j
)
29 if (lex_cmp(rat
->n
.power
[j
], rat
->n
.power
[m
]) < 0)
32 vec_ZZ tmp
= rat
->n
.power
[m
];
33 rat
->n
.power
[m
] = rat
->n
.power
[i
];
34 rat
->n
.power
[i
] = tmp
;
35 QQ tmp_coeff
= rat
->n
.coeff
[m
];
36 rat
->n
.coeff
[m
] = rat
->n
.coeff
[i
];
37 rat
->n
.coeff
[i
] = tmp_coeff
;
42 short_rat::short_rat(const short_rat
& r
)
49 short_rat::short_rat(Value c
)
52 value2zz(c
, n
.coeff
[0].n
);
54 n
.power
.SetDims(1, 0);
55 d
.power
.SetDims(0, 0);
58 short_rat::short_rat(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den
)
64 n
.power
.SetDims(1, num
.length());
70 short_rat::short_rat(const vec_QQ
& c
, const mat_ZZ
& num
, const mat_ZZ
& den
)
78 void short_rat::normalize()
80 /* Make all powers in denominator reverse-lexico-positive */
81 for (int i
= 0; i
< d
.power
.NumRows(); ++i
) {
83 for (j
= d
.power
.NumCols()-1; j
>= 0; --j
)
84 if (!IsZero(d
.power
[i
][j
]))
87 if (sign(d
.power
[i
][j
]) < 0) {
88 negate(d
.power
[i
], d
.power
[i
]);
89 for (int k
= 0; k
< n
.coeff
.length(); ++k
) {
90 negate(n
.coeff
[k
].n
, n
.coeff
[k
].n
);
91 n
.power
[k
] += d
.power
[i
];
96 /* Order powers in denominator */
97 lex_order_rows(d
.power
);
100 void short_rat::add(const short_rat
*r
)
102 for (int i
= 0; i
< r
->n
.power
.NumRows(); ++i
) {
103 int len
= n
.coeff
.length();
105 for (j
= 0; j
< len
; ++j
)
106 if (r
->n
.power
[i
] == n
.power
[j
])
109 n
.coeff
[j
] += r
->n
.coeff
[i
];
110 if (n
.coeff
[j
].n
== 0) {
112 n
.power
[j
] = n
.power
[len
-1];
113 n
.coeff
[j
] = n
.coeff
[len
-1];
115 int dim
= n
.power
.NumCols();
116 n
.coeff
.SetLength(len
-1);
117 n
.power
.SetDims(len
-1, dim
);
120 int dim
= n
.power
.NumCols();
121 n
.coeff
.SetLength(len
+1);
122 n
.power
.SetDims(len
+1, dim
);
123 n
.coeff
[len
] = r
->n
.coeff
[i
];
124 n
.power
[len
] = r
->n
.power
[i
];
129 QQ
short_rat::coefficient(Value
* params
, barvinok_options
*options
) const
131 unsigned nvar
= d
.power
.NumRows();
132 unsigned nparam
= d
.power
.NumCols();
133 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ 1);
139 for (int j
= 0; j
< n
.coeff
.length(); ++j
) {
140 C
->NbRows
= nparam
+nvar
;
141 for (int r
= 0; r
< nparam
; ++r
) {
142 value_set_si(C
->p
[r
][0], 0);
143 for (int c
= 0; c
< nvar
; ++c
) {
144 zz2value(d
.power
[c
][r
], C
->p
[r
][1+c
]);
146 zz2value(n
.power
[j
][r
], C
->p
[r
][1+nvar
]);
147 value_subtract(C
->p
[r
][1+nvar
], C
->p
[r
][1+nvar
], params
[r
]);
149 for (int r
= 0; r
< nvar
; ++r
) {
150 value_set_si(C
->p
[nparam
+r
][0], 1);
151 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ 1);
152 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
154 Polyhedron
*P
= Constraints2Polyhedron(C
, options
->MaxRays
);
159 barvinok_count_with_options(P
, &tmp
, options
);
161 if (value_zero_p(tmp
))
173 bool short_rat::reduced()
175 int dim
= n
.power
.NumCols();
176 lex_order_terms(this);
177 if (n
.power
.NumRows() % 2 == 0) {
178 if (n
.coeff
[0].n
== -n
.coeff
[1].n
&&
179 n
.coeff
[0].d
== n
.coeff
[1].d
) {
180 vec_ZZ step
= n
.power
[1] - n
.power
[0];
182 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
183 if (n
.coeff
[2*k
].n
!= -n
.coeff
[2*k
+1].n
||
184 n
.coeff
[2*k
].d
!= n
.coeff
[2*k
+1].d
)
186 if (step
!= n
.power
[2*k
+1] - n
.power
[2*k
])
189 if (k
== n
.power
.NumRows()/2) {
190 for (k
= 0; k
< d
.power
.NumRows(); ++k
)
191 if (d
.power
[k
] == step
)
193 if (k
< d
.power
.NumRows()) {
194 for (++k
; k
< d
.power
.NumRows(); ++k
)
195 d
.power
[k
-1] = d
.power
[k
];
196 d
.power
.SetDims(k
-1, dim
);
197 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
198 n
.coeff
[k
] = n
.coeff
[2*k
];
199 n
.power
[k
] = n
.power
[2*k
];
201 n
.coeff
.SetLength(k
);
202 n
.power
.SetDims(k
, dim
);
211 gen_fun::gen_fun(Value c
)
213 short_rat
*r
= new short_rat(c
);
214 context
= Universe_Polyhedron(0);
218 void gen_fun::add(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den
)
223 add(new short_rat(c
, num
, den
));
226 void gen_fun::add(short_rat
*r
)
228 short_rat_list::iterator i
= term
.find(r
);
229 while (i
!= term
.end()) {
231 if ((*i
)->n
.coeff
.length() == 0) {
234 } else if ((*i
)->reduced()) {
236 /* we've modified term[i], so remove it
237 * and add it back again
251 void gen_fun::add(const QQ
& c
, const gen_fun
*gf
, barvinok_options
*options
)
253 Polyhedron
*U
= DomainUnion(context
, gf
->context
, options
->MaxRays
);
254 Polyhedron
*C
= DomainConvex(U
, options
->MaxRays
);
256 Domain_Free(context
);
262 void gen_fun::add(const QQ
& c
, const gen_fun
*gf
)
265 for (short_rat_list::iterator i
= gf
->term
.begin(); i
!= gf
->term
.end(); ++i
) {
266 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
) {
268 p
*= (*i
)->n
.coeff
[j
];
269 add(p
, (*i
)->n
.power
[j
], (*i
)->d
.power
);
274 static void split_param_compression(Matrix
*CP
, mat_ZZ
& map
, vec_ZZ
& offset
)
276 Matrix
*T
= Transpose(CP
);
277 matrix2zz(T
, map
, T
->NbRows
-1, T
->NbColumns
-1);
278 values2zz(T
->p
[T
->NbRows
-1], offset
, T
->NbColumns
-1);
283 * Perform the substitution specified by CP
285 * CP is a homogeneous matrix that maps a set of "compressed parameters"
286 * to the original set of parameters.
288 * This function is applied to a gen_fun computed with the compressed parameters
289 * and adapts it to refer to the original parameters.
291 * That is, if y are the compressed parameters and x = A y + b are the original
292 * parameters, then we want the coefficient of the monomial t^y in the original
293 * generating function to be the coefficient of the monomial u^x in the resulting
294 * generating function.
295 * The original generating function has the form
297 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
299 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
301 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
303 * = a u^{A m + b}/(1-u^{A n})
305 * Therefore, we multiply the powers m and n in both numerator and denominator by A
306 * and add b to the power in the numerator.
307 * Since the above powers are stored as row vectors m^T and n^T,
308 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
310 * The pair (map, offset) contains the same information as CP.
311 * map is the transpose of the linear part of CP, while offset is the constant part.
313 void gen_fun::substitute(Matrix
*CP
)
317 split_param_compression(CP
, map
, offset
);
318 Polyhedron
*C
= Polyhedron_Image(context
, CP
, 0);
319 Polyhedron_Free(context
);
322 short_rat_list new_term
;
323 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
327 for (int j
= 0; j
< r
->n
.power
.NumRows(); ++j
)
328 r
->n
.power
[j
] += offset
;
335 struct parallel_cones
{
337 vector
<pair
<Vector
*, QQ
> > vertices
;
338 parallel_cones(int *pos
) : pos(pos
) {}
341 struct parallel_polytopes
{
348 vector
<parallel_cones
> cones
;
349 barvinok_options
*options
;
351 parallel_polytopes(int n
, Polyhedron
*context
, int nparam
,
352 barvinok_options
*options
) :
353 context(context
), dim(-1), nparam(nparam
),
360 bool add(const QQ
& c
, Polyhedron
*P
) {
363 for (i
= 0; i
< P
->NbEq
; ++i
)
364 if (First_Non_Zero(P
->Constraint
[i
]+1,
365 P
->Dimension
-nparam
) == -1)
370 Polyhedron
*Q
= remove_equalities_p(Polyhedron_Copy(P
), P
->Dimension
-nparam
,
371 NULL
, options
->MaxRays
);
372 POL_ENSURE_VERTICES(Q
);
382 M
= Matrix_Alloc(Q
->NbEq
, Q
->Dimension
+2);
383 Vector_Copy(Q
->Constraint
[0], M
->p
[0], Q
->NbEq
* (Q
->Dimension
+2));
384 CP
= compress_parms(M
, nparam
);
385 T
= align_matrix(CP
, Q
->Dimension
+1);
388 R
= Polyhedron_Preimage(Q
, T
, options
->MaxRays
);
390 Q
= remove_equalities_p(R
, R
->Dimension
-nparam
, NULL
,
393 assert(Q
->NbEq
== 0);
395 if (First_Non_Zero(Q
->Constraint
[Q
->NbConstraints
-1]+1, Q
->Dimension
) == -1)
400 red
= gf_base::create(Polyhedron_Copy(context
), dim
, nparam
, options
);
402 Constraints
= Matrix_Alloc(Q
->NbConstraints
, Q
->Dimension
);
403 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
404 Vector_Copy(Q
->Constraint
[i
]+1, Constraints
->p
[i
], Q
->Dimension
);
407 assert(Q
->Dimension
== dim
);
408 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
410 for (j
= 0; j
< Constraints
->NbRows
; ++j
)
411 if (Vector_Equal(Q
->Constraint
[i
]+1, Constraints
->p
[j
],
414 assert(j
< Constraints
->NbRows
);
418 for (int i
= 0; i
< Q
->NbRays
; ++i
) {
419 if (!value_pos_p(Q
->Ray
[i
][dim
+1]))
422 Polyhedron
*C
= supporting_cone(Q
, i
);
424 if (First_Non_Zero(C
->Constraint
[C
->NbConstraints
-1]+1,
428 int *pos
= new int[1+C
->NbConstraints
];
429 pos
[0] = C
->NbConstraints
;
431 for (int k
= 0; k
< Constraints
->NbRows
; ++k
) {
432 for (int j
= 0; j
< C
->NbConstraints
; ++j
) {
433 if (Vector_Equal(C
->Constraint
[j
]+1, Constraints
->p
[k
],
440 assert(l
== C
->NbConstraints
);
443 for (j
= 0; j
< cones
.size(); ++j
)
444 if (!memcmp(pos
, cones
[j
].pos
, (1+C
->NbConstraints
)*sizeof(int)))
446 if (j
== cones
.size())
447 cones
.push_back(parallel_cones(pos
));
454 for (k
= 0; k
< cones
[j
].vertices
.size(); ++k
)
455 if (Vector_Equal(Q
->Ray
[i
]+1, cones
[j
].vertices
[k
].first
->p
,
459 if (k
== cones
[j
].vertices
.size()) {
460 Vector
*vertex
= Vector_Alloc(Q
->Dimension
+1);
461 Vector_Copy(Q
->Ray
[i
]+1, vertex
->p
, Q
->Dimension
+1);
462 cones
[j
].vertices
.push_back(pair
<Vector
*,QQ
>(vertex
, c
));
464 cones
[j
].vertices
[k
].second
+= c
;
465 if (cones
[j
].vertices
[k
].second
.n
== 0) {
466 int size
= cones
[j
].vertices
.size();
467 Vector_Free(cones
[j
].vertices
[k
].first
);
469 cones
[j
].vertices
[k
] = cones
[j
].vertices
[size
-1];
470 cones
[j
].vertices
.pop_back();
481 for (int i
= 0; i
< cones
.size(); ++i
) {
482 Matrix
*M
= Matrix_Alloc(cones
[i
].pos
[0], 1+Constraints
->NbColumns
+1);
484 for (int j
= 0; j
<cones
[i
].pos
[0]; ++j
) {
485 value_set_si(M
->p
[j
][0], 1);
486 Vector_Copy(Constraints
->p
[cones
[i
].pos
[1+j
]], M
->p
[j
]+1,
487 Constraints
->NbColumns
);
489 Cone
= Constraints2Polyhedron(M
, options
->MaxRays
);
491 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
492 red
->base
->do_vertex_cone(cones
[i
].vertices
[j
].second
,
493 Polyhedron_Copy(Cone
),
494 cones
[i
].vertices
[j
].first
->p
, options
);
496 Polyhedron_Free(Cone
);
499 red
->gf
->substitute(CP
);
502 void print(std::ostream
& os
) const {
503 for (int i
= 0; i
< cones
.size(); ++i
) {
505 for (int j
= 0; j
< cones
[i
].pos
[0]; ++j
) {
508 os
<< cones
[i
].pos
[1+j
];
511 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
512 Vector_Print(stderr
, P_VALUE_FMT
, cones
[i
].vertices
[j
].first
);
513 os
<< cones
[i
].vertices
[j
].second
<< endl
;
517 ~parallel_polytopes() {
518 for (int i
= 0; i
< cones
.size(); ++i
) {
519 delete [] cones
[i
].pos
;
520 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
)
521 Vector_Free(cones
[i
].vertices
[j
].first
);
524 Matrix_Free(Constraints
);
533 gen_fun
*gen_fun::Hadamard_product(const gen_fun
*gf
, barvinok_options
*options
)
536 Polyhedron
*C
= DomainIntersection(context
, gf
->context
, options
->MaxRays
);
537 Polyhedron
*U
= Universe_Polyhedron(C
->Dimension
);
538 gen_fun
*sum
= new gen_fun(C
);
539 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
540 for (short_rat_list::iterator i2
= gf
->term
.begin(); i2
!= gf
->term
.end();
542 int d
= (*i
)->d
.power
.NumCols();
543 int k1
= (*i
)->d
.power
.NumRows();
544 int k2
= (*i2
)->d
.power
.NumRows();
545 assert((*i
)->d
.power
.NumCols() == (*i2
)->d
.power
.NumCols());
547 parallel_polytopes
pp((*i
)->n
.power
.NumRows() *
548 (*i2
)->n
.power
.NumRows(),
549 sum
->context
, d
, options
);
551 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
) {
552 for (int j2
= 0; j2
< (*i2
)->n
.power
.NumRows(); ++j2
) {
553 Matrix
*M
= Matrix_Alloc(k1
+k2
+d
+d
, 1+k1
+k2
+d
+1);
554 for (int k
= 0; k
< k1
+k2
; ++k
) {
555 value_set_si(M
->p
[k
][0], 1);
556 value_set_si(M
->p
[k
][1+k
], 1);
558 for (int k
= 0; k
< d
; ++k
) {
559 value_set_si(M
->p
[k1
+k2
+k
][1+k1
+k2
+k
], -1);
560 zz2value((*i
)->n
.power
[j
][k
], M
->p
[k1
+k2
+k
][1+k1
+k2
+d
]);
561 for (int l
= 0; l
< k1
; ++l
)
562 zz2value((*i
)->d
.power
[l
][k
], M
->p
[k1
+k2
+k
][1+l
]);
564 for (int k
= 0; k
< d
; ++k
) {
565 value_set_si(M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+k
], -1);
566 zz2value((*i2
)->n
.power
[j2
][k
],
567 M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+d
]);
568 for (int l
= 0; l
< k2
; ++l
)
569 zz2value((*i2
)->d
.power
[l
][k
],
570 M
->p
[k1
+k2
+d
+k
][1+k1
+l
]);
572 Polyhedron
*P
= Constraints2Polyhedron(M
, options
->MaxRays
);
575 QQ c
= (*i
)->n
.coeff
[j
];
576 c
*= (*i2
)->n
.coeff
[j2
];
578 gen_fun
*t
= barvinok_series_with_options(P
, U
, options
);
587 gen_fun
*t
= pp
.compute();
598 void gen_fun::add_union(gen_fun
*gf
, barvinok_options
*options
)
600 QQ
one(1, 1), mone(-1, 1);
602 gen_fun
*hp
= Hadamard_product(gf
, options
);
608 static void Polyhedron_Shift(Polyhedron
*P
, Vector
*offset
)
612 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
613 Inner_Product(P
->Constraint
[i
]+1, offset
->p
, P
->Dimension
, &tmp
);
614 value_subtract(P
->Constraint
[i
][1+P
->Dimension
],
615 P
->Constraint
[i
][1+P
->Dimension
], tmp
);
617 for (int i
= 0; i
< P
->NbRays
; ++i
) {
618 if (value_notone_p(P
->Ray
[i
][0]))
620 if (value_zero_p(P
->Ray
[i
][1+P
->Dimension
]))
622 Vector_Combine(P
->Ray
[i
]+1, offset
->p
, P
->Ray
[i
]+1,
623 P
->Ray
[i
][0], P
->Ray
[i
][1+P
->Dimension
], P
->Dimension
);
628 void gen_fun::shift(const vec_ZZ
& offset
)
630 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
631 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
)
632 (*i
)->n
.power
[j
] += offset
;
634 Vector
*v
= Vector_Alloc(offset
.length());
635 zz2values(offset
, v
->p
);
636 Polyhedron_Shift(context
, v
);
640 /* Divide the generating functin by 1/(1-z^power).
641 * The effect on the corresponding explicit function f(x) is
642 * f'(x) = \sum_{i=0}^\infty f(x - i * power)
644 void gen_fun::divide(const vec_ZZ
& power
)
646 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
647 int r
= (*i
)->d
.power
.NumRows();
648 int c
= (*i
)->d
.power
.NumCols();
649 (*i
)->d
.power
.SetDims(r
+1, c
);
650 (*i
)->d
.power
[r
] = power
;
653 Vector
*v
= Vector_Alloc(1+power
.length()+1);
654 value_set_si(v
->p
[0], 1);
655 zz2values(power
, v
->p
+1);
656 Polyhedron
*C
= AddRays(v
->p
, 1, context
, context
->NbConstraints
+1);
658 Polyhedron_Free(context
);
662 static void print_power(std::ostream
& os
, const QQ
& c
, const vec_ZZ
& p
,
663 unsigned int nparam
, char **param_name
)
667 for (int i
= 0; i
< p
.length(); ++i
) {
671 if (c
.n
== -1 && c
.d
== 1)
673 else if (c
.n
!= 1 || c
.d
!= 1) {
689 os
<< "^(" << p
[i
] << ")";
700 void short_rat::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
704 for (int j
= 0; j
< n
.coeff
.length(); ++j
) {
705 if (j
!= 0 && n
.coeff
[j
].n
>= 0)
707 print_power(os
, n
.coeff
[j
], n
.power
[j
], nparam
, param_name
);
710 for (int j
= 0; j
< d
.power
.NumRows(); ++j
) {
714 print_power(os
, mone
, d
.power
[j
], nparam
, param_name
);
720 void gen_fun::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
722 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
723 if (i
!= term
.begin())
725 (*i
)->print(os
, nparam
, param_name
);
729 std::ostream
& operator<< (std::ostream
& os
, const short_rat
& r
)
731 os
<< r
.n
.coeff
<< endl
;
732 os
<< r
.n
.power
<< endl
;
733 os
<< r
.d
.power
<< endl
;
737 std::ostream
& operator<< (std::ostream
& os
, const Polyhedron
& P
)
740 void (*gmp_free
)(void *, size_t);
741 mp_get_memory_functions(NULL
, NULL
, &gmp_free
);
742 os
<< P
.NbConstraints
<< " " << P
.Dimension
+2 << endl
;
743 for (int i
= 0; i
< P
.NbConstraints
; ++i
) {
744 for (int j
= 0; j
< P
.Dimension
+2; ++j
) {
745 str
= mpz_get_str(0, 10, P
.Constraint
[i
][j
]);
746 os
<< std::setw(4) << str
<< " ";
747 (*gmp_free
)(str
, strlen(str
)+1);
754 std::ostream
& operator<< (std::ostream
& os
, const gen_fun
& gf
)
756 os
<< *gf
.context
<< endl
;
758 os
<< gf
.term
.size() << endl
;
759 for (short_rat_list::iterator i
= gf
.term
.begin(); i
!= gf
.term
.end(); ++i
)
764 gen_fun
*gen_fun::read(std::istream
& is
, barvinok_options
*options
)
766 Matrix
*M
= Matrix_Read(is
);
767 Polyhedron
*C
= Constraints2Polyhedron(M
, options
->MaxRays
);
770 gen_fun
*gf
= new gen_fun(C
);
778 for (int i
= 0; i
< n
; ++i
) {
779 is
>> c
>> num
>> den
;
780 gf
->add(new short_rat(c
, num
, den
));
786 gen_fun::operator evalue
*() const
790 value_init(factor
.d
);
791 value_init(factor
.x
.n
);
792 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
793 unsigned nvar
= (*i
)->d
.power
.NumRows();
794 unsigned nparam
= (*i
)->d
.power
.NumCols();
795 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ nparam
+ 1);
796 mat_ZZ
& d
= (*i
)->d
.power
;
797 Polyhedron
*U
= context
;
799 for (int j
= 0; j
< (*i
)->n
.coeff
.length(); ++j
) {
800 for (int r
= 0; r
< nparam
; ++r
) {
801 value_set_si(C
->p
[r
][0], 0);
802 for (int c
= 0; c
< nvar
; ++c
) {
803 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
805 Vector_Set(&C
->p
[r
][1+nvar
], 0, nparam
);
806 value_set_si(C
->p
[r
][1+nvar
+r
], -1);
807 zz2value((*i
)->n
.power
[j
][r
], C
->p
[r
][1+nvar
+nparam
]);
809 for (int r
= 0; r
< nvar
; ++r
) {
810 value_set_si(C
->p
[nparam
+r
][0], 1);
811 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ nparam
+ 1);
812 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
814 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
815 evalue
*E
= barvinok_enumerate_ev(P
, U
, 0);
817 if (EVALUE_IS_ZERO(*E
)) {
821 zz2value((*i
)->n
.coeff
[j
].n
, factor
.x
.n
);
822 zz2value((*i
)->n
.coeff
[j
].d
, factor
.d
);
833 value_clear(factor
.d
);
834 value_clear(factor
.x
.n
);
835 return EP
? EP
: evalue_zero();
838 ZZ
gen_fun::coefficient(Value
* params
, barvinok_options
*options
) const
840 if (!in_domain(context
, params
))
845 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
846 sum
+= (*i
)->coefficient(params
, options
);
852 void gen_fun::coefficient(Value
* params
, Value
* c
) const
854 barvinok_options
*options
= barvinok_options_new_with_defaults();
856 ZZ coeff
= coefficient(params
, options
);
860 barvinok_options_free(options
);
863 gen_fun
*gen_fun::summate(int nvar
, barvinok_options
*options
) const
865 int dim
= context
->Dimension
;
866 int nparam
= dim
- nvar
;
870 if (options
->incremental_specialization
== 1) {
871 red
= new partial_ireducer(Polyhedron_Project(context
, nparam
), dim
, nparam
);
873 red
= new partial_reducer(Polyhedron_Project(context
, nparam
), dim
, nparam
);
877 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
878 red
->reduce((*i
)->n
.coeff
, (*i
)->n
.power
, (*i
)->d
.power
);
880 } catch (OrthogonalException
&e
) {
889 /* returns true if the set was finite and false otherwise */
890 bool gen_fun::summate(Value
*sum
) const
892 if (term
.size() == 0) {
893 value_set_si(*sum
, 0);
898 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
899 if ((*i
)->d
.power
.NumRows() > maxlen
)
900 maxlen
= (*i
)->d
.power
.NumRows();
902 infinite_icounter
cnt((*term
.begin())->d
.power
.NumCols(), maxlen
);
903 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
904 cnt
.reduce((*i
)->n
.coeff
, (*i
)->n
.power
, (*i
)->d
.power
);
906 for (int i
= 1; i
<= maxlen
; ++i
)
907 if (value_notzero_p(mpq_numref(cnt
.count
[i
]))) {
908 value_set_si(*sum
, -1);
912 assert(value_one_p(mpq_denref(cnt
.count
[0])));
913 value_assign(*sum
, mpq_numref(cnt
.count
[0]));