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
)
254 for (short_rat_list::iterator i
= gf
->term
.begin(); i
!= gf
->term
.end(); ++i
) {
255 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
) {
257 p
*= (*i
)->n
.coeff
[j
];
258 add(p
, (*i
)->n
.power
[j
], (*i
)->d
.power
);
263 static void split_param_compression(Matrix
*CP
, mat_ZZ
& map
, vec_ZZ
& offset
)
265 Matrix
*T
= Transpose(CP
);
266 matrix2zz(T
, map
, T
->NbRows
-1, T
->NbColumns
-1);
267 values2zz(T
->p
[T
->NbRows
-1], offset
, T
->NbColumns
-1);
272 * Perform the substitution specified by CP
274 * CP is a homogeneous matrix that maps a set of "compressed parameters"
275 * to the original set of parameters.
277 * This function is applied to a gen_fun computed with the compressed parameters
278 * and adapts it to refer to the original parameters.
280 * That is, if y are the compressed parameters and x = A y + b are the original
281 * parameters, then we want the coefficient of the monomial t^y in the original
282 * generating function to be the coefficient of the monomial u^x in the resulting
283 * generating function.
284 * The original generating function has the form
286 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
288 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
290 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
292 * = a u^{A m + b}/(1-u^{A n})
294 * Therefore, we multiply the powers m and n in both numerator and denominator by A
295 * and add b to the power in the numerator.
296 * Since the above powers are stored as row vectors m^T and n^T,
297 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
299 * The pair (map, offset) contains the same information as CP.
300 * map is the transpose of the linear part of CP, while offset is the constant part.
302 void gen_fun::substitute(Matrix
*CP
)
306 split_param_compression(CP
, map
, offset
);
307 Polyhedron
*C
= Polyhedron_Image(context
, CP
, 0);
308 Polyhedron_Free(context
);
311 short_rat_list new_term
;
312 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
316 for (int j
= 0; j
< r
->n
.power
.NumRows(); ++j
)
317 r
->n
.power
[j
] += offset
;
324 struct parallel_cones
{
326 vector
<pair
<Vector
*, QQ
> > vertices
;
327 parallel_cones(int *pos
) : pos(pos
) {}
330 struct parallel_polytopes
{
337 vector
<parallel_cones
> cones
;
338 barvinok_options
*options
;
340 parallel_polytopes(int n
, Polyhedron
*context
, int nparam
,
341 barvinok_options
*options
) :
342 context(context
), dim(-1), nparam(nparam
),
349 bool add(const QQ
& c
, Polyhedron
*P
) {
352 for (i
= 0; i
< P
->NbEq
; ++i
)
353 if (First_Non_Zero(P
->Constraint
[i
]+1,
354 P
->Dimension
-nparam
) == -1)
359 Polyhedron
*Q
= remove_equalities_p(Polyhedron_Copy(P
), P
->Dimension
-nparam
,
360 NULL
, options
->MaxRays
);
361 POL_ENSURE_VERTICES(Q
);
371 M
= Matrix_Alloc(Q
->NbEq
, Q
->Dimension
+2);
372 Vector_Copy(Q
->Constraint
[0], M
->p
[0], Q
->NbEq
* (Q
->Dimension
+2));
373 CP
= compress_parms(M
, nparam
);
374 T
= align_matrix(CP
, Q
->Dimension
+1);
377 R
= Polyhedron_Preimage(Q
, T
, options
->MaxRays
);
379 Q
= remove_equalities_p(R
, R
->Dimension
-nparam
, NULL
,
382 assert(Q
->NbEq
== 0);
384 if (First_Non_Zero(Q
->Constraint
[Q
->NbConstraints
-1]+1, Q
->Dimension
) == -1)
389 red
= gf_base::create(Polyhedron_Copy(context
), dim
, nparam
, options
);
391 Constraints
= Matrix_Alloc(Q
->NbConstraints
, Q
->Dimension
);
392 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
393 Vector_Copy(Q
->Constraint
[i
]+1, Constraints
->p
[i
], Q
->Dimension
);
396 assert(Q
->Dimension
== dim
);
397 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
399 for (j
= 0; j
< Constraints
->NbRows
; ++j
)
400 if (Vector_Equal(Q
->Constraint
[i
]+1, Constraints
->p
[j
],
403 assert(j
< Constraints
->NbRows
);
407 for (int i
= 0; i
< Q
->NbRays
; ++i
) {
408 if (!value_pos_p(Q
->Ray
[i
][dim
+1]))
411 Polyhedron
*C
= supporting_cone(Q
, i
);
413 if (First_Non_Zero(C
->Constraint
[C
->NbConstraints
-1]+1,
417 int *pos
= new int[1+C
->NbConstraints
];
418 pos
[0] = C
->NbConstraints
;
420 for (int k
= 0; k
< Constraints
->NbRows
; ++k
) {
421 for (int j
= 0; j
< C
->NbConstraints
; ++j
) {
422 if (Vector_Equal(C
->Constraint
[j
]+1, Constraints
->p
[k
],
429 assert(l
== C
->NbConstraints
);
432 for (j
= 0; j
< cones
.size(); ++j
)
433 if (!memcmp(pos
, cones
[j
].pos
, (1+C
->NbConstraints
)*sizeof(int)))
435 if (j
== cones
.size())
436 cones
.push_back(parallel_cones(pos
));
443 for (k
= 0; k
< cones
[j
].vertices
.size(); ++k
)
444 if (Vector_Equal(Q
->Ray
[i
]+1, cones
[j
].vertices
[k
].first
->p
,
448 if (k
== cones
[j
].vertices
.size()) {
449 Vector
*vertex
= Vector_Alloc(Q
->Dimension
+1);
450 Vector_Copy(Q
->Ray
[i
]+1, vertex
->p
, Q
->Dimension
+1);
451 cones
[j
].vertices
.push_back(pair
<Vector
*,QQ
>(vertex
, c
));
453 cones
[j
].vertices
[k
].second
+= c
;
454 if (cones
[j
].vertices
[k
].second
.n
== 0) {
455 int size
= cones
[j
].vertices
.size();
456 Vector_Free(cones
[j
].vertices
[k
].first
);
458 cones
[j
].vertices
[k
] = cones
[j
].vertices
[size
-1];
459 cones
[j
].vertices
.pop_back();
470 for (int i
= 0; i
< cones
.size(); ++i
) {
471 Matrix
*M
= Matrix_Alloc(cones
[i
].pos
[0], 1+Constraints
->NbColumns
+1);
473 for (int j
= 0; j
<cones
[i
].pos
[0]; ++j
) {
474 value_set_si(M
->p
[j
][0], 1);
475 Vector_Copy(Constraints
->p
[cones
[i
].pos
[1+j
]], M
->p
[j
]+1,
476 Constraints
->NbColumns
);
478 Cone
= Constraints2Polyhedron(M
, options
->MaxRays
);
480 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
481 red
->base
->do_vertex_cone(cones
[i
].vertices
[j
].second
,
482 Polyhedron_Copy(Cone
),
483 cones
[i
].vertices
[j
].first
->p
, options
);
485 Polyhedron_Free(Cone
);
488 red
->gf
->substitute(CP
);
491 void print(std::ostream
& os
) const {
492 for (int i
= 0; i
< cones
.size(); ++i
) {
494 for (int j
= 0; j
< cones
[i
].pos
[0]; ++j
) {
497 os
<< cones
[i
].pos
[1+j
];
500 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
501 Vector_Print(stderr
, P_VALUE_FMT
, cones
[i
].vertices
[j
].first
);
502 os
<< cones
[i
].vertices
[j
].second
<< endl
;
506 ~parallel_polytopes() {
507 for (int i
= 0; i
< cones
.size(); ++i
) {
508 delete [] cones
[i
].pos
;
509 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
)
510 Vector_Free(cones
[i
].vertices
[j
].first
);
513 Matrix_Free(Constraints
);
522 gen_fun
*gen_fun::Hadamard_product(const gen_fun
*gf
, barvinok_options
*options
)
525 Polyhedron
*C
= DomainIntersection(context
, gf
->context
, options
->MaxRays
);
526 Polyhedron
*U
= Universe_Polyhedron(C
->Dimension
);
527 gen_fun
*sum
= new gen_fun(C
);
528 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
529 for (short_rat_list::iterator i2
= gf
->term
.begin(); i2
!= gf
->term
.end();
531 int d
= (*i
)->d
.power
.NumCols();
532 int k1
= (*i
)->d
.power
.NumRows();
533 int k2
= (*i2
)->d
.power
.NumRows();
534 assert((*i
)->d
.power
.NumCols() == (*i2
)->d
.power
.NumCols());
536 parallel_polytopes
pp((*i
)->n
.power
.NumRows() *
537 (*i2
)->n
.power
.NumRows(),
538 sum
->context
, d
, options
);
540 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
) {
541 for (int j2
= 0; j2
< (*i2
)->n
.power
.NumRows(); ++j2
) {
542 Matrix
*M
= Matrix_Alloc(k1
+k2
+d
+d
, 1+k1
+k2
+d
+1);
543 for (int k
= 0; k
< k1
+k2
; ++k
) {
544 value_set_si(M
->p
[k
][0], 1);
545 value_set_si(M
->p
[k
][1+k
], 1);
547 for (int k
= 0; k
< d
; ++k
) {
548 value_set_si(M
->p
[k1
+k2
+k
][1+k1
+k2
+k
], -1);
549 zz2value((*i
)->n
.power
[j
][k
], M
->p
[k1
+k2
+k
][1+k1
+k2
+d
]);
550 for (int l
= 0; l
< k1
; ++l
)
551 zz2value((*i
)->d
.power
[l
][k
], M
->p
[k1
+k2
+k
][1+l
]);
553 for (int k
= 0; k
< d
; ++k
) {
554 value_set_si(M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+k
], -1);
555 zz2value((*i2
)->n
.power
[j2
][k
],
556 M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+d
]);
557 for (int l
= 0; l
< k2
; ++l
)
558 zz2value((*i2
)->d
.power
[l
][k
],
559 M
->p
[k1
+k2
+d
+k
][1+k1
+l
]);
561 Polyhedron
*P
= Constraints2Polyhedron(M
, options
->MaxRays
);
564 QQ c
= (*i
)->n
.coeff
[j
];
565 c
*= (*i2
)->n
.coeff
[j2
];
567 gen_fun
*t
= barvinok_series_with_options(P
, U
, options
);
576 gen_fun
*t
= pp
.compute();
587 void gen_fun::add_union(gen_fun
*gf
, barvinok_options
*options
)
589 QQ
one(1, 1), mone(-1, 1);
591 gen_fun
*hp
= Hadamard_product(gf
, options
);
597 static void Polyhedron_Shift(Polyhedron
*P
, Vector
*offset
)
601 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
602 Inner_Product(P
->Constraint
[i
]+1, offset
->p
, P
->Dimension
, &tmp
);
603 value_subtract(P
->Constraint
[i
][1+P
->Dimension
],
604 P
->Constraint
[i
][1+P
->Dimension
], tmp
);
606 for (int i
= 0; i
< P
->NbRays
; ++i
) {
607 if (value_notone_p(P
->Ray
[i
][0]))
609 if (value_zero_p(P
->Ray
[i
][1+P
->Dimension
]))
611 Vector_Combine(P
->Ray
[i
]+1, offset
->p
, P
->Ray
[i
]+1,
612 P
->Ray
[i
][0], P
->Ray
[i
][1+P
->Dimension
], P
->Dimension
);
617 void gen_fun::shift(const vec_ZZ
& offset
)
619 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
620 for (int j
= 0; j
< (*i
)->n
.power
.NumRows(); ++j
)
621 (*i
)->n
.power
[j
] += offset
;
623 Vector
*v
= Vector_Alloc(offset
.length());
624 zz2values(offset
, v
->p
);
625 Polyhedron_Shift(context
, v
);
629 /* Divide the generating functin by 1/(1-z^power).
630 * The effect on the corresponding explicit function f(x) is
631 * f'(x) = \sum_{i=0}^\infty f(x - i * power)
633 void gen_fun::divide(const vec_ZZ
& power
)
635 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
636 int r
= (*i
)->d
.power
.NumRows();
637 int c
= (*i
)->d
.power
.NumCols();
638 (*i
)->d
.power
.SetDims(r
+1, c
);
639 (*i
)->d
.power
[r
] = power
;
642 Vector
*v
= Vector_Alloc(1+power
.length()+1);
643 value_set_si(v
->p
[0], 1);
644 zz2values(power
, v
->p
+1);
645 Polyhedron
*C
= AddRays(v
->p
, 1, context
, context
->NbConstraints
+1);
647 Polyhedron_Free(context
);
651 static void print_power(std::ostream
& os
, const QQ
& c
, const vec_ZZ
& p
,
652 unsigned int nparam
, char **param_name
)
656 for (int i
= 0; i
< p
.length(); ++i
) {
660 if (c
.n
== -1 && c
.d
== 1)
662 else if (c
.n
!= 1 || c
.d
!= 1) {
678 os
<< "^(" << p
[i
] << ")";
689 void short_rat::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
693 for (int j
= 0; j
< n
.coeff
.length(); ++j
) {
694 if (j
!= 0 && n
.coeff
[j
].n
>= 0)
696 print_power(os
, n
.coeff
[j
], n
.power
[j
], nparam
, param_name
);
699 for (int j
= 0; j
< d
.power
.NumRows(); ++j
) {
703 print_power(os
, mone
, d
.power
[j
], nparam
, param_name
);
709 void gen_fun::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
711 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
712 if (i
!= term
.begin())
714 (*i
)->print(os
, nparam
, param_name
);
718 std::ostream
& operator<< (std::ostream
& os
, const short_rat
& r
)
720 os
<< r
.n
.coeff
<< endl
;
721 os
<< r
.n
.power
<< endl
;
722 os
<< r
.d
.power
<< endl
;
726 std::ostream
& operator<< (std::ostream
& os
, const Polyhedron
& P
)
729 void (*gmp_free
)(void *, size_t);
730 mp_get_memory_functions(NULL
, NULL
, &gmp_free
);
731 os
<< P
.NbConstraints
<< " " << P
.Dimension
+2 << endl
;
732 for (int i
= 0; i
< P
.NbConstraints
; ++i
) {
733 for (int j
= 0; j
< P
.Dimension
+2; ++j
) {
734 str
= mpz_get_str(0, 10, P
.Constraint
[i
][j
]);
735 os
<< std::setw(4) << str
<< " ";
736 (*gmp_free
)(str
, strlen(str
)+1);
743 std::ostream
& operator<< (std::ostream
& os
, const gen_fun
& gf
)
745 os
<< *gf
.context
<< endl
;
747 os
<< gf
.term
.size() << endl
;
748 for (short_rat_list::iterator i
= gf
.term
.begin(); i
!= gf
.term
.end(); ++i
)
753 gen_fun
*gen_fun::read(std::istream
& is
, barvinok_options
*options
)
755 Matrix
*M
= Matrix_Read(is
);
756 Polyhedron
*C
= Constraints2Polyhedron(M
, options
->MaxRays
);
759 gen_fun
*gf
= new gen_fun(C
);
767 for (int i
= 0; i
< n
; ++i
) {
768 is
>> c
>> num
>> den
;
769 gf
->add(new short_rat(c
, num
, den
));
775 gen_fun::operator evalue
*() const
779 value_init(factor
.d
);
780 value_init(factor
.x
.n
);
781 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
) {
782 unsigned nvar
= (*i
)->d
.power
.NumRows();
783 unsigned nparam
= (*i
)->d
.power
.NumCols();
784 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ nparam
+ 1);
785 mat_ZZ
& d
= (*i
)->d
.power
;
786 Polyhedron
*U
= context
? context
: Universe_Polyhedron(nparam
);
788 for (int j
= 0; j
< (*i
)->n
.coeff
.length(); ++j
) {
789 for (int r
= 0; r
< nparam
; ++r
) {
790 value_set_si(C
->p
[r
][0], 0);
791 for (int c
= 0; c
< nvar
; ++c
) {
792 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
794 Vector_Set(&C
->p
[r
][1+nvar
], 0, nparam
);
795 value_set_si(C
->p
[r
][1+nvar
+r
], -1);
796 zz2value((*i
)->n
.power
[j
][r
], C
->p
[r
][1+nvar
+nparam
]);
798 for (int r
= 0; r
< nvar
; ++r
) {
799 value_set_si(C
->p
[nparam
+r
][0], 1);
800 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ nparam
+ 1);
801 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
803 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
804 evalue
*E
= barvinok_enumerate_ev(P
, U
, 0);
806 if (EVALUE_IS_ZERO(*E
)) {
810 zz2value((*i
)->n
.coeff
[j
].n
, factor
.x
.n
);
811 zz2value((*i
)->n
.coeff
[j
].d
, factor
.d
);
814 Matrix_Print(stdout, P_VALUE_FMT, C);
815 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
816 print_evalue(stdout, E, test);
829 value_clear(factor
.d
);
830 value_clear(factor
.x
.n
);
834 ZZ
gen_fun::coefficient(Value
* params
, barvinok_options
*options
) const
836 if (context
&& !in_domain(context
, params
))
841 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
842 sum
+= (*i
)->coefficient(params
, options
);
848 void gen_fun::coefficient(Value
* params
, Value
* c
) const
850 barvinok_options
*options
= barvinok_options_new_with_defaults();
852 ZZ coeff
= coefficient(params
, options
);
856 barvinok_options_free(options
);
859 gen_fun
*gen_fun::summate(int nvar
, barvinok_options
*options
) const
861 int dim
= context
->Dimension
;
862 int nparam
= dim
- nvar
;
866 if (options
->incremental_specialization
== 1) {
867 red
= new partial_ireducer(Polyhedron_Project(context
, nparam
), dim
, nparam
);
869 red
= new partial_reducer(Polyhedron_Project(context
, nparam
), dim
, nparam
);
873 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
874 red
->reduce((*i
)->n
.coeff
, (*i
)->n
.power
, (*i
)->d
.power
);
876 } catch (OrthogonalException
&e
) {
885 /* returns true if the set was finite and false otherwise */
886 bool gen_fun::summate(Value
*sum
) const
888 if (term
.size() == 0) {
889 value_set_si(*sum
, 0);
894 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
895 if ((*i
)->d
.power
.NumRows() > maxlen
)
896 maxlen
= (*i
)->d
.power
.NumRows();
898 infinite_icounter
cnt((*term
.begin())->d
.power
.NumCols(), maxlen
);
899 for (short_rat_list::iterator i
= term
.begin(); i
!= term
.end(); ++i
)
900 cnt
.reduce((*i
)->n
.coeff
, (*i
)->n
.power
, (*i
)->d
.power
);
902 for (int i
= 1; i
<= maxlen
; ++i
)
903 if (value_notzero_p(mpq_numref(cnt
.count
[i
]))) {
904 value_set_si(*sum
, -1);
908 assert(value_one_p(mpq_denref(cnt
.count
[0])));
909 value_assign(*sum
, mpq_numref(cnt
.count
[0]));