5 #include <barvinok/genfun.h>
6 #include <barvinok/barvinok.h>
7 #include "conversion.h"
8 #include "genfun_constructor.h"
17 static int lex_cmp(mat_ZZ
& a
, mat_ZZ
& b
)
19 assert(a
.NumCols() == b
.NumCols());
20 int alen
= a
.NumRows();
21 int blen
= b
.NumRows();
22 int len
= alen
< blen
? alen
: blen
;
24 for (int i
= 0; i
< len
; ++i
) {
25 int s
= lex_cmp(a
[i
], b
[i
]);
32 static void lex_order_terms(struct short_rat
* rat
)
34 for (int i
= 0; i
< rat
->n
.power
.NumRows(); ++i
) {
36 for (int j
= i
+1; j
< rat
->n
.power
.NumRows(); ++j
)
37 if (lex_cmp(rat
->n
.power
[j
], rat
->n
.power
[m
]) < 0)
40 vec_ZZ tmp
= rat
->n
.power
[m
];
41 rat
->n
.power
[m
] = rat
->n
.power
[i
];
42 rat
->n
.power
[i
] = tmp
;
43 QQ tmp_coeff
= rat
->n
.coeff
[m
];
44 rat
->n
.coeff
[m
] = rat
->n
.coeff
[i
];
45 rat
->n
.coeff
[i
] = tmp_coeff
;
50 void short_rat::add(short_rat
*r
)
52 for (int i
= 0; i
< r
->n
.power
.NumRows(); ++i
) {
53 int len
= n
.coeff
.length();
55 for (j
= 0; j
< len
; ++j
)
56 if (r
->n
.power
[i
] == n
.power
[j
])
59 n
.coeff
[j
] += r
->n
.coeff
[i
];
60 if (n
.coeff
[j
].n
== 0) {
62 n
.power
[j
] = n
.power
[len
-1];
63 n
.coeff
[j
] = n
.coeff
[len
-1];
65 int dim
= n
.power
.NumCols();
66 n
.coeff
.SetLength(len
-1);
67 n
.power
.SetDims(len
-1, dim
);
70 int dim
= n
.power
.NumCols();
71 n
.coeff
.SetLength(len
+1);
72 n
.power
.SetDims(len
+1, dim
);
73 n
.coeff
[len
] = r
->n
.coeff
[i
];
74 n
.power
[len
] = r
->n
.power
[i
];
79 bool short_rat::reduced()
81 int dim
= n
.power
.NumCols();
82 lex_order_terms(this);
83 if (n
.power
.NumRows() % 2 == 0) {
84 if (n
.coeff
[0].n
== -n
.coeff
[1].n
&&
85 n
.coeff
[0].d
== n
.coeff
[1].d
) {
86 vec_ZZ step
= n
.power
[1] - n
.power
[0];
88 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
89 if (n
.coeff
[2*k
].n
!= -n
.coeff
[2*k
+1].n
||
90 n
.coeff
[2*k
].d
!= n
.coeff
[2*k
+1].d
)
92 if (step
!= n
.power
[2*k
+1] - n
.power
[2*k
])
95 if (k
== n
.power
.NumRows()/2) {
96 for (k
= 0; k
< d
.power
.NumRows(); ++k
)
97 if (d
.power
[k
] == step
)
99 if (k
< d
.power
.NumRows()) {
100 for (++k
; k
< d
.power
.NumRows(); ++k
)
101 d
.power
[k
-1] = d
.power
[k
];
102 d
.power
.SetDims(k
-1, dim
);
103 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
104 n
.coeff
[k
] = n
.coeff
[2*k
];
105 n
.power
[k
] = n
.power
[2*k
];
107 n
.coeff
.SetLength(k
);
108 n
.power
.SetDims(k
, dim
);
117 void gen_fun::add(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den
)
122 short_rat
* r
= new short_rat
;
123 r
->n
.coeff
.SetLength(1);
124 ZZ g
= GCD(c
.n
, c
.d
);
125 r
->n
.coeff
[0].n
= c
.n
/g
;
126 r
->n
.coeff
[0].d
= c
.d
/g
;
127 r
->n
.power
.SetDims(1, num
.length());
131 /* Make all powers in denominator lexico-positive */
132 for (int i
= 0; i
< r
->d
.power
.NumRows(); ++i
) {
134 for (j
= 0; j
< r
->d
.power
.NumCols(); ++j
)
135 if (r
->d
.power
[i
][j
] != 0)
137 if (r
->d
.power
[i
][j
] < 0) {
138 r
->d
.power
[i
] = -r
->d
.power
[i
];
139 r
->n
.coeff
[0].n
= -r
->n
.coeff
[0].n
;
140 r
->n
.power
[0] += r
->d
.power
[i
];
144 /* Order powers in denominator */
145 lex_order_rows(r
->d
.power
);
147 for (int i
= 0; i
< term
.size(); ++i
)
148 if (lex_cmp(term
[i
]->d
.power
, r
->d
.power
) == 0) {
150 if (term
[i
]->n
.coeff
.length() == 0) {
152 if (i
!= term
.size()-1)
153 term
[i
] = term
[term
.size()-1];
155 } else if (term
[i
]->reduced()) {
157 /* we've modified term[i], so removed it
158 * and add it back again
161 if (i
!= term
.size()-1)
162 term
[i
] = term
[term
.size()-1];
174 void gen_fun::add(const QQ
& c
, const gen_fun
*gf
)
177 for (int i
= 0; i
< gf
->term
.size(); ++i
) {
178 for (int j
= 0; j
< gf
->term
[i
]->n
.power
.NumRows(); ++j
) {
180 p
*= gf
->term
[i
]->n
.coeff
[j
];
181 add(p
, gf
->term
[i
]->n
.power
[j
], gf
->term
[i
]->d
.power
);
186 static void split_param_compression(Matrix
*CP
, mat_ZZ
& map
, vec_ZZ
& offset
)
188 Matrix
*T
= Transpose(CP
);
189 matrix2zz(T
, map
, T
->NbRows
-1, T
->NbColumns
-1);
190 values2zz(T
->p
[T
->NbRows
-1], offset
, T
->NbColumns
-1);
195 * Perform the substitution specified by CP
197 * CP is a homogeneous matrix that maps a set of "compressed parameters"
198 * to the original set of parameters.
200 * This function is applied to a gen_fun computed with the compressed parameters
201 * and adapts it to refer to the original parameters.
203 * That is, if y are the compressed parameters and x = A y + b are the original
204 * parameters, then we want the coefficient of the monomial t^y in the original
205 * generating function to be the coefficient of the monomial u^x in the resulting
206 * generating function.
207 * The original generating function has the form
209 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
211 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
213 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
215 * = a u^{A m + b}/(1-u^{A n})
217 * Therefore, we multiply the powers m and n in both numerator and denominator by A
218 * and add b to the power in the numerator.
219 * Since the above powers are stored as row vectors m^T and n^T,
220 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
222 * The pair (map, offset) contains the same information as CP.
223 * map is the transpose of the linear part of CP, while offset is the constant part.
225 void gen_fun::substitute(Matrix
*CP
)
229 split_param_compression(CP
, map
, offset
);
230 Polyhedron
*C
= Polyhedron_Image(context
, CP
, 0);
231 Polyhedron_Free(context
);
233 for (int i
= 0; i
< term
.size(); ++i
) {
234 term
[i
]->d
.power
*= map
;
235 term
[i
]->n
.power
*= map
;
236 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
237 term
[i
]->n
.power
[j
] += offset
;
243 vector
<pair
<Vector
*, QQ
> > vertices
;
244 cone(int *pos
) : pos(pos
) {}
247 #ifndef HAVE_COMPRESS_PARMS
248 static Matrix
*compress_parms(Matrix
*M
, unsigned nparam
)
254 struct parallel_polytopes
{
262 parallel_polytopes(int n
, Polyhedron
*context
, int dim
, int nparam
) :
263 dim(dim
), nparam(nparam
) {
264 red
= gf_base::create(Polyhedron_Copy(context
), dim
, nparam
);
269 void add(const QQ
& c
, Polyhedron
*P
, unsigned MaxRays
) {
270 Polyhedron
*Q
= remove_equalities_p(Polyhedron_Copy(P
), P
->Dimension
-nparam
,
272 POL_ENSURE_VERTICES(Q
);
282 M
= Matrix_Alloc(Q
->NbEq
, Q
->Dimension
+2);
283 Vector_Copy(Q
->Constraint
[0], M
->p
[0], Q
->NbEq
* (Q
->Dimension
+2));
284 CP
= compress_parms(M
, nparam
);
285 T
= align_matrix(CP
, Q
->Dimension
+1);
288 R
= Polyhedron_Preimage(Q
, T
, MaxRays
);
290 Q
= remove_equalities_p(R
, R
->Dimension
-nparam
, NULL
);
292 assert(Q
->NbEq
== 0);
293 assert(Q
->Dimension
== dim
);
295 if (First_Non_Zero(Q
->Constraint
[Q
->NbConstraints
-1]+1, Q
->Dimension
) == -1)
300 Constraints
= Matrix_Alloc(Q
->NbConstraints
, Q
->Dimension
);
301 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
302 Vector_Copy(Q
->Constraint
[i
]+1, Constraints
->p
[i
], Q
->Dimension
);
305 for (int i
= 0; i
< Q
->NbConstraints
; ++i
) {
307 for (j
= 0; j
< Constraints
->NbRows
; ++j
)
308 if (Vector_Equal(Q
->Constraint
[i
]+1, Constraints
->p
[j
],
311 assert(j
< Constraints
->NbRows
);
315 for (int i
= 0; i
< Q
->NbRays
; ++i
) {
316 if (!value_pos_p(Q
->Ray
[i
][dim
+1]))
319 Polyhedron
*C
= supporting_cone(Q
, i
);
321 if (First_Non_Zero(C
->Constraint
[C
->NbConstraints
-1]+1,
325 int *pos
= new int[1+C
->NbConstraints
];
326 pos
[0] = C
->NbConstraints
;
328 for (int k
= 0; k
< Constraints
->NbRows
; ++k
) {
329 for (int j
= 0; j
< C
->NbConstraints
; ++j
) {
330 if (Vector_Equal(C
->Constraint
[j
]+1, Constraints
->p
[k
],
337 assert(l
== C
->NbConstraints
);
340 for (j
= 0; j
< cones
.size(); ++j
)
341 if (!memcmp(pos
, cones
[j
].pos
, (1+C
->NbConstraints
)*sizeof(int)))
343 if (j
== cones
.size())
344 cones
.push_back(cone(pos
));
351 for (k
= 0; k
< cones
[j
].vertices
.size(); ++k
)
352 if (Vector_Equal(Q
->Ray
[i
]+1, cones
[j
].vertices
[k
].first
->p
,
356 if (k
== cones
[j
].vertices
.size()) {
357 Vector
*vertex
= Vector_Alloc(Q
->Dimension
+1);
358 Vector_Copy(Q
->Ray
[i
]+1, vertex
->p
, Q
->Dimension
+1);
359 cones
[j
].vertices
.push_back(pair
<Vector
*,QQ
>(vertex
, c
));
361 cones
[j
].vertices
[k
].second
+= c
;
362 if (cones
[j
].vertices
[k
].second
.n
== 0) {
363 int size
= cones
[j
].vertices
.size();
364 Vector_Free(cones
[j
].vertices
[k
].first
);
366 cones
[j
].vertices
[k
] = cones
[j
].vertices
[size
-1];
367 cones
[j
].vertices
.pop_back();
374 gen_fun
*compute(unsigned MaxRays
) {
375 for (int i
= 0; i
< cones
.size(); ++i
) {
376 Matrix
*M
= Matrix_Alloc(cones
[i
].pos
[0], 1+Constraints
->NbColumns
+1);
378 for (int j
= 0; j
<cones
[i
].pos
[0]; ++j
) {
379 value_set_si(M
->p
[j
][0], 1);
380 Vector_Copy(Constraints
->p
[cones
[i
].pos
[1+j
]], M
->p
[j
]+1,
381 Constraints
->NbColumns
);
383 Cone
= Constraints2Polyhedron(M
, MaxRays
);
385 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
386 red
->base
->do_vertex_cone(cones
[i
].vertices
[j
].second
,
387 Polyhedron_Copy(Cone
),
388 cones
[i
].vertices
[j
].first
->p
,
391 Polyhedron_Free(Cone
);
394 red
->gf
->substitute(CP
);
397 void print(std::ostream
& os
) const {
398 for (int i
= 0; i
< cones
.size(); ++i
) {
400 for (int j
= 0; j
< cones
[i
].pos
[0]; ++j
) {
403 os
<< cones
[i
].pos
[1+j
];
406 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
) {
407 Vector_Print(stderr
, P_VALUE_FMT
, cones
[i
].vertices
[j
].first
);
408 os
<< cones
[i
].vertices
[j
].second
<< endl
;
412 ~parallel_polytopes() {
413 for (int i
= 0; i
< cones
.size(); ++i
) {
414 delete [] cones
[i
].pos
;
415 for (int j
= 0; j
< cones
[i
].vertices
.size(); ++j
)
416 Vector_Free(cones
[i
].vertices
[j
].first
);
419 Matrix_Free(Constraints
);
428 gen_fun
*gen_fun::Hadamard_product(const gen_fun
*gf
, unsigned MaxRays
)
431 Polyhedron
*C
= DomainIntersection(context
, gf
->context
, MaxRays
);
432 Polyhedron
*U
= Universe_Polyhedron(C
->Dimension
);
433 gen_fun
*sum
= new gen_fun(C
);
434 for (int i
= 0; i
< term
.size(); ++i
) {
435 for (int i2
= 0; i2
< gf
->term
.size(); ++i2
) {
436 int d
= term
[i
]->d
.power
.NumCols();
437 int k1
= term
[i
]->d
.power
.NumRows();
438 int k2
= gf
->term
[i2
]->d
.power
.NumRows();
439 assert(term
[i
]->d
.power
.NumCols() == gf
->term
[i2
]->d
.power
.NumCols());
441 parallel_polytopes
pp(term
[i
]->n
.power
.NumRows() *
442 gf
->term
[i2
]->n
.power
.NumRows(),
443 sum
->context
, k1
+k2
-d
, d
);
445 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
) {
446 for (int j2
= 0; j2
< gf
->term
[i2
]->n
.power
.NumRows(); ++j2
) {
447 Matrix
*M
= Matrix_Alloc(k1
+k2
+d
+d
, 1+k1
+k2
+d
+1);
448 for (int k
= 0; k
< k1
+k2
; ++k
) {
449 value_set_si(M
->p
[k
][0], 1);
450 value_set_si(M
->p
[k
][1+k
], 1);
452 for (int k
= 0; k
< d
; ++k
) {
453 value_set_si(M
->p
[k1
+k2
+k
][1+k1
+k2
+k
], -1);
454 zz2value(term
[i
]->n
.power
[j
][k
], M
->p
[k1
+k2
+k
][1+k1
+k2
+d
]);
455 for (int l
= 0; l
< k1
; ++l
)
456 zz2value(term
[i
]->d
.power
[l
][k
], M
->p
[k1
+k2
+k
][1+l
]);
458 for (int k
= 0; k
< d
; ++k
) {
459 value_set_si(M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+k
], -1);
460 zz2value(gf
->term
[i2
]->n
.power
[j2
][k
],
461 M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+d
]);
462 for (int l
= 0; l
< k2
; ++l
)
463 zz2value(gf
->term
[i2
]->d
.power
[l
][k
],
464 M
->p
[k1
+k2
+d
+k
][1+k1
+l
]);
466 Polyhedron
*P
= Constraints2Polyhedron(M
, MaxRays
);
469 QQ c
= term
[i
]->n
.coeff
[j
];
470 c
*= gf
->term
[i2
]->n
.coeff
[j2
];
471 pp
.add(c
, P
, MaxRays
);
477 gen_fun
*t
= pp
.compute(MaxRays
);
486 void gen_fun::add_union(gen_fun
*gf
, unsigned MaxRays
)
488 QQ
one(1, 1), mone(-1, 1);
490 gen_fun
*hp
= Hadamard_product(gf
, MaxRays
);
496 static void Polyhedron_Shift(Polyhedron
*P
, Vector
*offset
)
500 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
501 Inner_Product(P
->Constraint
[i
]+1, offset
->p
, P
->Dimension
, &tmp
);
502 value_subtract(P
->Constraint
[i
][1+P
->Dimension
],
503 P
->Constraint
[i
][1+P
->Dimension
], tmp
);
505 for (int i
= 0; i
< P
->NbRays
; ++i
) {
506 if (value_notone_p(P
->Ray
[i
][0]))
508 if (value_zero_p(P
->Ray
[i
][1+P
->Dimension
]))
510 Vector_Combine(P
->Ray
[i
]+1, offset
->p
, P
->Ray
[i
]+1,
511 P
->Ray
[i
][0], P
->Ray
[i
][1+P
->Dimension
], P
->Dimension
);
516 void gen_fun::shift(const vec_ZZ
& offset
)
518 for (int i
= 0; i
< term
.size(); ++i
)
519 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
520 term
[i
]->n
.power
[j
] += offset
;
522 Vector
*v
= Vector_Alloc(offset
.length());
523 zz2values(offset
, v
->p
);
524 Polyhedron_Shift(context
, v
);
528 /* Divide the generating functin by 1/(1-z^power).
529 * The effect on the corresponding explicit function f(x) is
530 * f'(x) = \sum_{i=0}^\infty f(x - i * power)
532 void gen_fun::divide(const vec_ZZ
& power
)
534 for (int i
= 0; i
< term
.size(); ++i
) {
535 int r
= term
[i
]->d
.power
.NumRows();
536 int c
= term
[i
]->d
.power
.NumCols();
537 term
[i
]->d
.power
.SetDims(r
+1, c
);
538 term
[i
]->d
.power
[r
] = power
;
541 Vector
*v
= Vector_Alloc(1+power
.length()+1);
542 value_set_si(v
->p
[0], 1);
543 zz2values(power
, v
->p
+1);
544 Polyhedron
*C
= AddRays(v
->p
, 1, context
, context
->NbConstraints
+1);
546 Polyhedron_Free(context
);
550 static void print_power(std::ostream
& os
, QQ
& c
, vec_ZZ
& p
,
551 unsigned int nparam
, char **param_name
)
555 for (int i
= 0; i
< p
.length(); ++i
) {
559 if (c
.n
== -1 && c
.d
== 1)
561 else if (c
.n
!= 1 || c
.d
!= 1) {
577 os
<< "^(" << p
[i
] << ")";
588 void gen_fun::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
591 for (int i
= 0; i
< term
.size(); ++i
) {
595 for (int j
= 0; j
< term
[i
]->n
.coeff
.length(); ++j
) {
596 if (j
!= 0 && term
[i
]->n
.coeff
[j
].n
> 0)
598 print_power(os
, term
[i
]->n
.coeff
[j
], term
[i
]->n
.power
[j
],
602 for (int j
= 0; j
< term
[i
]->d
.power
.NumRows(); ++j
) {
606 print_power(os
, mone
, term
[i
]->d
.power
[j
], nparam
, param_name
);
613 gen_fun::operator evalue
*() const
617 value_init(factor
.d
);
618 value_init(factor
.x
.n
);
619 for (int i
= 0; i
< term
.size(); ++i
) {
620 unsigned nvar
= term
[i
]->d
.power
.NumRows();
621 unsigned nparam
= term
[i
]->d
.power
.NumCols();
622 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ nparam
+ 1);
623 mat_ZZ
& d
= term
[i
]->d
.power
;
624 Polyhedron
*U
= context
? context
: Universe_Polyhedron(nparam
);
626 for (int j
= 0; j
< term
[i
]->n
.coeff
.length(); ++j
) {
627 for (int r
= 0; r
< nparam
; ++r
) {
628 value_set_si(C
->p
[r
][0], 0);
629 for (int c
= 0; c
< nvar
; ++c
) {
630 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
632 Vector_Set(&C
->p
[r
][1+nvar
], 0, nparam
);
633 value_set_si(C
->p
[r
][1+nvar
+r
], -1);
634 zz2value(term
[i
]->n
.power
[j
][r
], C
->p
[r
][1+nvar
+nparam
]);
636 for (int r
= 0; r
< nvar
; ++r
) {
637 value_set_si(C
->p
[nparam
+r
][0], 1);
638 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ nparam
+ 1);
639 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
641 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
642 evalue
*E
= barvinok_enumerate_ev(P
, U
, 0);
644 if (EVALUE_IS_ZERO(*E
)) {
649 zz2value(term
[i
]->n
.coeff
[j
].n
, factor
.x
.n
);
650 zz2value(term
[i
]->n
.coeff
[j
].d
, factor
.d
);
653 Matrix_Print(stdout, P_VALUE_FMT, C);
654 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
655 print_evalue(stdout, E, test);
669 value_clear(factor
.d
);
670 value_clear(factor
.x
.n
);
674 void gen_fun::coefficient(Value
* params
, Value
* c
) const
676 if (context
&& !in_domain(context
, params
)) {
683 value_init(part
.x
.n
);
686 evalue_set_si(&sum
, 0, 1);
690 for (int i
= 0; i
< term
.size(); ++i
) {
691 unsigned nvar
= term
[i
]->d
.power
.NumRows();
692 unsigned nparam
= term
[i
]->d
.power
.NumCols();
693 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ 1);
694 mat_ZZ
& d
= term
[i
]->d
.power
;
696 for (int j
= 0; j
< term
[i
]->n
.coeff
.length(); ++j
) {
697 for (int r
= 0; r
< nparam
; ++r
) {
698 value_set_si(C
->p
[r
][0], 0);
699 for (int c
= 0; c
< nvar
; ++c
) {
700 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
702 zz2value(term
[i
]->n
.power
[j
][r
], C
->p
[r
][1+nvar
]);
703 value_subtract(C
->p
[r
][1+nvar
], C
->p
[r
][1+nvar
], params
[r
]);
705 for (int r
= 0; r
< nvar
; ++r
) {
706 value_set_si(C
->p
[nparam
+r
][0], 1);
707 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ 1);
708 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
710 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
715 barvinok_count(P
, &tmp
, 0);
717 if (value_zero_p(tmp
))
719 zz2value(term
[i
]->n
.coeff
[j
].n
, part
.x
.n
);
720 zz2value(term
[i
]->n
.coeff
[j
].d
, part
.d
);
721 value_multiply(part
.x
.n
, part
.x
.n
, tmp
);
727 assert(value_one_p(sum
.d
));
728 value_assign(*c
, sum
.x
.n
);
732 value_clear(part
.x
.n
);
734 value_clear(sum
.x
.n
);
737 gen_fun
*gen_fun::summate(int nvar
) const
739 int dim
= context
->Dimension
;
740 int nparam
= dim
- nvar
;
742 #ifdef USE_INCREMENTAL_DF
743 partial_ireducer
red(Polyhedron_Project(context
, nparam
), dim
, nparam
);
745 partial_reducer
red(Polyhedron_Project(context
, nparam
), dim
, nparam
);
748 for (int i
= 0; i
< term
.size(); ++i
)
749 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
750 red
.reduce(term
[i
]->n
.coeff
[j
], term
[i
]->n
.power
[j
], term
[i
]->d
.power
);
754 /* returns true if the set was finite and false otherwise */
755 bool gen_fun::summate(Value
*sum
) const
757 if (term
.size() == 0) {
758 value_set_si(*sum
, 0);
763 for (int i
= 0; i
< term
.size(); ++i
)
764 if (term
[i
]->d
.power
.NumRows() > maxlen
)
765 maxlen
= term
[i
]->d
.power
.NumRows();
767 infinite_icounter
cnt(term
[0]->d
.power
.NumCols(), maxlen
);
768 for (int i
= 0; i
< term
.size(); ++i
)
769 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
770 cnt
.reduce(term
[i
]->n
.coeff
[j
], term
[i
]->n
.power
[j
], term
[i
]->d
.power
);
772 for (int i
= 1; i
<= maxlen
; ++i
)
773 if (value_notzero_p(mpq_numref(cnt
.count
[i
]))) {
774 value_set_si(*sum
, -1);
778 assert(value_one_p(mpq_denref(cnt
.count
[0])));
779 value_assign(*sum
, mpq_numref(cnt
.count
[0]));