8 #include <NTL/mat_ZZ.h>
10 #include <barvinok/util.h>
12 #include <polylib/polylibgmp.h>
13 #include <barvinok/evalue.h>
17 #include <barvinok/barvinok.h>
18 #include <barvinok/genfun.h>
19 #include "conversion.h"
20 #include "decomposer.h"
21 #include "lattice_point.h"
22 #include "reduce_domain.h"
34 using std::ostringstream
;
36 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
38 static void rays(mat_ZZ
& r
, Polyhedron
*C
)
40 unsigned dim
= C
->NbRays
- 1; /* don't count zero vertex */
41 assert(C
->NbRays
- 1 == C
->Dimension
);
46 for (i
= 0, c
= 0; i
< dim
; ++i
)
47 if (value_zero_p(C
->Ray
[i
][dim
+1])) {
48 for (int j
= 0; j
< dim
; ++j
) {
49 value2zz(C
->Ray
[i
][j
+1], tmp
);
59 dpoly(int d
, ZZ
& degree
, int offset
= 0) {
63 if (degree
>= 0 && degree
< ZZ(INIT_VAL
, min
))
66 ZZ c
= ZZ(INIT_VAL
, 1);
69 for (int i
= 1; i
<= min
; ++i
) {
75 void operator *= (dpoly
& f
) {
76 assert(coeff
.length() == f
.coeff
.length());
78 coeff
= f
.coeff
[0] * coeff
;
79 for (int i
= 1; i
< coeff
.length(); ++i
)
80 for (int j
= 0; i
+j
< coeff
.length(); ++j
)
81 coeff
[i
+j
] += f
.coeff
[i
] * old
[j
];
83 void div(dpoly
& d
, mpq_t count
, ZZ
& sign
) {
84 int len
= coeff
.length();
87 mpq_t
* c
= new mpq_t
[coeff
.length()];
90 for (int i
= 0; i
< len
; ++i
) {
92 zz2value(coeff
[i
], tmp
);
95 for (int j
= 1; j
<= i
; ++j
) {
96 zz2value(d
.coeff
[j
], tmp
);
98 mpq_mul(qtmp
, qtmp
, c
[i
-j
]);
99 mpq_sub(c
[i
], c
[i
], qtmp
);
102 zz2value(d
.coeff
[0], tmp
);
103 mpq_set_z(qtmp
, tmp
);
104 mpq_div(c
[i
], c
[i
], qtmp
);
107 mpq_sub(count
, count
, c
[len
-1]);
109 mpq_add(count
, count
, c
[len
-1]);
113 for (int i
= 0; i
< len
; ++i
)
125 dpoly_n(int d
, ZZ
& degree_0
, ZZ
& degree_1
, int offset
= 0) {
129 zz2value(degree_0
, d0
);
130 zz2value(degree_1
, d1
);
131 coeff
= Matrix_Alloc(d
+1, d
+1+1);
132 value_set_si(coeff
->p
[0][0], 1);
133 value_set_si(coeff
->p
[0][d
+1], 1);
134 for (int i
= 1; i
<= d
; ++i
) {
135 value_multiply(coeff
->p
[i
][0], coeff
->p
[i
-1][0], d0
);
136 Vector_Combine(coeff
->p
[i
-1], coeff
->p
[i
-1]+1, coeff
->p
[i
]+1,
138 value_set_si(coeff
->p
[i
][d
+1], i
);
139 value_multiply(coeff
->p
[i
][d
+1], coeff
->p
[i
][d
+1], coeff
->p
[i
-1][d
+1]);
140 value_decrement(d0
, d0
);
145 void div(dpoly
& d
, Vector
*count
, ZZ
& sign
) {
146 int len
= coeff
->NbRows
;
147 Matrix
* c
= Matrix_Alloc(coeff
->NbRows
, coeff
->NbColumns
);
150 for (int i
= 0; i
< len
; ++i
) {
151 Vector_Copy(coeff
->p
[i
], c
->p
[i
], len
+1);
152 for (int j
= 1; j
<= i
; ++j
) {
153 zz2value(d
.coeff
[j
], tmp
);
154 value_multiply(tmp
, tmp
, c
->p
[i
][len
]);
155 value_oppose(tmp
, tmp
);
156 Vector_Combine(c
->p
[i
], c
->p
[i
-j
], c
->p
[i
],
157 c
->p
[i
-j
][len
], tmp
, len
);
158 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], c
->p
[i
-j
][len
]);
160 zz2value(d
.coeff
[0], tmp
);
161 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], tmp
);
164 value_set_si(tmp
, -1);
165 Vector_Scale(c
->p
[len
-1], count
->p
, tmp
, len
);
166 value_assign(count
->p
[len
], c
->p
[len
-1][len
]);
168 Vector_Copy(c
->p
[len
-1], count
->p
, len
+1);
169 Vector_Normalize(count
->p
, len
+1);
175 struct dpoly_r_term
{
180 /* len: number of elements in c
181 * each element in c is the coefficient of a power of t
182 * in the MacLaurin expansion
185 vector
< dpoly_r_term
* > *c
;
190 void add_term(int i
, int * powers
, ZZ
& coeff
) {
193 for (int k
= 0; k
< c
[i
].size(); ++k
) {
194 if (memcmp(c
[i
][k
]->powers
, powers
, dim
* sizeof(int)) == 0) {
195 c
[i
][k
]->coeff
+= coeff
;
199 dpoly_r_term
*t
= new dpoly_r_term
;
200 t
->powers
= new int[dim
];
201 memcpy(t
->powers
, powers
, dim
* sizeof(int));
205 dpoly_r(int len
, int dim
) {
209 c
= new vector
< dpoly_r_term
* > [len
];
211 dpoly_r(dpoly
& num
, int dim
) {
213 len
= num
.coeff
.length();
214 c
= new vector
< dpoly_r_term
* > [len
];
217 memset(powers
, 0, dim
* sizeof(int));
219 for (int i
= 0; i
< len
; ++i
) {
220 ZZ coeff
= num
.coeff
[i
];
221 add_term(i
, powers
, coeff
);
224 dpoly_r(dpoly
& num
, dpoly
& den
, int pos
, int dim
) {
226 len
= num
.coeff
.length();
227 c
= new vector
< dpoly_r_term
* > [len
];
231 for (int i
= 0; i
< len
; ++i
) {
232 ZZ coeff
= num
.coeff
[i
];
233 memset(powers
, 0, dim
* sizeof(int));
236 add_term(i
, powers
, coeff
);
238 for (int j
= 1; j
<= i
; ++j
) {
239 for (int k
= 0; k
< c
[i
-j
].size(); ++k
) {
240 memcpy(powers
, c
[i
-j
][k
]->powers
, dim
*sizeof(int));
242 coeff
= -den
.coeff
[j
-1] * c
[i
-j
][k
]->coeff
;
243 add_term(i
, powers
, coeff
);
249 dpoly_r(dpoly_r
* num
, dpoly
& den
, int pos
, int dim
) {
252 c
= new vector
< dpoly_r_term
* > [len
];
257 for (int i
= 0 ; i
< len
; ++i
) {
258 for (int k
= 0; k
< num
->c
[i
].size(); ++k
) {
259 memcpy(powers
, num
->c
[i
][k
]->powers
, dim
*sizeof(int));
261 add_term(i
, powers
, num
->c
[i
][k
]->coeff
);
264 for (int j
= 1; j
<= i
; ++j
) {
265 for (int k
= 0; k
< c
[i
-j
].size(); ++k
) {
266 memcpy(powers
, c
[i
-j
][k
]->powers
, dim
*sizeof(int));
268 coeff
= -den
.coeff
[j
-1] * c
[i
-j
][k
]->coeff
;
269 add_term(i
, powers
, coeff
);
275 for (int i
= 0 ; i
< len
; ++i
)
276 for (int k
= 0; k
< c
[i
].size(); ++k
) {
277 delete [] c
[i
][k
]->powers
;
282 dpoly_r
*div(dpoly
& d
) {
283 dpoly_r
*rc
= new dpoly_r(len
, dim
);
284 rc
->denom
= power(d
.coeff
[0], len
);
285 ZZ inv_d
= rc
->denom
/ d
.coeff
[0];
288 for (int i
= 0; i
< len
; ++i
) {
289 for (int k
= 0; k
< c
[i
].size(); ++k
) {
290 coeff
= c
[i
][k
]->coeff
* inv_d
;
291 rc
->add_term(i
, c
[i
][k
]->powers
, coeff
);
294 for (int j
= 1; j
<= i
; ++j
) {
295 for (int k
= 0; k
< rc
->c
[i
-j
].size(); ++k
) {
296 coeff
= - d
.coeff
[j
] * rc
->c
[i
-j
][k
]->coeff
/ d
.coeff
[0];
297 rc
->add_term(i
, rc
->c
[i
-j
][k
]->powers
, coeff
);
304 for (int i
= 0; i
< len
; ++i
) {
307 cerr
<< c
[i
].size() << endl
;
308 for (int j
= 0; j
< c
[i
].size(); ++j
) {
309 for (int k
= 0; k
< dim
; ++k
) {
310 cerr
<< c
[i
][j
]->powers
[k
] << " ";
312 cerr
<< ": " << c
[i
][j
]->coeff
<< "/" << denom
<< endl
;
319 const int MAX_TRY
=10;
321 * Searches for a vector that is not orthogonal to any
322 * of the rays in rays.
324 static void nonorthog(mat_ZZ
& rays
, vec_ZZ
& lambda
)
326 int dim
= rays
.NumCols();
328 lambda
.SetLength(dim
);
332 for (int i
= 2; !found
&& i
<= 50*dim
; i
+=4) {
333 for (int j
= 0; j
< MAX_TRY
; ++j
) {
334 for (int k
= 0; k
< dim
; ++k
) {
335 int r
= random_int(i
)+2;
336 int v
= (2*(r
%2)-1) * (r
>> 1);
340 for (; k
< rays
.NumRows(); ++k
)
341 if (lambda
* rays
[k
] == 0)
343 if (k
== rays
.NumRows()) {
352 static void randomvector(Polyhedron
*P
, vec_ZZ
& lambda
, int nvar
)
356 unsigned int dim
= P
->Dimension
;
359 for (int i
= 0; i
< P
->NbRays
; ++i
) {
360 for (int j
= 1; j
<= dim
; ++j
) {
361 value_absolute(tmp
, P
->Ray
[i
][j
]);
362 int t
= VALUE_TO_LONG(tmp
) * 16;
367 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
368 for (int j
= 1; j
<= dim
; ++j
) {
369 value_absolute(tmp
, P
->Constraint
[i
][j
]);
370 int t
= VALUE_TO_LONG(tmp
) * 16;
377 lambda
.SetLength(nvar
);
378 for (int k
= 0; k
< nvar
; ++k
) {
379 int r
= random_int(max
*dim
)+2;
380 int v
= (2*(r
%2)-1) * (max
/2*dim
+ (r
>> 1));
385 static void add_rays(mat_ZZ
& rays
, Polyhedron
*i
, int *r
, int nvar
= -1,
388 unsigned dim
= i
->Dimension
;
391 for (int k
= 0; k
< i
->NbRays
; ++k
) {
392 if (!value_zero_p(i
->Ray
[k
][dim
+1]))
394 if (!all
&& nvar
!= dim
&& First_Non_Zero(i
->Ray
[k
]+1, nvar
) == -1)
396 values2zz(i
->Ray
[k
]+1, rays
[(*r
)++], nvar
);
400 static void mask_r(Matrix
*f
, int nr
, Vector
*lcm
, int p
, Vector
*val
, evalue
*ev
)
402 unsigned nparam
= lcm
->Size
;
405 Vector
* prod
= Vector_Alloc(f
->NbRows
);
406 Matrix_Vector_Product(f
, val
->p
, prod
->p
);
408 for (int i
= 0; i
< nr
; ++i
) {
409 value_modulus(prod
->p
[i
], prod
->p
[i
], f
->p
[i
][nparam
+1]);
410 isint
&= value_zero_p(prod
->p
[i
]);
412 value_set_si(ev
->d
, 1);
414 value_set_si(ev
->x
.n
, isint
);
421 if (value_one_p(lcm
->p
[p
]))
422 mask_r(f
, nr
, lcm
, p
+1, val
, ev
);
424 value_assign(tmp
, lcm
->p
[p
]);
425 value_set_si(ev
->d
, 0);
426 ev
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
428 value_decrement(tmp
, tmp
);
429 value_assign(val
->p
[p
], tmp
);
430 mask_r(f
, nr
, lcm
, p
+1, val
, &ev
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
431 } while (value_pos_p(tmp
));
437 static void mask(Matrix
*f
, evalue
*factor
)
439 int nr
= f
->NbRows
, nc
= f
->NbColumns
;
442 for (n
= 0; n
< nr
&& value_notzero_p(f
->p
[n
][nc
-1]); ++n
)
443 if (value_notone_p(f
->p
[n
][nc
-1]) &&
444 value_notmone_p(f
->p
[n
][nc
-1]))
458 value_set_si(EV
.x
.n
, 1);
460 for (n
= 0; n
< nr
; ++n
) {
461 value_assign(m
, f
->p
[n
][nc
-1]);
462 if (value_one_p(m
) || value_mone_p(m
))
465 int j
= normal_mod(f
->p
[n
], nc
-1, &m
);
467 free_evalue_refs(factor
);
468 value_init(factor
->d
);
469 evalue_set_si(factor
, 0, 1);
473 values2zz(f
->p
[n
], row
, nc
-1);
476 if (j
< (nc
-1)-1 && row
[j
] > g
/2) {
477 for (int k
= j
; k
< (nc
-1); ++k
)
483 value_set_si(EP
.d
, 0);
484 EP
.x
.p
= new_enode(relation
, 2, 0);
485 value_clear(EP
.x
.p
->arr
[1].d
);
486 EP
.x
.p
->arr
[1] = *factor
;
487 evalue
*ev
= &EP
.x
.p
->arr
[0];
488 value_set_si(ev
->d
, 0);
489 ev
->x
.p
= new_enode(fractional
, 3, -1);
490 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
491 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
492 evalue
*E
= multi_monom(row
);
493 value_assign(EV
.d
, m
);
495 value_clear(ev
->x
.p
->arr
[0].d
);
496 ev
->x
.p
->arr
[0] = *E
;
502 free_evalue_refs(&EV
);
508 static void mask(Matrix
*f
, evalue
*factor
)
510 int nr
= f
->NbRows
, nc
= f
->NbColumns
;
513 for (n
= 0; n
< nr
&& value_notzero_p(f
->p
[n
][nc
-1]); ++n
)
514 if (value_notone_p(f
->p
[n
][nc
-1]) &&
515 value_notmone_p(f
->p
[n
][nc
-1]))
523 unsigned np
= nc
- 2;
524 Vector
*lcm
= Vector_Alloc(np
);
525 Vector
*val
= Vector_Alloc(nc
);
526 Vector_Set(val
->p
, 0, nc
);
527 value_set_si(val
->p
[np
], 1);
528 Vector_Set(lcm
->p
, 1, np
);
529 for (n
= 0; n
< nr
; ++n
) {
530 if (value_one_p(f
->p
[n
][nc
-1]) ||
531 value_mone_p(f
->p
[n
][nc
-1]))
533 for (int j
= 0; j
< np
; ++j
)
534 if (value_notzero_p(f
->p
[n
][j
])) {
535 Gcd(f
->p
[n
][j
], f
->p
[n
][nc
-1], &tmp
);
536 value_division(tmp
, f
->p
[n
][nc
-1], tmp
);
537 value_lcm(tmp
, lcm
->p
[j
], &lcm
->p
[j
]);
542 mask_r(f
, nr
, lcm
, 0, val
, &EP
);
547 free_evalue_refs(&EP
);
551 /* This structure encodes the power of the term in a rational generating function.
553 * Either E == NULL or constant = 0
554 * If E != NULL, then the power is E
555 * If E == NULL, then the power is coeff * param[pos] + constant
564 /* Returns the power of (t+1) in the term of a rational generating function,
565 * i.e., the scalar product of the actual lattice point and lambda.
566 * The lattice point is the unique lattice point in the fundamental parallelepiped
567 * of the unimodual cone i shifted to the parametric vertex V.
569 * PD is the parameter domain, which, if != NULL, may be used to simply the
570 * resulting expression.
572 * The result is returned in term.
575 Param_Vertices
* V
, Polyhedron
*i
, vec_ZZ
& lambda
, term_info
* term
,
578 unsigned nparam
= V
->Vertex
->NbColumns
- 2;
579 unsigned dim
= i
->Dimension
;
581 vertex
.SetDims(V
->Vertex
->NbRows
, nparam
+1);
585 value_set_si(lcm
, 1);
586 for (int j
= 0; j
< V
->Vertex
->NbRows
; ++j
) {
587 value_lcm(lcm
, V
->Vertex
->p
[j
][nparam
+1], &lcm
);
589 if (value_notone_p(lcm
)) {
590 Matrix
* mv
= Matrix_Alloc(dim
, nparam
+1);
591 for (int j
= 0 ; j
< dim
; ++j
) {
592 value_division(tmp
, lcm
, V
->Vertex
->p
[j
][nparam
+1]);
593 Vector_Scale(V
->Vertex
->p
[j
], mv
->p
[j
], tmp
, nparam
+1);
596 term
->E
= lattice_point(i
, lambda
, mv
, lcm
, PD
);
604 for (int i
= 0; i
< V
->Vertex
->NbRows
; ++i
) {
605 assert(value_one_p(V
->Vertex
->p
[i
][nparam
+1])); // for now
606 values2zz(V
->Vertex
->p
[i
], vertex
[i
], nparam
+1);
610 num
= lambda
* vertex
;
614 for (int j
= 0; j
< nparam
; ++j
)
620 term
->E
= multi_monom(num
);
624 term
->constant
= num
[nparam
];
627 term
->coeff
= num
[p
];
634 static void normalize(ZZ
& sign
, ZZ
& num
, vec_ZZ
& den
)
636 unsigned dim
= den
.length();
640 for (int j
= 0; j
< den
.length(); ++j
) {
644 den
[j
] = abs(den
[j
]);
653 * f: the powers in the denominator for the remaining vars
654 * each row refers to a factor
655 * den_s: for each factor, the power of (s+1)
657 * num_s: powers in the numerator corresponding to the summed vars
658 * num_p: powers in the numerator corresponding to the remaining vars
659 * number of rays in cone: "dim" = "k"
660 * length of each ray: "dim" = "d"
661 * for now, it is assumed: k == d
663 * den_p: for each factor
664 * 0: independent of remaining vars
665 * 1: power corresponds to corresponding row in f
667 * all inputs are subject to change
669 static void normalize(ZZ
& sign
,
670 ZZ
& num_s
, vec_ZZ
& num_p
, vec_ZZ
& den_s
, vec_ZZ
& den_p
,
673 unsigned dim
= f
.NumRows();
674 unsigned nparam
= num_p
.length();
675 unsigned nvar
= dim
- nparam
;
679 for (int j
= 0; j
< den_s
.length(); ++j
) {
685 for (k
= 0; k
< nparam
; ++k
)
699 den_s
[j
] = abs(den_s
[j
]);
708 struct counter
: public polar_decomposer
{
720 counter(Polyhedron
*P
) {
723 rays
.SetDims(dim
, dim
);
728 void start(unsigned MaxRays
);
734 virtual void handle_polar(Polyhedron
*P
, int sign
);
737 struct OrthogonalException
{} Orthogonal
;
739 void counter::handle_polar(Polyhedron
*C
, int s
)
742 assert(C
->NbRays
-1 == dim
);
743 add_rays(rays
, C
, &r
);
744 for (int k
= 0; k
< dim
; ++k
) {
745 if (lambda
* rays
[k
] == 0)
751 lattice_point(P
->Ray
[j
]+1, C
, vertex
);
752 num
= vertex
* lambda
;
754 normalize(sign
, num
, den
);
757 dpoly
n(dim
, den
[0], 1);
758 for (int k
= 1; k
< dim
; ++k
) {
759 dpoly
fact(dim
, den
[k
], 1);
762 d
.div(n
, count
, sign
);
765 void counter::start(unsigned MaxRays
)
769 randomvector(P
, lambda
, dim
);
770 for (j
= 0; j
< P
->NbRays
; ++j
) {
771 Polyhedron
*C
= supporting_cone(P
, j
);
772 decompose(C
, MaxRays
);
775 } catch (OrthogonalException
&e
) {
776 mpq_set_si(count
, 0, 0);
781 /* base for non-parametric counting */
782 struct np_base
: public polar_decomposer
{
787 np_base(Polyhedron
*P
, unsigned dim
) {
793 struct reducer
: public np_base
{
802 int lower
; // call base when only this many variables is left
804 reducer(Polyhedron
*P
) : np_base(P
, P
->Dimension
) {
805 //den.SetLength(dim);
812 void start(unsigned MaxRays
);
820 virtual void handle_polar(Polyhedron
*P
, int sign
);
821 void reduce(ZZ c
, ZZ cd
, vec_ZZ
& num
, mat_ZZ
& den_f
);
822 virtual void base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
) = 0;
823 virtual void split(vec_ZZ
& num
, ZZ
& num_s
, vec_ZZ
& num_p
,
824 mat_ZZ
& den_f
, vec_ZZ
& den_s
, mat_ZZ
& den_r
) = 0;
827 void reducer::reduce(ZZ c
, ZZ cd
, vec_ZZ
& num
, mat_ZZ
& den_f
)
829 unsigned len
= den_f
.NumRows(); // number of factors in den
831 if (num
.length() == lower
) {
832 base(c
, cd
, num
, den_f
);
835 assert(num
.length() > 1);
842 split(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
845 den_p
.SetLength(len
);
847 normalize(c
, num_s
, num_p
, den_s
, den_p
, den_r
);
849 int only_param
= 0; // k-r-s from text
850 int no_param
= 0; // r from text
851 for (int k
= 0; k
< len
; ++k
) {
854 else if (den_s
[k
] == 0)
858 reduce(c
, cd
, num_p
, den_r
);
862 pden
.SetDims(only_param
, den_r
.NumCols());
864 for (k
= 0, l
= 0; k
< len
; ++k
)
866 pden
[l
++] = den_r
[k
];
868 for (k
= 0; k
< len
; ++k
)
872 dpoly
n(no_param
, num_s
);
873 dpoly
D(no_param
, den_s
[k
], 1);
876 dpoly
fact(no_param
, den_s
[k
], 1);
880 if (no_param
+ only_param
== len
) {
881 mpq_set_si(tcount
, 0, 1);
882 n
.div(D
, tcount
, one
);
885 value2zz(mpq_numref(tcount
), qn
);
886 value2zz(mpq_denref(tcount
), qd
);
892 reduce(qn
, qd
, num_p
, pden
);
896 for (k
= 0; k
< len
; ++k
) {
897 if (den_s
[k
] == 0 || den_p
[k
] == 0)
900 dpoly
pd(no_param
-1, den_s
[k
], 1);
903 for (l
= 0; l
< k
; ++l
)
904 if (den_r
[l
] == den_r
[k
])
908 r
= new dpoly_r(n
, pd
, l
, len
);
910 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
916 dpoly_r
*rc
= r
->div(D
);
920 int common
= pden
.NumRows();
921 vector
< dpoly_r_term
* >& final
= rc
->c
[rc
->len
-1];
923 for (int j
= 0; j
< final
.size(); ++j
) {
924 if (final
[j
]->coeff
== 0)
927 pden
.SetDims(rows
, pden
.NumCols());
928 for (int k
= 0; k
< rc
->dim
; ++k
) {
929 int n
= final
[j
]->powers
[k
];
932 pden
.SetDims(rows
+n
, pden
.NumCols());
933 for (int l
= 0; l
< n
; ++l
)
934 pden
[rows
+l
] = den_r
[k
];
937 final
[j
]->coeff
*= c
;
938 reduce(final
[j
]->coeff
, rc
->denom
, num_p
, pden
);
947 void reducer::handle_polar(Polyhedron
*C
, int s
)
949 assert(C
->NbRays
-1 == dim
);
953 lattice_point(P
->Ray
[current_vertex
]+1, C
, vertex
);
956 den
.SetDims(dim
, dim
);
959 for (r
= 0; r
< dim
; ++r
)
960 values2zz(C
->Ray
[r
]+1, den
[r
], dim
);
962 reduce(sgn
, one
, vertex
, den
);
965 void reducer::start(unsigned MaxRays
)
967 for (current_vertex
= 0; current_vertex
< P
->NbRays
; ++current_vertex
) {
968 Polyhedron
*C
= supporting_cone(P
, current_vertex
);
969 decompose(C
, MaxRays
);
973 struct ireducer
: public reducer
{
974 ireducer(Polyhedron
*P
) : reducer(P
) {}
976 virtual void split(vec_ZZ
& num
, ZZ
& num_s
, vec_ZZ
& num_p
,
977 mat_ZZ
& den_f
, vec_ZZ
& den_s
, mat_ZZ
& den_r
) {
978 unsigned len
= den_f
.NumRows(); // number of factors in den
979 unsigned d
= num
.length() - 1;
981 den_s
.SetLength(len
);
982 den_r
.SetDims(len
, d
);
984 for (int r
= 0; r
< len
; ++r
) {
985 den_s
[r
] = den_f
[r
][0];
986 for (int k
= 1; k
<= d
; ++k
)
987 den_r
[r
][k
-1] = den_f
[r
][k
];
992 for (int k
= 1 ; k
<= d
; ++k
)
997 // incremental counter
998 struct icounter
: public ireducer
{
1001 icounter(Polyhedron
*P
) : ireducer(P
) {
1008 virtual void base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
);
1011 void icounter::base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
)
1014 unsigned len
= den_f
.NumRows(); // number of factors in den
1016 den_s
.SetLength(len
);
1018 for (r
= 0; r
< len
; ++r
)
1019 den_s
[r
] = den_f
[r
][0];
1020 normalize(c
, num_s
, den_s
);
1022 dpoly
n(len
, num_s
);
1023 dpoly
D(len
, den_s
[0], 1);
1024 for (int k
= 1; k
< len
; ++k
) {
1025 dpoly
fact(len
, den_s
[k
], 1);
1028 mpq_set_si(tcount
, 0, 1);
1029 n
.div(D
, tcount
, one
);
1032 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
1033 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
1034 mpq_canonicalize(tcount
);
1035 mpq_add(count
, count
, tcount
);
1038 /* base for generating function counting */
1043 gf_base(np_base
*npb
, unsigned nparam
) : base(npb
) {
1044 gf
= new gen_fun(Polyhedron_Project(base
->P
, nparam
));
1046 void start(unsigned MaxRays
);
1049 void gf_base::start(unsigned MaxRays
)
1051 for (int i
= 0; i
< base
->P
->NbRays
; ++i
) {
1052 if (!value_pos_p(base
->P
->Ray
[i
][base
->dim
+1]))
1055 Polyhedron
*C
= supporting_cone(base
->P
, i
);
1056 base
->current_vertex
= i
;
1057 base
->decompose(C
, MaxRays
);
1061 struct partial_ireducer
: public ireducer
, public gf_base
{
1062 partial_ireducer(Polyhedron
*P
, unsigned nparam
) :
1063 ireducer(P
), gf_base(this, nparam
) {
1066 ~partial_ireducer() {
1068 virtual void base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
);
1069 /* we want to override the start method from reducer with the one from gf_base */
1070 void start(unsigned MaxRays
) {
1071 gf_base::start(MaxRays
);
1075 void partial_ireducer::base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
)
1077 gf
->add(c
, cd
, num
, den_f
);
1080 struct partial_reducer
: public reducer
, public gf_base
{
1084 partial_reducer(Polyhedron
*P
, unsigned nparam
) :
1085 reducer(P
), gf_base(this, nparam
) {
1088 tmp
.SetLength(dim
- nparam
);
1089 randomvector(P
, lambda
, dim
- nparam
);
1091 ~partial_reducer() {
1093 virtual void base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
);
1094 /* we want to override the start method from reducer with the one from gf_base */
1095 void start(unsigned MaxRays
) {
1096 gf_base::start(MaxRays
);
1099 virtual void split(vec_ZZ
& num
, ZZ
& num_s
, vec_ZZ
& num_p
,
1100 mat_ZZ
& den_f
, vec_ZZ
& den_s
, mat_ZZ
& den_r
) {
1101 unsigned len
= den_f
.NumRows(); // number of factors in den
1102 unsigned nvar
= tmp
.length();
1104 den_s
.SetLength(len
);
1105 den_r
.SetDims(len
, lower
);
1107 for (int r
= 0; r
< len
; ++r
) {
1108 for (int k
= 0; k
< nvar
; ++k
)
1109 tmp
[k
] = den_f
[r
][k
];
1110 den_s
[r
] = tmp
* lambda
;
1112 for (int k
= nvar
; k
< dim
; ++k
)
1113 den_r
[r
][k
-nvar
] = den_f
[r
][k
];
1116 for (int k
= 0; k
< nvar
; ++k
)
1118 num_s
= tmp
*lambda
;
1119 num_p
.SetLength(lower
);
1120 for (int k
= nvar
; k
< dim
; ++k
)
1121 num_p
[k
-nvar
] = num
[k
];
1125 void partial_reducer::base(ZZ
& c
, ZZ
& cd
, vec_ZZ
& num
, mat_ZZ
& den_f
)
1127 gf
->add(c
, cd
, num
, den_f
);
1130 struct bfc_term_base
{
1131 // the number of times a given factor appears in the denominator
1135 bfc_term_base(int len
) {
1136 powers
= new int[len
];
1139 virtual ~bfc_term_base() {
1144 struct bfc_term
: public bfc_term_base
{
1148 bfc_term(int len
) : bfc_term_base(len
) {}
1151 struct bfe_term
: public bfc_term_base
{
1152 vector
<evalue
*> factors
;
1154 bfe_term(int len
) : bfc_term_base(len
) {
1158 for (int i
= 0; i
< factors
.size(); ++i
) {
1161 free_evalue_refs(factors
[i
]);
1167 typedef vector
< bfc_term_base
* > bfc_vec
;
1171 struct bf_base
: public np_base
{
1176 int lower
; // call base when only this many variables is left
1178 bf_base(Polyhedron
*P
, unsigned dim
) : np_base(P
, dim
) {
1191 void start(unsigned MaxRays
);
1192 virtual void handle_polar(Polyhedron
*P
, int sign
);
1193 int setup_factors(Polyhedron
*P
, mat_ZZ
& factors
, bfc_term_base
* t
, int s
);
1195 bfc_term_base
* find_bfc_term(bfc_vec
& v
, int *powers
, int len
);
1196 void add_term(bfc_term_base
*t
, vec_ZZ
& num1
, vec_ZZ
& num
);
1197 void add_term(bfc_term_base
*t
, vec_ZZ
& num
);
1199 void reduce(mat_ZZ
& factors
, bfc_vec
& v
);
1200 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
) = 0;
1202 virtual bfc_term_base
* new_bf_term(int len
) = 0;
1203 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) = 0;
1204 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&f
, int change
) = 0;
1205 virtual void set_factor(bfc_term_base
*t
, int k
, ZZ
& n
, ZZ
& d
, int change
) = 0;
1206 virtual void update_term(bfc_term_base
*t
, int i
) = 0;
1207 virtual void insert_term(bfc_term_base
*t
, int i
) = 0;
1208 virtual bool constant_vertex(int dim
) = 0;
1209 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
,
1213 static int lex_cmp(vec_ZZ
& a
, vec_ZZ
& b
)
1215 assert(a
.length() == b
.length());
1217 for (int j
= 0; j
< a
.length(); ++j
)
1219 return a
[j
] < b
[j
] ? -1 : 1;
1223 void bf_base::add_term(bfc_term_base
*t
, vec_ZZ
& num_orig
, vec_ZZ
& extra_num
)
1226 int d
= num_orig
.length();
1228 for (int l
= 0; l
< d
-1; ++l
)
1229 num
[l
] = num_orig
[l
+1] + extra_num
[l
];
1234 void bf_base::add_term(bfc_term_base
*t
, vec_ZZ
& num
)
1236 int len
= t
->terms
.NumRows();
1238 for (i
= 0; i
< len
; ++i
) {
1239 r
= lex_cmp(t
->terms
[i
], num
);
1243 if (i
== len
|| r
> 0) {
1244 t
->terms
.SetDims(len
+1, num
.length());
1248 // i < len && r == 0
1253 static void print_int_vector(int *v
, int len
, char *name
)
1255 cerr
<< name
<< endl
;
1256 for (int j
= 0; j
< len
; ++j
) {
1257 cerr
<< v
[j
] << " ";
1262 static void print_bfc_terms(mat_ZZ
& factors
, bfc_vec
& v
)
1265 cerr
<< "factors" << endl
;
1266 cerr
<< factors
<< endl
;
1267 for (int i
= 0; i
< v
.size(); ++i
) {
1268 cerr
<< "term: " << i
<< endl
;
1269 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
1270 cerr
<< "terms" << endl
;
1271 cerr
<< v
[i
]->terms
<< endl
;
1272 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
1273 cerr
<< bfct
->cn
<< endl
;
1274 cerr
<< bfct
->cd
<< endl
;
1278 static void print_bfe_terms(mat_ZZ
& factors
, bfc_vec
& v
)
1281 cerr
<< "factors" << endl
;
1282 cerr
<< factors
<< endl
;
1283 for (int i
= 0; i
< v
.size(); ++i
) {
1284 cerr
<< "term: " << i
<< endl
;
1285 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
1286 cerr
<< "terms" << endl
;
1287 cerr
<< v
[i
]->terms
<< endl
;
1288 bfe_term
* bfet
= static_cast<bfe_term
*>(v
[i
]);
1289 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
) {
1290 char * test
[] = {"a", "b"};
1291 print_evalue(stderr
, bfet
->factors
[j
], test
);
1292 fprintf(stderr
, "\n");
1297 bfc_term_base
* bf_base::find_bfc_term(bfc_vec
& v
, int *powers
, int len
)
1299 bfc_vec::iterator i
;
1300 for (i
= v
.begin(); i
!= v
.end(); ++i
) {
1302 for (j
= 0; j
< len
; ++j
)
1303 if ((*i
)->powers
[j
] != powers
[j
])
1307 if ((*i
)->powers
[j
] > powers
[j
])
1311 bfc_term_base
* t
= new_bf_term(len
);
1313 memcpy(t
->powers
, powers
, len
* sizeof(int));
1334 int no_param
; // r from text
1335 int only_param
; // k-r-s from text
1336 int total_power
; // k from text
1338 // created in compute_reduced_factors
1340 // set in update_powers
1345 bf_reducer(mat_ZZ
& factors
, bfc_vec
& v
, bf_base
*bf
)
1346 : factors(factors
), v(v
), bf(bf
) {
1347 nf
= factors
.NumRows();
1348 d
= factors
.NumCols();
1349 old2new
= new int[nf
];
1352 extra_num
.SetLength(d
-1);
1361 void compute_reduced_factors();
1362 void compute_extra_num(int i
);
1366 void update_powers(int *powers
, int len
);
1369 void bf_reducer::compute_extra_num(int i
)
1373 no_param
= 0; // r from text
1374 only_param
= 0; // k-r-s from text
1375 total_power
= 0; // k from text
1377 for (int j
= 0; j
< nf
; ++j
) {
1378 if (v
[i
]->powers
[j
] == 0)
1381 total_power
+= v
[i
]->powers
[j
];
1382 if (factors
[j
][0] == 0) {
1383 only_param
+= v
[i
]->powers
[j
];
1387 if (old2new
[j
] == -1)
1388 no_param
+= v
[i
]->powers
[j
];
1390 extra_num
+= -sign
[j
] * v
[i
]->powers
[j
] * nfactors
[old2new
[j
]];
1391 changes
+= v
[i
]->powers
[j
];
1395 void bf_reducer::update_powers(int *powers
, int len
)
1397 for (int l
= 0; l
< nnf
; ++l
)
1398 npowers
[l
] = bpowers
[l
];
1400 l_extra_num
= extra_num
;
1401 l_changes
= changes
;
1403 for (int l
= 0; l
< len
; ++l
) {
1407 assert(old2new
[l
] != -1);
1409 npowers
[old2new
[l
]] += n
;
1410 // interpretation of sign has been inverted
1411 // since we inverted the power for specialization
1413 l_extra_num
+= n
* nfactors
[old2new
[l
]];
1420 void bf_reducer::compute_reduced_factors()
1422 unsigned nf
= factors
.NumRows();
1423 unsigned d
= factors
.NumCols();
1425 nfactors
.SetDims(nnf
, d
-1);
1427 for (int i
= 0; i
< nf
; ++i
) {
1430 for (j
= 0; j
< nnf
; ++j
) {
1432 for (k
= 1; k
< d
; ++k
)
1433 if (factors
[i
][k
] != 0 || nfactors
[j
][k
-1] != 0)
1435 if (k
< d
&& factors
[i
][k
] == -nfactors
[j
][k
-1])
1438 if (factors
[i
][k
] != s
* nfactors
[j
][k
-1])
1446 for (k
= 1; k
< d
; ++k
)
1447 if (factors
[i
][k
] != 0)
1450 if (factors
[i
][k
] < 0)
1452 nfactors
.SetDims(++nnf
, d
-1);
1453 for (int k
= 1; k
< d
; ++k
)
1454 nfactors
[j
][k
-1] = s
* factors
[i
][k
];
1460 npowers
= new int[nnf
];
1461 bpowers
= new int[nnf
];
1464 void bf_reducer::reduce()
1466 compute_reduced_factors();
1468 for (int i
= 0; i
< v
.size(); ++i
) {
1469 compute_extra_num(i
);
1471 if (no_param
== 0) {
1473 extra_num
.SetLength(d
-1);
1476 for (int k
= 0; k
< nnf
; ++k
)
1478 for (int k
= 0; k
< nf
; ++k
) {
1479 assert(old2new
[k
] != -1);
1480 npowers
[old2new
[k
]] += v
[i
]->powers
[k
];
1481 if (sign
[k
] == -1) {
1482 extra_num
+= v
[i
]->powers
[k
] * nfactors
[old2new
[k
]];
1483 changes
+= v
[i
]->powers
[k
];
1487 bfc_term_base
* t
= bf
->find_bfc_term(vn
, npowers
, nnf
);
1488 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
1489 bf
->set_factor(v
[i
], k
, changes
% 2);
1490 bf
->add_term(t
, v
[i
]->terms
[k
], extra_num
);
1493 // powers of "constant" part
1494 for (int k
= 0; k
< nnf
; ++k
)
1496 for (int k
= 0; k
< nf
; ++k
) {
1497 if (factors
[k
][0] != 0)
1499 assert(old2new
[k
] != -1);
1500 bpowers
[old2new
[k
]] += v
[i
]->powers
[k
];
1501 if (sign
[k
] == -1) {
1502 extra_num
+= v
[i
]->powers
[k
] * nfactors
[old2new
[k
]];
1503 changes
+= v
[i
]->powers
[k
];
1508 for (j
= 0; j
< nf
; ++j
)
1509 if (old2new
[j
] == -1 && v
[i
]->powers
[j
] > 0)
1512 dpoly
D(no_param
, factors
[j
][0], 1);
1513 for (int k
= 1; k
< v
[i
]->powers
[j
]; ++k
) {
1514 dpoly
fact(no_param
, factors
[j
][0], 1);
1518 if (old2new
[j
] == -1)
1519 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
) {
1520 dpoly
fact(no_param
, factors
[j
][0], 1);
1524 if (no_param
+ only_param
== total_power
&&
1525 bf
->constant_vertex(d
)) {
1526 bfc_term_base
* t
= NULL
;
1531 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
1532 dpoly
n(no_param
, v
[i
]->terms
[k
][0]);
1533 mpq_set_si(bf
->tcount
, 0, 1);
1534 n
.div(D
, bf
->tcount
, bf
->one
);
1536 if (value_zero_p(mpq_numref(bf
->tcount
)))
1540 t
= bf
->find_bfc_term(vn
, bpowers
, nnf
);
1541 bf
->set_factor(v
[i
], k
, bf
->tcount
, changes
% 2);
1542 bf
->add_term(t
, v
[i
]->terms
[k
], extra_num
);
1545 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
) {
1546 dpoly
n(no_param
, v
[i
]->terms
[j
][0]);
1549 if (no_param
+ only_param
== total_power
)
1550 r
= new dpoly_r(n
, nf
);
1552 for (int k
= 0; k
< nf
; ++k
) {
1553 if (v
[i
]->powers
[k
] == 0)
1555 if (factors
[k
][0] == 0 || old2new
[k
] == -1)
1558 dpoly
pd(no_param
-1, factors
[k
][0], 1);
1560 for (int l
= 0; l
< v
[i
]->powers
[k
]; ++l
) {
1562 for (q
= 0; q
< k
; ++q
)
1563 if (old2new
[q
] == old2new
[k
] &&
1568 r
= new dpoly_r(n
, pd
, q
, nf
);
1570 dpoly_r
*nr
= new dpoly_r(r
, pd
, q
, nf
);
1577 dpoly_r
*rc
= r
->div(D
);
1580 if (bf
->constant_vertex(d
)) {
1581 vector
< dpoly_r_term
* >& final
= rc
->c
[rc
->len
-1];
1583 for (int k
= 0; k
< final
.size(); ++k
) {
1584 if (final
[k
]->coeff
== 0)
1587 update_powers(final
[k
]->powers
, rc
->dim
);
1589 bfc_term_base
* t
= bf
->find_bfc_term(vn
, npowers
, nnf
);
1590 bf
->set_factor(v
[i
], j
, final
[k
]->coeff
, rc
->denom
, l_changes
% 2);
1591 bf
->add_term(t
, v
[i
]->terms
[j
], l_extra_num
);
1594 bf
->cum(this, v
[i
], j
, rc
);
1605 void bf_base::reduce(mat_ZZ
& factors
, bfc_vec
& v
)
1607 assert(v
.size() > 0);
1608 unsigned nf
= factors
.NumRows();
1609 unsigned d
= factors
.NumCols();
1612 return base(factors
, v
);
1614 bf_reducer
bfr(factors
, v
, this);
1618 if (bfr
.vn
.size() > 0)
1619 reduce(bfr
.nfactors
, bfr
.vn
);
1622 int bf_base::setup_factors(Polyhedron
*C
, mat_ZZ
& factors
,
1623 bfc_term_base
* t
, int s
)
1625 factors
.SetDims(dim
, dim
);
1629 for (r
= 0; r
< dim
; ++r
)
1632 for (r
= 0; r
< dim
; ++r
) {
1633 values2zz(C
->Ray
[r
]+1, factors
[r
], dim
);
1635 for (k
= 0; k
< dim
; ++k
)
1636 if (factors
[r
][k
] != 0)
1638 if (factors
[r
][k
] < 0) {
1639 factors
[r
] = -factors
[r
];
1640 t
->terms
[0] += factors
[r
];
1648 void bf_base::handle_polar(Polyhedron
*C
, int s
)
1650 bfc_term
* t
= new bfc_term(dim
);
1651 vector
< bfc_term_base
* > v
;
1654 assert(C
->NbRays
-1 == dim
);
1659 t
->terms
.SetDims(1, dim
);
1660 lattice_point(P
->Ray
[current_vertex
]+1, C
, t
->terms
[0]);
1662 // the elements of factors are always lexpositive
1664 s
= setup_factors(C
, factors
, t
, s
);
1672 void bf_base::start(unsigned MaxRays
)
1674 for (current_vertex
= 0; current_vertex
< P
->NbRays
; ++current_vertex
) {
1675 Polyhedron
*C
= supporting_cone(P
, current_vertex
);
1676 decompose(C
, MaxRays
);
1680 struct bfcounter_base
: public bf_base
{
1684 bfcounter_base(Polyhedron
*P
) : bf_base(P
, P
->Dimension
) {
1687 bfc_term_base
* new_bf_term(int len
) {
1688 bfc_term
* t
= new bfc_term(len
);
1694 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) {
1695 bfc_term
* bfct
= static_cast<bfc_term
*>(t
);
1702 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&f
, int change
) {
1703 bfc_term
* bfct
= static_cast<bfc_term
*>(t
);
1704 value2zz(mpq_numref(f
), cn
);
1705 value2zz(mpq_denref(f
), cd
);
1712 virtual void set_factor(bfc_term_base
*t
, int k
, ZZ
& n
, ZZ
& d
, int change
) {
1713 bfc_term
* bfct
= static_cast<bfc_term
*>(t
);
1714 cn
= bfct
->cn
[k
] * n
;
1717 cd
= bfct
->cd
[k
] * d
;
1720 virtual void insert_term(bfc_term_base
*t
, int i
) {
1721 bfc_term
* bfct
= static_cast<bfc_term
*>(t
);
1722 int len
= t
->terms
.NumRows()-1; // already increased by one
1724 bfct
->cn
.SetLength(len
+1);
1725 bfct
->cd
.SetLength(len
+1);
1726 for (int j
= len
; j
> i
; --j
) {
1727 bfct
->cn
[j
] = bfct
->cn
[j
-1];
1728 bfct
->cd
[j
] = bfct
->cd
[j
-1];
1729 t
->terms
[j
] = t
->terms
[j
-1];
1735 virtual void update_term(bfc_term_base
*t
, int i
) {
1736 bfc_term
* bfct
= static_cast<bfc_term
*>(t
);
1738 ZZ g
= GCD(bfct
->cd
[i
], cd
);
1739 ZZ n
= cn
* bfct
->cd
[i
]/g
+ bfct
->cn
[i
] * cd
/g
;
1740 ZZ d
= bfct
->cd
[i
] * cd
/ g
;
1745 virtual bool constant_vertex(int dim
) { return true; }
1746 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
, dpoly_r
*r
) {
1751 struct bfcounter
: public bfcounter_base
{
1754 bfcounter(Polyhedron
*P
) : bfcounter_base(P
) {
1761 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
1764 void bfcounter::base(mat_ZZ
& factors
, bfc_vec
& v
)
1766 unsigned nf
= factors
.NumRows();
1768 for (int i
= 0; i
< v
.size(); ++i
) {
1769 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
1770 int total_power
= 0;
1771 // factor is always positive, so we always
1773 for (int k
= 0; k
< nf
; ++k
)
1774 total_power
+= v
[i
]->powers
[k
];
1777 for (j
= 0; j
< nf
; ++j
)
1778 if (v
[i
]->powers
[j
] > 0)
1781 dpoly
D(total_power
, factors
[j
][0], 1);
1782 for (int k
= 1; k
< v
[i
]->powers
[j
]; ++k
) {
1783 dpoly
fact(total_power
, factors
[j
][0], 1);
1787 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
) {
1788 dpoly
fact(total_power
, factors
[j
][0], 1);
1792 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
1793 dpoly
n(total_power
, v
[i
]->terms
[k
][0]);
1794 mpq_set_si(tcount
, 0, 1);
1795 n
.div(D
, tcount
, one
);
1796 if (total_power
% 2)
1797 bfct
->cn
[k
] = -bfct
->cn
[k
];
1798 zz2value(bfct
->cn
[k
], tn
);
1799 zz2value(bfct
->cd
[k
], td
);
1801 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
1802 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
1803 mpq_canonicalize(tcount
);
1804 mpq_add(count
, count
, tcount
);
1810 struct partial_bfcounter
: public bfcounter_base
, public gf_base
{
1811 partial_bfcounter(Polyhedron
*P
, unsigned nparam
) :
1812 bfcounter_base(P
), gf_base(this, nparam
) {
1815 ~partial_bfcounter() {
1817 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
1818 /* we want to override the start method from bf_base with the one from gf_base */
1819 void start(unsigned MaxRays
) {
1820 gf_base::start(MaxRays
);
1824 void partial_bfcounter::base(mat_ZZ
& factors
, bfc_vec
& v
)
1827 unsigned nf
= factors
.NumRows();
1829 for (int i
= 0; i
< v
.size(); ++i
) {
1830 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
1831 den
.SetDims(0, lower
);
1832 int total_power
= 0;
1834 for (int j
= 0; j
< nf
; ++j
) {
1835 total_power
+= v
[i
]->powers
[j
];
1836 den
.SetDims(total_power
, lower
);
1837 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
)
1838 den
[p
++] = factors
[j
];
1840 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
)
1841 gf
->add(bfct
->cn
[j
], bfct
->cd
[j
], v
[i
]->terms
[j
], den
);
1847 /* Check whether the polyhedron is unbounded and if so,
1848 * check whether it has any (and therefore an infinite number of)
1850 * If one of the vertices is integer, then we are done.
1851 * Otherwise, transform the polyhedron such that one of the rays
1852 * is the first unit vector and cut it off at a height that ensures
1853 * that if the whole polyhedron has any points, then the remaining part
1854 * has integer points. In particular we add the largest coefficient
1855 * of a ray to the highest vertex (rounded up).
1857 static bool Polyhedron_is_infinite(Polyhedron
*P
, Value
* result
, unsigned MaxRays
)
1869 for (; r
< P
->NbRays
; ++r
)
1870 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
1872 if (P
->NbBid
== 0 && r
== P
->NbRays
)
1878 sample
= Polyhedron_Sample(P
, MaxRays
);
1880 value_set_si(*result
, 0);
1882 value_set_si(*result
, -1);
1883 Vector_Free(sample
);
1888 for (int i
= 0; i
< P
->NbRays
; ++i
)
1889 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
1890 value_set_si(*result
, -1);
1895 v
= Vector_Alloc(P
->Dimension
+1);
1896 Vector_Gcd(P
->Ray
[r
]+1, P
->Dimension
, &g
);
1897 Vector_AntiScale(P
->Ray
[r
]+1, v
->p
, g
, P
->Dimension
+1);
1898 M
= unimodular_complete(v
);
1899 value_set_si(M
->p
[P
->Dimension
][P
->Dimension
], 1);
1902 P
= Polyhedron_Preimage(P
, M2
, 0);
1911 value_set_si(size
, 0);
1913 for (int i
= 0; i
< P
->NbBid
; ++i
) {
1914 value_absolute(tmp
, P
->Ray
[i
][1]);
1915 if (value_gt(tmp
, size
))
1916 value_assign(size
, tmp
);
1918 for (int i
= P
->NbBid
; i
< P
->NbRays
; ++i
) {
1919 if (value_zero_p(P
->Ray
[i
][P
->Dimension
+1])) {
1920 if (value_gt(P
->Ray
[i
][1], size
))
1921 value_assign(size
, P
->Ray
[i
][1]);
1924 mpz_cdiv_q(tmp
, P
->Ray
[i
][1], P
->Ray
[i
][P
->Dimension
+1]);
1925 if (first
|| value_gt(tmp
, offset
)) {
1926 value_assign(offset
, tmp
);
1930 value_addto(offset
, offset
, size
);
1934 v
= Vector_Alloc(P
->Dimension
+2);
1935 value_set_si(v
->p
[0], 1);
1936 value_set_si(v
->p
[1], -1);
1937 value_assign(v
->p
[1+P
->Dimension
], offset
);
1938 R
= AddConstraints(v
->p
, 1, P
, MaxRays
);
1942 value_clear(offset
);
1946 barvinok_count(P
, &c
, MaxRays
);
1948 if (value_zero_p(c
))
1949 value_set_si(*result
, 0);
1951 value_set_si(*result
, -1);
1957 typedef Polyhedron
* Polyhedron_p
;
1959 static void barvinok_count_f(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
);
1961 void barvinok_count(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
)
1966 bool infinite
= false;
1969 value_set_si(*result
, 0);
1973 P
= remove_equalities(P
);
1976 value_set_si(*result
, 0);
1981 if (Polyhedron_is_infinite(P
, result
, NbMaxCons
)) {
1986 if (P
->Dimension
== 0) {
1987 /* Test whether the constraints are satisfied */
1988 POL_ENSURE_VERTICES(P
);
1989 value_set_si(*result
, !emptyQ(P
));
1994 Q
= Polyhedron_Factor(P
, 0, NbMaxCons
);
2002 barvinok_count_f(P
, result
, NbMaxCons
);
2003 if (value_neg_p(*result
))
2005 if (Q
&& P
->next
&& value_notzero_p(*result
)) {
2009 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
2010 barvinok_count_f(Q
, &factor
, NbMaxCons
);
2011 if (value_neg_p(factor
)) {
2014 } else if (Q
->next
&& value_zero_p(factor
)) {
2015 value_set_si(*result
, 0);
2018 value_multiply(*result
, *result
, factor
);
2021 value_clear(factor
);
2027 value_set_si(*result
, -1);
2030 static void barvinok_count_f(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
)
2033 value_set_si(*result
, 0);
2037 if (P
->Dimension
== 1)
2038 return Line_Length(P
, result
);
2040 int c
= P
->NbConstraints
;
2041 POL_ENSURE_FACETS(P
);
2042 if (c
!= P
->NbConstraints
|| P
->NbEq
!= 0)
2043 return barvinok_count(P
, result
, NbMaxCons
);
2045 POL_ENSURE_VERTICES(P
);
2047 #ifdef USE_INCREMENTAL_BF
2049 #elif defined USE_INCREMENTAL_DF
2054 cnt
.start(NbMaxCons
);
2056 assert(value_one_p(&cnt
.count
[0]._mp_den
));
2057 value_assign(*result
, &cnt
.count
[0]._mp_num
);
2060 static void uni_polynom(int param
, Vector
*c
, evalue
*EP
)
2062 unsigned dim
= c
->Size
-2;
2064 value_set_si(EP
->d
,0);
2065 EP
->x
.p
= new_enode(polynomial
, dim
+1, param
+1);
2066 for (int j
= 0; j
<= dim
; ++j
)
2067 evalue_set(&EP
->x
.p
->arr
[j
], c
->p
[j
], c
->p
[dim
+1]);
2070 static void multi_polynom(Vector
*c
, evalue
* X
, evalue
*EP
)
2072 unsigned dim
= c
->Size
-2;
2076 evalue_set(&EC
, c
->p
[dim
], c
->p
[dim
+1]);
2079 evalue_set(EP
, c
->p
[dim
], c
->p
[dim
+1]);
2081 for (int i
= dim
-1; i
>= 0; --i
) {
2083 value_assign(EC
.x
.n
, c
->p
[i
]);
2086 free_evalue_refs(&EC
);
2089 Polyhedron
*unfringe (Polyhedron
*P
, unsigned MaxRays
)
2091 int len
= P
->Dimension
+2;
2092 Polyhedron
*T
, *R
= P
;
2095 Vector
*row
= Vector_Alloc(len
);
2096 value_set_si(row
->p
[0], 1);
2098 R
= DomainConstraintSimplify(Polyhedron_Copy(P
), MaxRays
);
2100 Matrix
*M
= Matrix_Alloc(2, len
-1);
2101 value_set_si(M
->p
[1][len
-2], 1);
2102 for (int v
= 0; v
< P
->Dimension
; ++v
) {
2103 value_set_si(M
->p
[0][v
], 1);
2104 Polyhedron
*I
= Polyhedron_Image(P
, M
, 2+1);
2105 value_set_si(M
->p
[0][v
], 0);
2106 for (int r
= 0; r
< I
->NbConstraints
; ++r
) {
2107 if (value_zero_p(I
->Constraint
[r
][0]))
2109 if (value_zero_p(I
->Constraint
[r
][1]))
2111 if (value_one_p(I
->Constraint
[r
][1]))
2113 if (value_mone_p(I
->Constraint
[r
][1]))
2115 value_absolute(g
, I
->Constraint
[r
][1]);
2116 Vector_Set(row
->p
+1, 0, len
-2);
2117 value_division(row
->p
[1+v
], I
->Constraint
[r
][1], g
);
2118 mpz_fdiv_q(row
->p
[len
-1], I
->Constraint
[r
][2], g
);
2120 R
= AddConstraints(row
->p
, 1, R
, MaxRays
);
2132 /* this procedure may have false negatives */
2133 static bool Polyhedron_is_infinite_param(Polyhedron
*P
, unsigned nparam
)
2136 for (r
= 0; r
< P
->NbRays
; ++r
) {
2137 if (!value_zero_p(P
->Ray
[r
][0]) &&
2138 !value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
2140 if (First_Non_Zero(P
->Ray
[r
]+1+P
->Dimension
-nparam
, nparam
) == -1)
2146 /* Check whether all rays point in the positive directions
2147 * for the parameters
2149 static bool Polyhedron_has_positive_rays(Polyhedron
*P
, unsigned nparam
)
2152 for (r
= 0; r
< P
->NbRays
; ++r
)
2153 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1])) {
2155 for (i
= P
->Dimension
- nparam
; i
< P
->Dimension
; ++i
)
2156 if (value_neg_p(P
->Ray
[r
][i
+1]))
2162 typedef evalue
* evalue_p
;
2164 struct enumerator
: public polar_decomposer
{
2178 enumerator(Polyhedron
*P
, unsigned dim
, unsigned nbV
) {
2182 randomvector(P
, lambda
, dim
);
2183 rays
.SetDims(dim
, dim
);
2185 c
= Vector_Alloc(dim
+2);
2187 vE
= new evalue_p
[nbV
];
2188 for (int j
= 0; j
< nbV
; ++j
)
2194 void decompose_at(Param_Vertices
*V
, int _i
, unsigned MaxRays
) {
2195 Polyhedron
*C
= supporting_cone_p(P
, V
);
2199 vE
[_i
] = new evalue
;
2200 value_init(vE
[_i
]->d
);
2201 evalue_set_si(vE
[_i
], 0, 1);
2203 decompose(C
, MaxRays
);
2210 for (int j
= 0; j
< nbV
; ++j
)
2212 free_evalue_refs(vE
[j
]);
2218 virtual void handle_polar(Polyhedron
*P
, int sign
);
2221 void enumerator::handle_polar(Polyhedron
*C
, int s
)
2224 assert(C
->NbRays
-1 == dim
);
2225 add_rays(rays
, C
, &r
);
2226 for (int k
= 0; k
< dim
; ++k
) {
2227 if (lambda
* rays
[k
] == 0)
2233 lattice_point(V
, C
, lambda
, &num
, 0);
2234 den
= rays
* lambda
;
2235 normalize(sign
, num
.constant
, den
);
2237 dpoly
n(dim
, den
[0], 1);
2238 for (int k
= 1; k
< dim
; ++k
) {
2239 dpoly
fact(dim
, den
[k
], 1);
2242 if (num
.E
!= NULL
) {
2243 ZZ
one(INIT_VAL
, 1);
2244 dpoly_n
d(dim
, num
.constant
, one
);
2247 multi_polynom(c
, num
.E
, &EV
);
2249 free_evalue_refs(&EV
);
2250 free_evalue_refs(num
.E
);
2252 } else if (num
.pos
!= -1) {
2253 dpoly_n
d(dim
, num
.constant
, num
.coeff
);
2256 uni_polynom(num
.pos
, c
, &EV
);
2258 free_evalue_refs(&EV
);
2260 mpq_set_si(count
, 0, 1);
2261 dpoly
d(dim
, num
.constant
);
2262 d
.div(n
, count
, sign
);
2265 evalue_set(&EV
, &count
[0]._mp_num
, &count
[0]._mp_den
);
2267 free_evalue_refs(&EV
);
2271 struct enumerator_base
{
2276 vertex_decomposer
*vpd
;
2278 enumerator_base(unsigned dim
, vertex_decomposer
*vpd
)
2283 vE
= new evalue_p
[vpd
->nbV
];
2284 for (int j
= 0; j
< vpd
->nbV
; ++j
)
2287 E_vertex
= new evalue_p
[dim
];
2290 evalue_set_si(&mone
, -1, 1);
2293 void decompose_at(Param_Vertices
*V
, int _i
, unsigned MaxRays
/*, Polyhedron *pVD*/) {
2296 vE
[_i
] = new evalue
;
2297 value_init(vE
[_i
]->d
);
2298 evalue_set_si(vE
[_i
], 0, 1);
2300 vpd
->decompose_at_vertex(V
, _i
, MaxRays
);
2303 ~enumerator_base() {
2304 for (int j
= 0; j
< vpd
->nbV
; ++j
)
2306 free_evalue_refs(vE
[j
]);
2313 free_evalue_refs(&mone
);
2316 evalue
*E_num(int i
, int d
) {
2317 return E_vertex
[i
+ (dim
-d
)];
2326 cumulator(evalue
*factor
, evalue
*v
, dpoly_r
*r
) :
2327 factor(factor
), v(v
), r(r
) {}
2331 virtual void add_term(int *powers
, int len
, evalue
*f2
) = 0;
2334 void cumulator::cumulate()
2336 evalue cum
; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
2338 evalue t
; // E_num[0] - (m-1)
2344 evalue_set_si(&mone
, -1, 1);
2348 evalue_copy(&cum
, factor
);
2351 value_set_si(f
.d
, 1);
2352 value_set_si(f
.x
.n
, 1);
2357 for (cst
= &t
; value_zero_p(cst
->d
); ) {
2358 if (cst
->x
.p
->type
== fractional
)
2359 cst
= &cst
->x
.p
->arr
[1];
2361 cst
= &cst
->x
.p
->arr
[0];
2365 for (int m
= 0; m
< r
->len
; ++m
) {
2368 value_set_si(f
.d
, m
);
2371 value_subtract(cst
->x
.n
, cst
->x
.n
, cst
->d
);
2378 vector
< dpoly_r_term
* >& current
= r
->c
[r
->len
-1-m
];
2379 for (int j
= 0; j
< current
.size(); ++j
) {
2380 if (current
[j
]->coeff
== 0)
2382 evalue
*f2
= new evalue
;
2384 value_init(f2
->x
.n
);
2385 zz2value(current
[j
]->coeff
, f2
->x
.n
);
2386 zz2value(r
->denom
, f2
->d
);
2389 add_term(current
[j
]->powers
, r
->dim
, f2
);
2392 free_evalue_refs(&f
);
2393 free_evalue_refs(&t
);
2394 free_evalue_refs(&cum
);
2396 free_evalue_refs(&mone
);
2400 struct E_poly_term
{
2405 struct ie_cum
: public cumulator
{
2406 vector
<E_poly_term
*> terms
;
2408 ie_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
) : cumulator(factor
, v
, r
) {}
2410 virtual void add_term(int *powers
, int len
, evalue
*f2
);
2413 void ie_cum::add_term(int *powers
, int len
, evalue
*f2
)
2416 for (k
= 0; k
< terms
.size(); ++k
) {
2417 if (memcmp(terms
[k
]->powers
, powers
, len
* sizeof(int)) == 0) {
2418 eadd(f2
, terms
[k
]->E
);
2419 free_evalue_refs(f2
);
2424 if (k
>= terms
.size()) {
2425 E_poly_term
*ET
= new E_poly_term
;
2426 ET
->powers
= new int[len
];
2427 memcpy(ET
->powers
, powers
, len
* sizeof(int));
2429 terms
.push_back(ET
);
2433 struct ienumerator
: public polar_decomposer
, public vertex_decomposer
,
2434 public enumerator_base
{
2440 ienumerator(Polyhedron
*P
, unsigned dim
, unsigned nbV
) :
2441 vertex_decomposer(P
, nbV
, this), enumerator_base(dim
, this) {
2442 vertex
.SetLength(dim
);
2444 den
.SetDims(dim
, dim
);
2452 virtual void handle_polar(Polyhedron
*P
, int sign
);
2453 void reduce(evalue
*factor
, vec_ZZ
& num
, mat_ZZ
& den_f
);
2456 void ienumerator::reduce(
2457 evalue
*factor
, vec_ZZ
& num
, mat_ZZ
& den_f
)
2459 unsigned len
= den_f
.NumRows(); // number of factors in den
2460 unsigned dim
= num
.length();
2463 eadd(factor
, vE
[vert
]);
2468 den_s
.SetLength(len
);
2470 den_r
.SetDims(len
, dim
-1);
2474 for (r
= 0; r
< len
; ++r
) {
2475 den_s
[r
] = den_f
[r
][0];
2476 for (k
= 0; k
<= dim
-1; ++k
)
2478 den_r
[r
][k
-(k
>0)] = den_f
[r
][k
];
2483 num_p
.SetLength(dim
-1);
2484 for (k
= 0 ; k
<= dim
-1; ++k
)
2486 num_p
[k
-(k
>0)] = num
[k
];
2489 den_p
.SetLength(len
);
2493 normalize(one
, num_s
, num_p
, den_s
, den_p
, den_r
);
2495 emul(&mone
, factor
);
2499 for (int k
= 0; k
< len
; ++k
) {
2502 else if (den_s
[k
] == 0)
2505 if (no_param
== 0) {
2506 reduce(factor
, num_p
, den_r
);
2510 pden
.SetDims(only_param
, dim
-1);
2512 for (k
= 0, l
= 0; k
< len
; ++k
)
2514 pden
[l
++] = den_r
[k
];
2516 for (k
= 0; k
< len
; ++k
)
2520 dpoly
n(no_param
, num_s
);
2521 dpoly
D(no_param
, den_s
[k
], 1);
2522 for ( ; ++k
< len
; )
2523 if (den_p
[k
] == 0) {
2524 dpoly
fact(no_param
, den_s
[k
], 1);
2529 // if no_param + only_param == len then all powers
2530 // below will be all zero
2531 if (no_param
+ only_param
== len
) {
2532 if (E_num(0, dim
) != 0)
2533 r
= new dpoly_r(n
, len
);
2535 mpq_set_si(tcount
, 0, 1);
2537 n
.div(D
, tcount
, one
);
2539 if (value_notzero_p(mpq_numref(tcount
))) {
2543 value_assign(f
.x
.n
, mpq_numref(tcount
));
2544 value_assign(f
.d
, mpq_denref(tcount
));
2546 reduce(factor
, num_p
, pden
);
2547 free_evalue_refs(&f
);
2552 for (k
= 0; k
< len
; ++k
) {
2553 if (den_s
[k
] == 0 || den_p
[k
] == 0)
2556 dpoly
pd(no_param
-1, den_s
[k
], 1);
2559 for (l
= 0; l
< k
; ++l
)
2560 if (den_r
[l
] == den_r
[k
])
2564 r
= new dpoly_r(n
, pd
, l
, len
);
2566 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
2572 dpoly_r
*rc
= r
->div(D
);
2575 if (E_num(0, dim
) == 0) {
2576 int common
= pden
.NumRows();
2577 vector
< dpoly_r_term
* >& final
= r
->c
[r
->len
-1];
2583 zz2value(r
->denom
, f
.d
);
2584 for (int j
= 0; j
< final
.size(); ++j
) {
2585 if (final
[j
]->coeff
== 0)
2588 for (int k
= 0; k
< r
->dim
; ++k
) {
2589 int n
= final
[j
]->powers
[k
];
2592 pden
.SetDims(rows
+n
, pden
.NumCols());
2593 for (int l
= 0; l
< n
; ++l
)
2594 pden
[rows
+l
] = den_r
[k
];
2598 evalue_copy(&t
, factor
);
2599 zz2value(final
[j
]->coeff
, f
.x
.n
);
2601 reduce(&t
, num_p
, pden
);
2602 free_evalue_refs(&t
);
2604 free_evalue_refs(&f
);
2606 ie_cum
cum(factor
, E_num(0, dim
), r
);
2609 int common
= pden
.NumRows();
2611 for (int j
= 0; j
< cum
.terms
.size(); ++j
) {
2613 pden
.SetDims(rows
, pden
.NumCols());
2614 for (int k
= 0; k
< r
->dim
; ++k
) {
2615 int n
= cum
.terms
[j
]->powers
[k
];
2618 pden
.SetDims(rows
+n
, pden
.NumCols());
2619 for (int l
= 0; l
< n
; ++l
)
2620 pden
[rows
+l
] = den_r
[k
];
2623 reduce(cum
.terms
[j
]->E
, num_p
, pden
);
2624 free_evalue_refs(cum
.terms
[j
]->E
);
2625 delete cum
.terms
[j
]->E
;
2626 delete [] cum
.terms
[j
]->powers
;
2627 delete cum
.terms
[j
];
2634 static int type_offset(enode
*p
)
2636 return p
->type
== fractional
? 1 :
2637 p
->type
== flooring
? 1 : 0;
2640 static int edegree(evalue
*e
)
2645 if (value_notzero_p(e
->d
))
2649 int i
= type_offset(p
);
2650 if (p
->size
-i
-1 > d
)
2651 d
= p
->size
- i
- 1;
2652 for (; i
< p
->size
; i
++) {
2653 int d2
= edegree(&p
->arr
[i
]);
2660 void ienumerator::handle_polar(Polyhedron
*C
, int s
)
2662 assert(C
->NbRays
-1 == dim
);
2664 lattice_point(V
, C
, vertex
, E_vertex
);
2667 for (r
= 0; r
< dim
; ++r
)
2668 values2zz(C
->Ray
[r
]+1, den
[r
], dim
);
2672 evalue_set_si(&one
, s
, 1);
2673 reduce(&one
, vertex
, den
);
2674 free_evalue_refs(&one
);
2676 for (int i
= 0; i
< dim
; ++i
)
2678 free_evalue_refs(E_vertex
[i
]);
2683 struct bfenumerator
: public vertex_decomposer
, public bf_base
,
2684 public enumerator_base
{
2687 bfenumerator(Polyhedron
*P
, unsigned dim
, unsigned nbV
) :
2688 vertex_decomposer(P
, nbV
, this),
2689 bf_base(P
, dim
), enumerator_base(dim
, this) {
2697 virtual void handle_polar(Polyhedron
*P
, int sign
);
2698 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
2700 bfc_term_base
* new_bf_term(int len
) {
2701 bfe_term
* t
= new bfe_term(len
);
2705 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) {
2706 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
2707 factor
= bfet
->factors
[k
];
2708 assert(factor
!= NULL
);
2709 bfet
->factors
[k
] = NULL
;
2711 emul(&mone
, factor
);
2714 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&q
, int change
) {
2715 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
2716 factor
= bfet
->factors
[k
];
2717 assert(factor
!= NULL
);
2718 bfet
->factors
[k
] = NULL
;
2724 value_oppose(f
.x
.n
, mpq_numref(q
));
2726 value_assign(f
.x
.n
, mpq_numref(q
));
2727 value_assign(f
.d
, mpq_denref(q
));
2731 virtual void set_factor(bfc_term_base
*t
, int k
, ZZ
& n
, ZZ
& d
, int change
) {
2732 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
2734 factor
= new evalue
;
2741 value_oppose(f
.x
.n
, f
.x
.n
);
2744 value_init(factor
->d
);
2745 evalue_copy(factor
, bfet
->factors
[k
]);
2749 void set_factor(evalue
*f
, int change
) {
2755 virtual void insert_term(bfc_term_base
*t
, int i
) {
2756 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
2757 int len
= t
->terms
.NumRows()-1; // already increased by one
2759 bfet
->factors
.resize(len
+1);
2760 for (int j
= len
; j
> i
; --j
) {
2761 bfet
->factors
[j
] = bfet
->factors
[j
-1];
2762 t
->terms
[j
] = t
->terms
[j
-1];
2764 bfet
->factors
[i
] = factor
;
2768 virtual void update_term(bfc_term_base
*t
, int i
) {
2769 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
2771 eadd(factor
, bfet
->factors
[i
]);
2772 free_evalue_refs(factor
);
2776 virtual bool constant_vertex(int dim
) { return E_num(0, dim
) == 0; }
2778 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
, dpoly_r
*r
);
2781 struct bfe_cum
: public cumulator
{
2783 bfc_term_base
*told
;
2787 bfe_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
, bf_reducer
*bfr
,
2788 bfc_term_base
*t
, int k
, bfenumerator
*e
) :
2789 cumulator(factor
, v
, r
), told(t
), k(k
),
2793 virtual void add_term(int *powers
, int len
, evalue
*f2
);
2796 void bfe_cum::add_term(int *powers
, int len
, evalue
*f2
)
2798 bfr
->update_powers(powers
, len
);
2800 bfc_term_base
* t
= bfe
->find_bfc_term(bfr
->vn
, bfr
->npowers
, bfr
->nnf
);
2801 bfe
->set_factor(f2
, bfr
->l_changes
% 2);
2802 bfe
->add_term(t
, told
->terms
[k
], bfr
->l_extra_num
);
2805 void bfenumerator::cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
,
2808 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
2809 bfe_cum
cum(bfet
->factors
[k
], E_num(0, bfr
->d
), r
, bfr
, t
, k
, this);
2813 void bfenumerator::base(mat_ZZ
& factors
, bfc_vec
& v
)
2815 for (int i
= 0; i
< v
.size(); ++i
) {
2816 assert(v
[i
]->terms
.NumRows() == 1);
2817 evalue
*factor
= static_cast<bfe_term
*>(v
[i
])->factors
[0];
2818 eadd(factor
, vE
[vert
]);
2823 void bfenumerator::handle_polar(Polyhedron
*C
, int s
)
2825 assert(C
->NbRays
-1 == enumerator_base::dim
);
2827 bfe_term
* t
= new bfe_term(enumerator_base::dim
);
2828 vector
< bfc_term_base
* > v
;
2831 t
->factors
.resize(1);
2833 t
->terms
.SetDims(1, enumerator_base::dim
);
2834 lattice_point(V
, C
, t
->terms
[0], E_vertex
);
2836 // the elements of factors are always lexpositive
2838 s
= setup_factors(C
, factors
, t
, s
);
2840 t
->factors
[0] = new evalue
;
2841 value_init(t
->factors
[0]->d
);
2842 evalue_set_si(t
->factors
[0], s
, 1);
2845 for (int i
= 0; i
< enumerator_base::dim
; ++i
)
2847 free_evalue_refs(E_vertex
[i
]);
2852 #ifdef HAVE_CORRECT_VERTICES
2853 static inline Param_Polyhedron
*Polyhedron2Param_SD(Polyhedron
**Din
,
2854 Polyhedron
*Cin
,int WS
,Polyhedron
**CEq
,Matrix
**CT
)
2856 if (WS
& POL_NO_DUAL
)
2858 return Polyhedron2Param_SimplifiedDomain(Din
, Cin
, WS
, CEq
, CT
);
2861 static Param_Polyhedron
*Polyhedron2Param_SD(Polyhedron
**Din
,
2862 Polyhedron
*Cin
,int WS
,Polyhedron
**CEq
,Matrix
**CT
)
2864 static char data
[] = " 1 0 0 0 0 1 -18 "
2865 " 1 0 0 -20 0 19 1 "
2866 " 1 0 1 20 0 -20 16 "
2869 " 1 4 -20 0 0 -1 23 "
2870 " 1 -4 20 0 0 1 -22 "
2871 " 1 0 1 0 20 -20 16 "
2872 " 1 0 0 0 -20 19 1 ";
2873 static int checked
= 0;
2878 Matrix
*M
= Matrix_Alloc(9, 7);
2879 for (i
= 0; i
< 9; ++i
)
2880 for (int j
= 0; j
< 7; ++j
) {
2881 sscanf(p
, "%d%n", &v
, &n
);
2883 value_set_si(M
->p
[i
][j
], v
);
2885 Polyhedron
*P
= Constraints2Polyhedron(M
, 1024);
2887 Polyhedron
*U
= Universe_Polyhedron(1);
2888 Param_Polyhedron
*PP
= Polyhedron2Param_Domain(P
, U
, 1024);
2892 for (i
= 0, V
= PP
->V
; V
; ++i
, V
= V
->next
)
2895 Param_Polyhedron_Free(PP
);
2897 fprintf(stderr
, "WARNING: results may be incorrect\n");
2899 "WARNING: use latest version of PolyLib to remove this warning\n");
2903 return Polyhedron2Param_SimplifiedDomain(Din
, Cin
, WS
, CEq
, CT
);
2907 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
2911 static evalue
* barvinok_enumerate_cst(Polyhedron
*P
, Polyhedron
* C
,
2916 ALLOC(evalue
, eres
);
2917 value_init(eres
->d
);
2918 value_set_si(eres
->d
, 0);
2919 eres
->x
.p
= new_enode(partition
, 2, C
->Dimension
);
2920 EVALUE_SET_DOMAIN(eres
->x
.p
->arr
[0], DomainConstraintSimplify(C
, MaxRays
));
2921 value_set_si(eres
->x
.p
->arr
[1].d
, 1);
2922 value_init(eres
->x
.p
->arr
[1].x
.n
);
2924 value_set_si(eres
->x
.p
->arr
[1].x
.n
, 0);
2926 barvinok_count(P
, &eres
->x
.p
->arr
[1].x
.n
, MaxRays
);
2931 evalue
* barvinok_enumerate_ev(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
2933 //P = unfringe(P, MaxRays);
2934 Polyhedron
*Corig
= C
;
2935 Polyhedron
*CEq
= NULL
, *rVD
, *CA
;
2937 unsigned nparam
= C
->Dimension
;
2941 value_init(factor
.d
);
2942 evalue_set_si(&factor
, 1, 1);
2944 CA
= align_context(C
, P
->Dimension
, MaxRays
);
2945 P
= DomainIntersection(P
, CA
, MaxRays
);
2946 Polyhedron_Free(CA
);
2949 POL_ENSURE_FACETS(P
);
2950 POL_ENSURE_VERTICES(P
);
2951 POL_ENSURE_FACETS(C
);
2952 POL_ENSURE_VERTICES(C
);
2954 if (C
->Dimension
== 0 || emptyQ(P
)) {
2956 eres
= barvinok_enumerate_cst(P
, CEq
? CEq
: Polyhedron_Copy(C
),
2959 emul(&factor
, eres
);
2960 reduce_evalue(eres
);
2961 free_evalue_refs(&factor
);
2968 if (Polyhedron_is_infinite_param(P
, nparam
))
2973 P
= remove_equalities_p(P
, P
->Dimension
-nparam
, &f
);
2977 if (P
->Dimension
== nparam
) {
2979 P
= Universe_Polyhedron(0);
2983 Polyhedron
*T
= Polyhedron_Factor(P
, nparam
, MaxRays
);
2984 if (T
|| (P
->Dimension
== nparam
+1)) {
2987 for (Q
= T
? T
: P
; Q
; Q
= Q
->next
) {
2988 Polyhedron
*next
= Q
->next
;
2992 if (Q
->Dimension
!= C
->Dimension
)
2993 QC
= Polyhedron_Project(Q
, nparam
);
2996 C
= DomainIntersection(C
, QC
, MaxRays
);
2998 Polyhedron_Free(C2
);
3000 Polyhedron_Free(QC
);
3008 if (T
->Dimension
== C
->Dimension
) {
3015 Polyhedron
*next
= P
->next
;
3017 eres
= barvinok_enumerate_ev_f(P
, C
, MaxRays
);
3024 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
3025 Polyhedron
*next
= Q
->next
;
3028 f
= barvinok_enumerate_ev_f(Q
, C
, MaxRays
);
3030 free_evalue_refs(f
);
3040 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
3043 unsigned nparam
= C
->Dimension
;
3045 if (P
->Dimension
- nparam
== 1)
3046 return ParamLine_Length(P
, C
, MaxRays
);
3048 Param_Polyhedron
*PP
= NULL
;
3049 Polyhedron
*CEq
= NULL
, *pVD
;
3051 Param_Domain
*D
, *next
;
3054 Polyhedron
*Porig
= P
;
3056 PP
= Polyhedron2Param_SD(&P
,C
,MaxRays
,&CEq
,&CT
);
3058 if (isIdentity(CT
)) {
3062 assert(CT
->NbRows
!= CT
->NbColumns
);
3063 if (CT
->NbRows
== 1) { // no more parameters
3064 eres
= barvinok_enumerate_cst(P
, CEq
, MaxRays
);
3069 Param_Polyhedron_Free(PP
);
3075 nparam
= CT
->NbRows
- 1;
3078 unsigned dim
= P
->Dimension
- nparam
;
3080 ALLOC(evalue
, eres
);
3081 value_init(eres
->d
);
3082 value_set_si(eres
->d
, 0);
3085 for (nd
= 0, D
=PP
->D
; D
; ++nd
, D
=D
->next
);
3086 struct section
{ Polyhedron
*D
; evalue E
; };
3087 section
*s
= new section
[nd
];
3088 Polyhedron
**fVD
= new Polyhedron_p
[nd
];
3091 #ifdef USE_INCREMENTAL_BF
3092 bfenumerator
et(P
, dim
, PP
->nbV
);
3093 #elif defined USE_INCREMENTAL_DF
3094 ienumerator
et(P
, dim
, PP
->nbV
);
3096 enumerator
et(P
, dim
, PP
->nbV
);
3099 for(nd
= 0, D
=PP
->D
; D
; D
=next
) {
3102 Polyhedron
*rVD
= reduce_domain(D
->Domain
, CT
, CEq
,
3107 pVD
= CT
? DomainImage(rVD
,CT
,MaxRays
) : rVD
;
3109 value_init(s
[nd
].E
.d
);
3110 evalue_set_si(&s
[nd
].E
, 0, 1);
3113 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
3116 et
.decompose_at(V
, _i
, MaxRays
);
3117 } catch (OrthogonalException
&e
) {
3120 for (; nd
>= 0; --nd
) {
3121 free_evalue_refs(&s
[nd
].E
);
3122 Domain_Free(s
[nd
].D
);
3123 Domain_Free(fVD
[nd
]);
3127 eadd(et
.vE
[_i
] , &s
[nd
].E
);
3128 END_FORALL_PVertex_in_ParamPolyhedron
;
3129 evalue_range_reduction_in_domain(&s
[nd
].E
, pVD
);
3132 addeliminatedparams_evalue(&s
[nd
].E
, CT
);
3139 evalue_set_si(eres
, 0, 1);
3141 eres
->x
.p
= new_enode(partition
, 2*nd
, C
->Dimension
);
3142 for (int j
= 0; j
< nd
; ++j
) {
3143 EVALUE_SET_DOMAIN(eres
->x
.p
->arr
[2*j
], s
[j
].D
);
3144 value_clear(eres
->x
.p
->arr
[2*j
+1].d
);
3145 eres
->x
.p
->arr
[2*j
+1] = s
[j
].E
;
3146 Domain_Free(fVD
[j
]);
3153 Polyhedron_Free(CEq
);
3157 Enumeration
* barvinok_enumerate(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
3159 evalue
*EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
3161 return partition2enumeration(EP
);
3164 static void SwapColumns(Value
**V
, int n
, int i
, int j
)
3166 for (int r
= 0; r
< n
; ++r
)
3167 value_swap(V
[r
][i
], V
[r
][j
]);
3170 static void SwapColumns(Polyhedron
*P
, int i
, int j
)
3172 SwapColumns(P
->Constraint
, P
->NbConstraints
, i
, j
);
3173 SwapColumns(P
->Ray
, P
->NbRays
, i
, j
);
3176 /* Construct a constraint c from constraints l and u such that if
3177 * if constraint c holds then for each value of the other variables
3178 * there is at most one value of variable pos (position pos+1 in the constraints).
3180 * Given a lower and an upper bound
3181 * n_l v_i + <c_l,x> + c_l >= 0
3182 * -n_u v_i + <c_u,x> + c_u >= 0
3183 * the constructed constraint is
3185 * -(n_l<c_u,x> + n_u<c_l,x>) + (-n_l c_u - n_u c_l + n_l n_u - 1)
3187 * which is then simplified to remove the content of the non-constant coefficients
3189 * len is the total length of the constraints.
3190 * v is a temporary variable that can be used by this procedure
3192 static void negative_test_constraint(Value
*l
, Value
*u
, Value
*c
, int pos
,
3195 value_oppose(*v
, u
[pos
+1]);
3196 Vector_Combine(l
+1, u
+1, c
+1, *v
, l
[pos
+1], len
-1);
3197 value_multiply(*v
, *v
, l
[pos
+1]);
3198 value_subtract(c
[len
-1], c
[len
-1], *v
);
3199 value_set_si(*v
, -1);
3200 Vector_Scale(c
+1, c
+1, *v
, len
-1);
3201 value_decrement(c
[len
-1], c
[len
-1]);
3202 ConstraintSimplify(c
, c
, len
, v
);
3205 static bool parallel_constraints(Value
*l
, Value
*u
, Value
*c
, int pos
,
3214 Vector_Gcd(&l
[1+pos
], len
, &g1
);
3215 Vector_Gcd(&u
[1+pos
], len
, &g2
);
3216 Vector_Combine(l
+1+pos
, u
+1+pos
, c
+1, g2
, g1
, len
);
3217 parallel
= First_Non_Zero(c
+1, len
) == -1;
3225 static void negative_test_constraint7(Value
*l
, Value
*u
, Value
*c
, int pos
,
3226 int exist
, int len
, Value
*v
)
3231 Vector_Gcd(&u
[1+pos
], exist
, v
);
3232 Vector_Gcd(&l
[1+pos
], exist
, &g
);
3233 Vector_Combine(l
+1, u
+1, c
+1, *v
, g
, len
-1);
3234 value_multiply(*v
, *v
, g
);
3235 value_subtract(c
[len
-1], c
[len
-1], *v
);
3236 value_set_si(*v
, -1);
3237 Vector_Scale(c
+1, c
+1, *v
, len
-1);
3238 value_decrement(c
[len
-1], c
[len
-1]);
3239 ConstraintSimplify(c
, c
, len
, v
);
3244 /* Turns a x + b >= 0 into a x + b <= -1
3246 * len is the total length of the constraint.
3247 * v is a temporary variable that can be used by this procedure
3249 static void oppose_constraint(Value
*c
, int len
, Value
*v
)
3251 value_set_si(*v
, -1);
3252 Vector_Scale(c
+1, c
+1, *v
, len
-1);
3253 value_decrement(c
[len
-1], c
[len
-1]);
3256 /* Split polyhedron P into two polyhedra *pos and *neg, where
3257 * existential variable i has at most one solution for each
3258 * value of the other variables in *neg.
3260 * The splitting is performed using constraints l and u.
3262 * nvar: number of set variables
3263 * row: temporary vector that can be used by this procedure
3264 * f: temporary value that can be used by this procedure
3266 static bool SplitOnConstraint(Polyhedron
*P
, int i
, int l
, int u
,
3267 int nvar
, int MaxRays
, Vector
*row
, Value
& f
,
3268 Polyhedron
**pos
, Polyhedron
**neg
)
3270 negative_test_constraint(P
->Constraint
[l
], P
->Constraint
[u
],
3271 row
->p
, nvar
+i
, P
->Dimension
+2, &f
);
3272 *neg
= AddConstraints(row
->p
, 1, P
, MaxRays
);
3274 /* We found an independent, but useless constraint
3275 * Maybe we should detect this earlier and not
3276 * mark the variable as INDEPENDENT
3278 if (emptyQ((*neg
))) {
3279 Polyhedron_Free(*neg
);
3283 oppose_constraint(row
->p
, P
->Dimension
+2, &f
);
3284 *pos
= AddConstraints(row
->p
, 1, P
, MaxRays
);
3286 if (emptyQ((*pos
))) {
3287 Polyhedron_Free(*neg
);
3288 Polyhedron_Free(*pos
);
3296 * unimodularly transform P such that constraint r is transformed
3297 * into a constraint that involves only a single (the first)
3298 * existential variable
3301 static Polyhedron
*rotate_along(Polyhedron
*P
, int r
, int nvar
, int exist
,
3307 Vector
*row
= Vector_Alloc(exist
);
3308 Vector_Copy(P
->Constraint
[r
]+1+nvar
, row
->p
, exist
);
3309 Vector_Gcd(row
->p
, exist
, &g
);
3310 if (value_notone_p(g
))
3311 Vector_AntiScale(row
->p
, row
->p
, g
, exist
);
3314 Matrix
*M
= unimodular_complete(row
);
3315 Matrix
*M2
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
3316 for (r
= 0; r
< nvar
; ++r
)
3317 value_set_si(M2
->p
[r
][r
], 1);
3318 for ( ; r
< nvar
+exist
; ++r
)
3319 Vector_Copy(M
->p
[r
-nvar
], M2
->p
[r
]+nvar
, exist
);
3320 for ( ; r
< P
->Dimension
+1; ++r
)
3321 value_set_si(M2
->p
[r
][r
], 1);
3322 Polyhedron
*T
= Polyhedron_Image(P
, M2
, MaxRays
);
3331 /* Split polyhedron P into two polyhedra *pos and *neg, where
3332 * existential variable i has at most one solution for each
3333 * value of the other variables in *neg.
3335 * If independent is set, then the two constraints on which the
3336 * split will be performed need to be independent of the other
3337 * existential variables.
3339 * Return true if an appropriate split could be performed.
3341 * nvar: number of set variables
3342 * exist: number of existential variables
3343 * row: temporary vector that can be used by this procedure
3344 * f: temporary value that can be used by this procedure
3346 static bool SplitOnVar(Polyhedron
*P
, int i
,
3347 int nvar
, int exist
, int MaxRays
,
3348 Vector
*row
, Value
& f
, bool independent
,
3349 Polyhedron
**pos
, Polyhedron
**neg
)
3353 for (int l
= P
->NbEq
; l
< P
->NbConstraints
; ++l
) {
3354 if (value_negz_p(P
->Constraint
[l
][nvar
+i
+1]))
3358 for (j
= 0; j
< exist
; ++j
)
3359 if (j
!= i
&& value_notzero_p(P
->Constraint
[l
][nvar
+j
+1]))
3365 for (int u
= P
->NbEq
; u
< P
->NbConstraints
; ++u
) {
3366 if (value_posz_p(P
->Constraint
[u
][nvar
+i
+1]))
3370 for (j
= 0; j
< exist
; ++j
)
3371 if (j
!= i
&& value_notzero_p(P
->Constraint
[u
][nvar
+j
+1]))
3377 if (SplitOnConstraint(P
, i
, l
, u
, nvar
, MaxRays
, row
, f
, pos
, neg
)) {
3380 SwapColumns(*neg
, nvar
+1, nvar
+1+i
);
3390 static bool double_bound_pair(Polyhedron
*P
, int nvar
, int exist
,
3391 int i
, int l1
, int l2
,
3392 Polyhedron
**pos
, Polyhedron
**neg
)
3396 Vector
*row
= Vector_Alloc(P
->Dimension
+2);
3397 value_set_si(row
->p
[0], 1);
3398 value_oppose(f
, P
->Constraint
[l1
][nvar
+i
+1]);
3399 Vector_Combine(P
->Constraint
[l1
]+1, P
->Constraint
[l2
]+1,
3401 P
->Constraint
[l2
][nvar
+i
+1], f
,
3403 ConstraintSimplify(row
->p
, row
->p
, P
->Dimension
+2, &f
);
3404 *pos
= AddConstraints(row
->p
, 1, P
, 0);
3405 value_set_si(f
, -1);
3406 Vector_Scale(row
->p
+1, row
->p
+1, f
, P
->Dimension
+1);
3407 value_decrement(row
->p
[P
->Dimension
+1], row
->p
[P
->Dimension
+1]);
3408 *neg
= AddConstraints(row
->p
, 1, P
, 0);
3412 return !emptyQ((*pos
)) && !emptyQ((*neg
));
3415 static bool double_bound(Polyhedron
*P
, int nvar
, int exist
,
3416 Polyhedron
**pos
, Polyhedron
**neg
)
3418 for (int i
= 0; i
< exist
; ++i
) {
3420 for (l1
= P
->NbEq
; l1
< P
->NbConstraints
; ++l1
) {
3421 if (value_negz_p(P
->Constraint
[l1
][nvar
+i
+1]))
3423 for (l2
= l1
+ 1; l2
< P
->NbConstraints
; ++l2
) {
3424 if (value_negz_p(P
->Constraint
[l2
][nvar
+i
+1]))
3426 if (double_bound_pair(P
, nvar
, exist
, i
, l1
, l2
, pos
, neg
))
3430 for (l1
= P
->NbEq
; l1
< P
->NbConstraints
; ++l1
) {
3431 if (value_posz_p(P
->Constraint
[l1
][nvar
+i
+1]))
3433 if (l1
< P
->NbConstraints
)
3434 for (l2
= l1
+ 1; l2
< P
->NbConstraints
; ++l2
) {
3435 if (value_posz_p(P
->Constraint
[l2
][nvar
+i
+1]))
3437 if (double_bound_pair(P
, nvar
, exist
, i
, l1
, l2
, pos
, neg
))
3449 INDEPENDENT
= 1 << 2,
3453 static evalue
* enumerate_or(Polyhedron
*D
,
3454 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3457 fprintf(stderr
, "\nER: Or\n");
3458 #endif /* DEBUG_ER */
3460 Polyhedron
*N
= D
->next
;
3463 barvinok_enumerate_e(D
, exist
, nparam
, MaxRays
);
3466 for (D
= N
; D
; D
= N
) {
3471 barvinok_enumerate_e(D
, exist
, nparam
, MaxRays
);
3474 free_evalue_refs(EN
);
3484 static evalue
* enumerate_sum(Polyhedron
*P
,
3485 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3487 int nvar
= P
->Dimension
- exist
- nparam
;
3488 int toswap
= nvar
< exist
? nvar
: exist
;
3489 for (int i
= 0; i
< toswap
; ++i
)
3490 SwapColumns(P
, 1 + i
, nvar
+exist
- i
);
3494 fprintf(stderr
, "\nER: Sum\n");
3495 #endif /* DEBUG_ER */
3497 evalue
*EP
= barvinok_enumerate_e(P
, exist
, nparam
, MaxRays
);
3499 for (int i
= 0; i
< /* nvar */ nparam
; ++i
) {
3500 Matrix
*C
= Matrix_Alloc(1, 1 + nparam
+ 1);
3501 value_set_si(C
->p
[0][0], 1);
3503 value_init(split
.d
);
3504 value_set_si(split
.d
, 0);
3505 split
.x
.p
= new_enode(partition
, 4, nparam
);
3506 value_set_si(C
->p
[0][1+i
], 1);
3507 Matrix
*C2
= Matrix_Copy(C
);
3508 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[0],
3509 Constraints2Polyhedron(C2
, MaxRays
));
3511 evalue_set_si(&split
.x
.p
->arr
[1], 1, 1);
3512 value_set_si(C
->p
[0][1+i
], -1);
3513 value_set_si(C
->p
[0][1+nparam
], -1);
3514 EVALUE_SET_DOMAIN(split
.x
.p
->arr
[2],
3515 Constraints2Polyhedron(C
, MaxRays
));
3516 evalue_set_si(&split
.x
.p
->arr
[3], 1, 1);
3518 free_evalue_refs(&split
);
3522 evalue_range_reduction(EP
);
3524 evalue_frac2floor(EP
);
3526 evalue
*sum
= esum(EP
, nvar
);
3528 free_evalue_refs(EP
);
3532 evalue_range_reduction(EP
);
3537 static evalue
* split_sure(Polyhedron
*P
, Polyhedron
*S
,
3538 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3540 int nvar
= P
->Dimension
- exist
- nparam
;
3542 Matrix
*M
= Matrix_Alloc(exist
, S
->Dimension
+2);
3543 for (int i
= 0; i
< exist
; ++i
)
3544 value_set_si(M
->p
[i
][nvar
+i
+1], 1);
3546 S
= DomainAddRays(S
, M
, MaxRays
);
3548 Polyhedron
*F
= DomainAddRays(P
, M
, MaxRays
);
3549 Polyhedron
*D
= DomainDifference(F
, S
, MaxRays
);
3551 D
= Disjoint_Domain(D
, 0, MaxRays
);
3556 M
= Matrix_Alloc(P
->Dimension
+1-exist
, P
->Dimension
+1);
3557 for (int j
= 0; j
< nvar
; ++j
)
3558 value_set_si(M
->p
[j
][j
], 1);
3559 for (int j
= 0; j
< nparam
+1; ++j
)
3560 value_set_si(M
->p
[nvar
+j
][nvar
+exist
+j
], 1);
3561 Polyhedron
*T
= Polyhedron_Image(S
, M
, MaxRays
);
3562 evalue
*EP
= barvinok_enumerate_e(T
, 0, nparam
, MaxRays
);
3567 for (Polyhedron
*Q
= D
; Q
; Q
= Q
->next
) {
3568 Polyhedron
*N
= Q
->next
;
3570 T
= DomainIntersection(P
, Q
, MaxRays
);
3571 evalue
*E
= barvinok_enumerate_e(T
, exist
, nparam
, MaxRays
);
3573 free_evalue_refs(E
);
3582 static evalue
* enumerate_sure(Polyhedron
*P
,
3583 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3587 int nvar
= P
->Dimension
- exist
- nparam
;
3593 for (i
= 0; i
< exist
; ++i
) {
3594 Matrix
*M
= Matrix_Alloc(S
->NbConstraints
, S
->Dimension
+2);
3596 value_set_si(lcm
, 1);
3597 for (int j
= 0; j
< S
->NbConstraints
; ++j
) {
3598 if (value_negz_p(S
->Constraint
[j
][1+nvar
+i
]))
3600 if (value_one_p(S
->Constraint
[j
][1+nvar
+i
]))
3602 value_lcm(lcm
, S
->Constraint
[j
][1+nvar
+i
], &lcm
);
3605 for (int j
= 0; j
< S
->NbConstraints
; ++j
) {
3606 if (value_negz_p(S
->Constraint
[j
][1+nvar
+i
]))
3608 if (value_one_p(S
->Constraint
[j
][1+nvar
+i
]))
3610 value_division(f
, lcm
, S
->Constraint
[j
][1+nvar
+i
]);
3611 Vector_Scale(S
->Constraint
[j
], M
->p
[c
], f
, S
->Dimension
+2);
3612 value_subtract(M
->p
[c
][S
->Dimension
+1],
3613 M
->p
[c
][S
->Dimension
+1],
3615 value_increment(M
->p
[c
][S
->Dimension
+1],
3616 M
->p
[c
][S
->Dimension
+1]);
3620 S
= AddConstraints(M
->p
[0], c
, S
, MaxRays
);
3635 fprintf(stderr
, "\nER: Sure\n");
3636 #endif /* DEBUG_ER */
3638 return split_sure(P
, S
, exist
, nparam
, MaxRays
);
3641 static evalue
* enumerate_sure2(Polyhedron
*P
,
3642 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3644 int nvar
= P
->Dimension
- exist
- nparam
;
3646 for (r
= 0; r
< P
->NbRays
; ++r
)
3647 if (value_one_p(P
->Ray
[r
][0]) &&
3648 value_one_p(P
->Ray
[r
][P
->Dimension
+1]))
3654 Matrix
*M
= Matrix_Alloc(nvar
+ 1 + nparam
, P
->Dimension
+2);
3655 for (int i
= 0; i
< nvar
; ++i
)
3656 value_set_si(M
->p
[i
][1+i
], 1);
3657 for (int i
= 0; i
< nparam
; ++i
)
3658 value_set_si(M
->p
[i
+nvar
][1+nvar
+exist
+i
], 1);
3659 Vector_Copy(P
->Ray
[r
]+1+nvar
, M
->p
[nvar
+nparam
]+1+nvar
, exist
);
3660 value_set_si(M
->p
[nvar
+nparam
][0], 1);
3661 value_set_si(M
->p
[nvar
+nparam
][P
->Dimension
+1], 1);
3662 Polyhedron
* F
= Rays2Polyhedron(M
, MaxRays
);
3665 Polyhedron
*I
= DomainIntersection(F
, P
, MaxRays
);
3669 fprintf(stderr
, "\nER: Sure2\n");
3670 #endif /* DEBUG_ER */
3672 return split_sure(P
, I
, exist
, nparam
, MaxRays
);
3675 static evalue
* enumerate_cyclic(Polyhedron
*P
,
3676 unsigned exist
, unsigned nparam
,
3677 evalue
* EP
, int r
, int p
, unsigned MaxRays
)
3679 int nvar
= P
->Dimension
- exist
- nparam
;
3681 /* If EP in its fractional maps only contains references
3682 * to the remainder parameter with appropriate coefficients
3683 * then we could in principle avoid adding existentially
3684 * quantified variables to the validity domains.
3685 * We'd have to replace the remainder by m { p/m }
3686 * and multiply with an appropriate factor that is one
3687 * only in the appropriate range.
3688 * This last multiplication can be avoided if EP
3689 * has a single validity domain with no (further)
3690 * constraints on the remainder parameter
3693 Matrix
*CT
= Matrix_Alloc(nparam
+1, nparam
+3);
3694 Matrix
*M
= Matrix_Alloc(1, 1+nparam
+3);
3695 for (int j
= 0; j
< nparam
; ++j
)
3697 value_set_si(CT
->p
[j
][j
], 1);
3698 value_set_si(CT
->p
[p
][nparam
+1], 1);
3699 value_set_si(CT
->p
[nparam
][nparam
+2], 1);
3700 value_set_si(M
->p
[0][1+p
], -1);
3701 value_absolute(M
->p
[0][1+nparam
], P
->Ray
[0][1+nvar
+exist
+p
]);
3702 value_set_si(M
->p
[0][1+nparam
+1], 1);
3703 Polyhedron
*CEq
= Constraints2Polyhedron(M
, 1);
3705 addeliminatedparams_enum(EP
, CT
, CEq
, MaxRays
, nparam
);
3706 Polyhedron_Free(CEq
);
3712 static void enumerate_vd_add_ray(evalue
*EP
, Matrix
*Rays
, unsigned MaxRays
)
3714 if (value_notzero_p(EP
->d
))
3717 assert(EP
->x
.p
->type
== partition
);
3718 assert(EP
->x
.p
->pos
== EVALUE_DOMAIN(EP
->x
.p
->arr
[0])->Dimension
);
3719 for (int i
= 0; i
< EP
->x
.p
->size
/2; ++i
) {
3720 Polyhedron
*D
= EVALUE_DOMAIN(EP
->x
.p
->arr
[2*i
]);
3721 Polyhedron
*N
= DomainAddRays(D
, Rays
, MaxRays
);
3722 EVALUE_SET_DOMAIN(EP
->x
.p
->arr
[2*i
], N
);
3727 static evalue
* enumerate_line(Polyhedron
*P
,
3728 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3734 fprintf(stderr
, "\nER: Line\n");
3735 #endif /* DEBUG_ER */
3737 int nvar
= P
->Dimension
- exist
- nparam
;
3739 for (i
= 0; i
< nparam
; ++i
)
3740 if (value_notzero_p(P
->Ray
[0][1+nvar
+exist
+i
]))
3743 for (j
= i
+1; j
< nparam
; ++j
)
3744 if (value_notzero_p(P
->Ray
[0][1+nvar
+exist
+i
]))
3746 assert(j
>= nparam
); // for now
3748 Matrix
*M
= Matrix_Alloc(2, P
->Dimension
+2);
3749 value_set_si(M
->p
[0][0], 1);
3750 value_set_si(M
->p
[0][1+nvar
+exist
+i
], 1);
3751 value_set_si(M
->p
[1][0], 1);
3752 value_set_si(M
->p
[1][1+nvar
+exist
+i
], -1);
3753 value_absolute(M
->p
[1][1+P
->Dimension
], P
->Ray
[0][1+nvar
+exist
+i
]);
3754 value_decrement(M
->p
[1][1+P
->Dimension
], M
->p
[1][1+P
->Dimension
]);
3755 Polyhedron
*S
= AddConstraints(M
->p
[0], 2, P
, MaxRays
);
3756 evalue
*EP
= barvinok_enumerate_e(S
, exist
, nparam
, MaxRays
);
3760 return enumerate_cyclic(P
, exist
, nparam
, EP
, 0, i
, MaxRays
);
3763 static int single_param_pos(Polyhedron
*P
, unsigned exist
, unsigned nparam
,
3766 int nvar
= P
->Dimension
- exist
- nparam
;
3767 if (First_Non_Zero(P
->Ray
[r
]+1, nvar
) != -1)
3769 int i
= First_Non_Zero(P
->Ray
[r
]+1+nvar
+exist
, nparam
);
3772 if (First_Non_Zero(P
->Ray
[r
]+1+nvar
+exist
+1, nparam
-i
-1) != -1)
3777 static evalue
* enumerate_remove_ray(Polyhedron
*P
, int r
,
3778 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3781 fprintf(stderr
, "\nER: RedundantRay\n");
3782 #endif /* DEBUG_ER */
3786 value_set_si(one
, 1);
3787 int len
= P
->NbRays
-1;
3788 Matrix
*M
= Matrix_Alloc(2 * len
, P
->Dimension
+2);
3789 Vector_Copy(P
->Ray
[0], M
->p
[0], r
* (P
->Dimension
+2));
3790 Vector_Copy(P
->Ray
[r
+1], M
->p
[r
], (len
-r
) * (P
->Dimension
+2));
3791 for (int j
= 0; j
< P
->NbRays
; ++j
) {
3794 Vector_Combine(P
->Ray
[j
], P
->Ray
[r
], M
->p
[len
+j
-(j
>r
)],
3795 one
, P
->Ray
[j
][P
->Dimension
+1], P
->Dimension
+2);
3798 P
= Rays2Polyhedron(M
, MaxRays
);
3800 evalue
*EP
= barvinok_enumerate_e(P
, exist
, nparam
, MaxRays
);
3807 static evalue
* enumerate_redundant_ray(Polyhedron
*P
,
3808 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3810 assert(P
->NbBid
== 0);
3811 int nvar
= P
->Dimension
- exist
- nparam
;
3815 for (int r
= 0; r
< P
->NbRays
; ++r
) {
3816 if (value_notzero_p(P
->Ray
[r
][P
->Dimension
+1]))
3818 int i1
= single_param_pos(P
, exist
, nparam
, r
);
3821 for (int r2
= r
+1; r2
< P
->NbRays
; ++r2
) {
3822 if (value_notzero_p(P
->Ray
[r2
][P
->Dimension
+1]))
3824 int i2
= single_param_pos(P
, exist
, nparam
, r2
);
3830 value_division(m
, P
->Ray
[r
][1+nvar
+exist
+i1
],
3831 P
->Ray
[r2
][1+nvar
+exist
+i1
]);
3832 value_multiply(m
, m
, P
->Ray
[r2
][1+nvar
+exist
+i1
]);
3833 /* r2 divides r => r redundant */
3834 if (value_eq(m
, P
->Ray
[r
][1+nvar
+exist
+i1
])) {
3836 return enumerate_remove_ray(P
, r
, exist
, nparam
, MaxRays
);
3839 value_division(m
, P
->Ray
[r2
][1+nvar
+exist
+i1
],
3840 P
->Ray
[r
][1+nvar
+exist
+i1
]);
3841 value_multiply(m
, m
, P
->Ray
[r
][1+nvar
+exist
+i1
]);
3842 /* r divides r2 => r2 redundant */
3843 if (value_eq(m
, P
->Ray
[r2
][1+nvar
+exist
+i1
])) {
3845 return enumerate_remove_ray(P
, r2
, exist
, nparam
, MaxRays
);
3853 static Polyhedron
*upper_bound(Polyhedron
*P
,
3854 int pos
, Value
*max
, Polyhedron
**R
)
3863 for (Polyhedron
*Q
= P
; Q
; Q
= N
) {
3865 for (r
= 0; r
< P
->NbRays
; ++r
) {
3866 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]) &&
3867 value_pos_p(P
->Ray
[r
][1+pos
]))
3870 if (r
< P
->NbRays
) {
3878 for (r
= 0; r
< P
->NbRays
; ++r
) {
3879 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
3881 mpz_fdiv_q(v
, P
->Ray
[r
][1+pos
], P
->Ray
[r
][1+P
->Dimension
]);
3882 if ((!Q
->next
&& r
== 0) || value_gt(v
, *max
))
3883 value_assign(*max
, v
);
3890 static evalue
* enumerate_ray(Polyhedron
*P
,
3891 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3893 assert(P
->NbBid
== 0);
3894 int nvar
= P
->Dimension
- exist
- nparam
;
3897 for (r
= 0; r
< P
->NbRays
; ++r
)
3898 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
3904 for (r2
= r
+1; r2
< P
->NbRays
; ++r2
)
3905 if (value_zero_p(P
->Ray
[r2
][P
->Dimension
+1]))
3907 if (r2
< P
->NbRays
) {
3909 return enumerate_sum(P
, exist
, nparam
, MaxRays
);
3913 fprintf(stderr
, "\nER: Ray\n");
3914 #endif /* DEBUG_ER */
3920 value_set_si(one
, 1);
3921 int i
= single_param_pos(P
, exist
, nparam
, r
);
3922 assert(i
!= -1); // for now;
3924 Matrix
*M
= Matrix_Alloc(P
->NbRays
, P
->Dimension
+2);
3925 for (int j
= 0; j
< P
->NbRays
; ++j
) {
3926 Vector_Combine(P
->Ray
[j
], P
->Ray
[r
], M
->p
[j
],
3927 one
, P
->Ray
[j
][P
->Dimension
+1], P
->Dimension
+2);
3929 Polyhedron
*S
= Rays2Polyhedron(M
, MaxRays
);
3931 Polyhedron
*D
= DomainDifference(P
, S
, MaxRays
);
3933 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
3934 assert(value_pos_p(P
->Ray
[r
][1+nvar
+exist
+i
])); // for now
3936 D
= upper_bound(D
, nvar
+exist
+i
, &m
, &R
);
3940 M
= Matrix_Alloc(2, P
->Dimension
+2);
3941 value_set_si(M
->p
[0][0], 1);
3942 value_set_si(M
->p
[1][0], 1);
3943 value_set_si(M
->p
[0][1+nvar
+exist
+i
], -1);
3944 value_set_si(M
->p
[1][1+nvar
+exist
+i
], 1);
3945 value_assign(M
->p
[0][1+P
->Dimension
], m
);
3946 value_oppose(M
->p
[1][1+P
->Dimension
], m
);
3947 value_addto(M
->p
[1][1+P
->Dimension
], M
->p
[1][1+P
->Dimension
],
3948 P
->Ray
[r
][1+nvar
+exist
+i
]);
3949 value_decrement(M
->p
[1][1+P
->Dimension
], M
->p
[1][1+P
->Dimension
]);
3950 // Matrix_Print(stderr, P_VALUE_FMT, M);
3951 D
= AddConstraints(M
->p
[0], 2, P
, MaxRays
);
3952 // Polyhedron_Print(stderr, P_VALUE_FMT, D);
3953 value_subtract(M
->p
[0][1+P
->Dimension
], M
->p
[0][1+P
->Dimension
],
3954 P
->Ray
[r
][1+nvar
+exist
+i
]);
3955 // Matrix_Print(stderr, P_VALUE_FMT, M);
3956 S
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
3957 // Polyhedron_Print(stderr, P_VALUE_FMT, S);
3960 evalue
*EP
= barvinok_enumerate_e(D
, exist
, nparam
, MaxRays
);
3965 if (value_notone_p(P
->Ray
[r
][1+nvar
+exist
+i
]))
3966 EP
= enumerate_cyclic(P
, exist
, nparam
, EP
, r
, i
, MaxRays
);
3968 M
= Matrix_Alloc(1, nparam
+2);
3969 value_set_si(M
->p
[0][0], 1);
3970 value_set_si(M
->p
[0][1+i
], 1);
3971 enumerate_vd_add_ray(EP
, M
, MaxRays
);
3976 evalue
*E
= barvinok_enumerate_e(S
, exist
, nparam
, MaxRays
);
3978 free_evalue_refs(E
);
3985 evalue
*ER
= enumerate_or(R
, exist
, nparam
, MaxRays
);
3987 free_evalue_refs(ER
);
3994 static evalue
* enumerate_vd(Polyhedron
**PA
,
3995 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
3997 Polyhedron
*P
= *PA
;
3998 int nvar
= P
->Dimension
- exist
- nparam
;
3999 Param_Polyhedron
*PP
= NULL
;
4000 Polyhedron
*C
= Universe_Polyhedron(nparam
);
4004 PP
= Polyhedron2Param_SimplifiedDomain(&PR
,C
,MaxRays
,&CEq
,&CT
);
4008 Param_Domain
*D
, *last
;
4011 for (nd
= 0, D
=PP
->D
; D
; D
=D
->next
, ++nd
)
4014 Polyhedron
**VD
= new Polyhedron_p
[nd
];
4015 Polyhedron
**fVD
= new Polyhedron_p
[nd
];
4016 for(nd
= 0, D
=PP
->D
; D
; D
=D
->next
) {
4017 Polyhedron
*rVD
= reduce_domain(D
->Domain
, CT
, CEq
,
4031 /* This doesn't seem to have any effect */
4033 Polyhedron
*CA
= align_context(VD
[0], P
->Dimension
, MaxRays
);
4035 P
= DomainIntersection(P
, CA
, MaxRays
);
4038 Polyhedron_Free(CA
);
4043 if (!EP
&& CT
->NbColumns
!= CT
->NbRows
) {
4044 Polyhedron
*CEqr
= DomainImage(CEq
, CT
, MaxRays
);
4045 Polyhedron
*CA
= align_context(CEqr
, PR
->Dimension
, MaxRays
);
4046 Polyhedron
*I
= DomainIntersection(PR
, CA
, MaxRays
);
4047 Polyhedron_Free(CEqr
);
4048 Polyhedron_Free(CA
);
4050 fprintf(stderr
, "\nER: Eliminate\n");
4051 #endif /* DEBUG_ER */
4052 nparam
-= CT
->NbColumns
- CT
->NbRows
;
4053 EP
= barvinok_enumerate_e(I
, exist
, nparam
, MaxRays
);
4054 nparam
+= CT
->NbColumns
- CT
->NbRows
;
4055 addeliminatedparams_enum(EP
, CT
, CEq
, MaxRays
, nparam
);
4059 Polyhedron_Free(PR
);
4062 if (!EP
&& nd
> 1) {
4064 fprintf(stderr
, "\nER: VD\n");
4065 #endif /* DEBUG_ER */
4066 for (int i
= 0; i
< nd
; ++i
) {
4067 Polyhedron
*CA
= align_context(VD
[i
], P
->Dimension
, MaxRays
);
4068 Polyhedron
*I
= DomainIntersection(P
, CA
, MaxRays
);
4071 EP
= barvinok_enumerate_e(I
, exist
, nparam
, MaxRays
);
4073 evalue
*E
= barvinok_enumerate_e(I
, exist
, nparam
, MaxRays
);
4075 free_evalue_refs(E
);
4079 Polyhedron_Free(CA
);
4083 for (int i
= 0; i
< nd
; ++i
) {
4084 Polyhedron_Free(VD
[i
]);
4085 Polyhedron_Free(fVD
[i
]);
4091 if (!EP
&& nvar
== 0) {
4094 Param_Vertices
*V
, *V2
;
4095 Matrix
* M
= Matrix_Alloc(1, P
->Dimension
+2);
4097 FORALL_PVertex_in_ParamPolyhedron(V
, last
, PP
) {
4099 FORALL_PVertex_in_ParamPolyhedron(V2
, last
, PP
) {
4106 for (int i
= 0; i
< exist
; ++i
) {
4107 value_oppose(f
, V
->Vertex
->p
[i
][nparam
+1]);
4108 Vector_Combine(V
->Vertex
->p
[i
],
4110 M
->p
[0] + 1 + nvar
+ exist
,
4111 V2
->Vertex
->p
[i
][nparam
+1],
4115 for (j
= 0; j
< nparam
; ++j
)
4116 if (value_notzero_p(M
->p
[0][1+nvar
+exist
+j
]))
4120 ConstraintSimplify(M
->p
[0], M
->p
[0],
4121 P
->Dimension
+2, &f
);
4122 value_set_si(M
->p
[0][0], 0);
4123 Polyhedron
*para
= AddConstraints(M
->p
[0], 1, P
,
4126 Polyhedron_Free(para
);
4129 Polyhedron
*pos
, *neg
;
4130 value_set_si(M
->p
[0][0], 1);
4131 value_decrement(M
->p
[0][P
->Dimension
+1],
4132 M
->p
[0][P
->Dimension
+1]);
4133 neg
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4134 value_set_si(f
, -1);
4135 Vector_Scale(M
->p
[0]+1, M
->p
[0]+1, f
,
4137 value_decrement(M
->p
[0][P
->Dimension
+1],
4138 M
->p
[0][P
->Dimension
+1]);
4139 value_decrement(M
->p
[0][P
->Dimension
+1],
4140 M
->p
[0][P
->Dimension
+1]);
4141 pos
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4142 if (emptyQ(neg
) && emptyQ(pos
)) {
4143 Polyhedron_Free(para
);
4144 Polyhedron_Free(pos
);
4145 Polyhedron_Free(neg
);
4149 fprintf(stderr
, "\nER: Order\n");
4150 #endif /* DEBUG_ER */
4151 EP
= barvinok_enumerate_e(para
, exist
, nparam
, MaxRays
);
4154 E
= barvinok_enumerate_e(pos
, exist
, nparam
, MaxRays
);
4156 free_evalue_refs(E
);
4160 E
= barvinok_enumerate_e(neg
, exist
, nparam
, MaxRays
);
4162 free_evalue_refs(E
);
4165 Polyhedron_Free(para
);
4166 Polyhedron_Free(pos
);
4167 Polyhedron_Free(neg
);
4172 } END_FORALL_PVertex_in_ParamPolyhedron
;
4175 } END_FORALL_PVertex_in_ParamPolyhedron
;
4178 /* Search for vertex coordinate to split on */
4179 /* First look for one independent of the parameters */
4180 FORALL_PVertex_in_ParamPolyhedron(V
, last
, PP
) {
4181 for (int i
= 0; i
< exist
; ++i
) {
4183 for (j
= 0; j
< nparam
; ++j
)
4184 if (value_notzero_p(V
->Vertex
->p
[i
][j
]))
4188 value_set_si(M
->p
[0][0], 1);
4189 Vector_Set(M
->p
[0]+1, 0, nvar
+exist
);
4190 Vector_Copy(V
->Vertex
->p
[i
],
4191 M
->p
[0] + 1 + nvar
+ exist
, nparam
+1);
4192 value_oppose(M
->p
[0][1+nvar
+i
],
4193 V
->Vertex
->p
[i
][nparam
+1]);
4195 Polyhedron
*pos
, *neg
;
4196 value_set_si(M
->p
[0][0], 1);
4197 value_decrement(M
->p
[0][P
->Dimension
+1],
4198 M
->p
[0][P
->Dimension
+1]);
4199 neg
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4200 value_set_si(f
, -1);
4201 Vector_Scale(M
->p
[0]+1, M
->p
[0]+1, f
,
4203 value_decrement(M
->p
[0][P
->Dimension
+1],
4204 M
->p
[0][P
->Dimension
+1]);
4205 value_decrement(M
->p
[0][P
->Dimension
+1],
4206 M
->p
[0][P
->Dimension
+1]);
4207 pos
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4208 if (emptyQ(neg
) || emptyQ(pos
)) {
4209 Polyhedron_Free(pos
);
4210 Polyhedron_Free(neg
);
4213 Polyhedron_Free(pos
);
4214 value_increment(M
->p
[0][P
->Dimension
+1],
4215 M
->p
[0][P
->Dimension
+1]);
4216 pos
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4218 fprintf(stderr
, "\nER: Vertex\n");
4219 #endif /* DEBUG_ER */
4221 EP
= enumerate_or(pos
, exist
, nparam
, MaxRays
);
4226 } END_FORALL_PVertex_in_ParamPolyhedron
;
4230 /* Search for vertex coordinate to split on */
4231 /* Now look for one that depends on the parameters */
4232 FORALL_PVertex_in_ParamPolyhedron(V
, last
, PP
) {
4233 for (int i
= 0; i
< exist
; ++i
) {
4234 value_set_si(M
->p
[0][0], 1);
4235 Vector_Set(M
->p
[0]+1, 0, nvar
+exist
);
4236 Vector_Copy(V
->Vertex
->p
[i
],
4237 M
->p
[0] + 1 + nvar
+ exist
, nparam
+1);
4238 value_oppose(M
->p
[0][1+nvar
+i
],
4239 V
->Vertex
->p
[i
][nparam
+1]);
4241 Polyhedron
*pos
, *neg
;
4242 value_set_si(M
->p
[0][0], 1);
4243 value_decrement(M
->p
[0][P
->Dimension
+1],
4244 M
->p
[0][P
->Dimension
+1]);
4245 neg
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4246 value_set_si(f
, -1);
4247 Vector_Scale(M
->p
[0]+1, M
->p
[0]+1, f
,
4249 value_decrement(M
->p
[0][P
->Dimension
+1],
4250 M
->p
[0][P
->Dimension
+1]);
4251 value_decrement(M
->p
[0][P
->Dimension
+1],
4252 M
->p
[0][P
->Dimension
+1]);
4253 pos
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4254 if (emptyQ(neg
) || emptyQ(pos
)) {
4255 Polyhedron_Free(pos
);
4256 Polyhedron_Free(neg
);
4259 Polyhedron_Free(pos
);
4260 value_increment(M
->p
[0][P
->Dimension
+1],
4261 M
->p
[0][P
->Dimension
+1]);
4262 pos
= AddConstraints(M
->p
[0], 1, P
, MaxRays
);
4264 fprintf(stderr
, "\nER: ParamVertex\n");
4265 #endif /* DEBUG_ER */
4267 EP
= enumerate_or(pos
, exist
, nparam
, MaxRays
);
4272 } END_FORALL_PVertex_in_ParamPolyhedron
;
4280 Polyhedron_Free(CEq
);
4284 Param_Polyhedron_Free(PP
);
4291 evalue
*barvinok_enumerate_pip(Polyhedron
*P
,
4292 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
4297 evalue
*barvinok_enumerate_pip(Polyhedron
*P
,
4298 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
4300 int nvar
= P
->Dimension
- exist
- nparam
;
4301 evalue
*EP
= evalue_zero();
4305 fprintf(stderr
, "\nER: PIP\n");
4306 #endif /* DEBUG_ER */
4308 Polyhedron
*D
= pip_projectout(P
, nvar
, exist
, nparam
);
4309 for (Q
= D
; Q
; Q
= N
) {
4313 exist
= Q
->Dimension
- nvar
- nparam
;
4314 E
= barvinok_enumerate_e(Q
, exist
, nparam
, MaxRays
);
4317 free_evalue_refs(E
);
4326 static bool is_single(Value
*row
, int pos
, int len
)
4328 return First_Non_Zero(row
, pos
) == -1 &&
4329 First_Non_Zero(row
+pos
+1, len
-pos
-1) == -1;
4332 static evalue
* barvinok_enumerate_e_r(Polyhedron
*P
,
4333 unsigned exist
, unsigned nparam
, unsigned MaxRays
);
4336 static int er_level
= 0;
4338 evalue
* barvinok_enumerate_e(Polyhedron
*P
,
4339 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
4341 fprintf(stderr
, "\nER: level %i\n", er_level
);
4343 Polyhedron_PrintConstraints(stderr
, P_VALUE_FMT
, P
);
4344 fprintf(stderr
, "\nE %d\nP %d\n", exist
, nparam
);
4346 P
= DomainConstraintSimplify(Polyhedron_Copy(P
), MaxRays
);
4347 evalue
*EP
= barvinok_enumerate_e_r(P
, exist
, nparam
, MaxRays
);
4353 evalue
* barvinok_enumerate_e(Polyhedron
*P
,
4354 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
4356 P
= DomainConstraintSimplify(Polyhedron_Copy(P
), MaxRays
);
4357 evalue
*EP
= barvinok_enumerate_e_r(P
, exist
, nparam
, MaxRays
);
4363 static evalue
* barvinok_enumerate_e_r(Polyhedron
*P
,
4364 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
4367 Polyhedron
*U
= Universe_Polyhedron(nparam
);
4368 evalue
*EP
= barvinok_enumerate_ev(P
, U
, MaxRays
);
4369 //char *param_name[] = {"P", "Q", "R", "S", "T" };
4370 //print_evalue(stdout, EP, param_name);
4375 int nvar
= P
->Dimension
- exist
- nparam
;
4376 int len
= P
->Dimension
+ 2;
4379 POL_ENSURE_FACETS(P
);
4380 POL_ENSURE_VERTICES(P
);
4383 return evalue_zero();
4385 if (nvar
== 0 && nparam
== 0) {
4386 evalue
*EP
= evalue_zero();
4387 barvinok_count(P
, &EP
->x
.n
, MaxRays
);
4388 if (value_pos_p(EP
->x
.n
))
4389 value_set_si(EP
->x
.n
, 1);
4394 for (r
= 0; r
< P
->NbRays
; ++r
)
4395 if (value_zero_p(P
->Ray
[r
][0]) ||
4396 value_zero_p(P
->Ray
[r
][P
->Dimension
+1])) {
4398 for (i
= 0; i
< nvar
; ++i
)
4399 if (value_notzero_p(P
->Ray
[r
][i
+1]))
4403 for (i
= nvar
+ exist
; i
< nvar
+ exist
+ nparam
; ++i
)
4404 if (value_notzero_p(P
->Ray
[r
][i
+1]))
4406 if (i
>= nvar
+ exist
+ nparam
)
4409 if (r
< P
->NbRays
) {
4410 evalue
*EP
= evalue_zero();
4411 value_set_si(EP
->x
.n
, -1);
4416 for (r
= 0; r
< P
->NbEq
; ++r
)
4417 if ((first
= First_Non_Zero(P
->Constraint
[r
]+1+nvar
, exist
)) != -1)
4420 if (First_Non_Zero(P
->Constraint
[r
]+1+nvar
+first
+1,
4421 exist
-first
-1) != -1) {
4422 Polyhedron
*T
= rotate_along(P
, r
, nvar
, exist
, MaxRays
);
4424 fprintf(stderr
, "\nER: Equality\n");
4425 #endif /* DEBUG_ER */
4426 evalue
*EP
= barvinok_enumerate_e(T
, exist
-1, nparam
, MaxRays
);
4431 fprintf(stderr
, "\nER: Fixed\n");
4432 #endif /* DEBUG_ER */
4434 return barvinok_enumerate_e(P
, exist
-1, nparam
, MaxRays
);
4436 Polyhedron
*T
= Polyhedron_Copy(P
);
4437 SwapColumns(T
, nvar
+1, nvar
+1+first
);
4438 evalue
*EP
= barvinok_enumerate_e(T
, exist
-1, nparam
, MaxRays
);
4445 Vector
*row
= Vector_Alloc(len
);
4446 value_set_si(row
->p
[0], 1);
4451 enum constraint
* info
= new constraint
[exist
];
4452 for (int i
= 0; i
< exist
; ++i
) {
4454 for (int l
= P
->NbEq
; l
< P
->NbConstraints
; ++l
) {
4455 if (value_negz_p(P
->Constraint
[l
][nvar
+i
+1]))
4457 bool l_parallel
= is_single(P
->Constraint
[l
]+nvar
+1, i
, exist
);
4458 for (int u
= P
->NbEq
; u
< P
->NbConstraints
; ++u
) {
4459 if (value_posz_p(P
->Constraint
[u
][nvar
+i
+1]))
4461 bool lu_parallel
= l_parallel
||
4462 is_single(P
->Constraint
[u
]+nvar
+1, i
, exist
);
4463 value_oppose(f
, P
->Constraint
[u
][nvar
+i
+1]);
4464 Vector_Combine(P
->Constraint
[l
]+1, P
->Constraint
[u
]+1, row
->p
+1,
4465 f
, P
->Constraint
[l
][nvar
+i
+1], len
-1);
4466 if (!(info
[i
] & INDEPENDENT
)) {
4468 for (j
= 0; j
< exist
; ++j
)
4469 if (j
!= i
&& value_notzero_p(row
->p
[nvar
+j
+1]))
4472 //printf("independent: i: %d, l: %d, u: %d\n", i, l, u);
4473 info
[i
] = (constraint
)(info
[i
] | INDEPENDENT
);
4476 if (info
[i
] & ALL_POS
) {
4477 value_addto(row
->p
[len
-1], row
->p
[len
-1],
4478 P
->Constraint
[l
][nvar
+i
+1]);
4479 value_addto(row
->p
[len
-1], row
->p
[len
-1], f
);
4480 value_multiply(f
, f
, P
->Constraint
[l
][nvar
+i
+1]);
4481 value_subtract(row
->p
[len
-1], row
->p
[len
-1], f
);
4482 value_decrement(row
->p
[len
-1], row
->p
[len
-1]);
4483 ConstraintSimplify(row
->p
, row
->p
, len
, &f
);
4484 value_set_si(f
, -1);
4485 Vector_Scale(row
->p
+1, row
->p
+1, f
, len
-1);
4486 value_decrement(row
->p
[len
-1], row
->p
[len
-1]);
4487 Polyhedron
*T
= AddConstraints(row
->p
, 1, P
, MaxRays
);
4489 //printf("not all_pos: i: %d, l: %d, u: %d\n", i, l, u);
4490 info
[i
] = (constraint
)(info
[i
] ^ ALL_POS
);
4492 //puts("pos remainder");
4493 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
4496 if (!(info
[i
] & ONE_NEG
)) {
4498 negative_test_constraint(P
->Constraint
[l
],
4500 row
->p
, nvar
+i
, len
, &f
);
4501 oppose_constraint(row
->p
, len
, &f
);
4502 Polyhedron
*T
= AddConstraints(row
->p
, 1, P
, MaxRays
);
4504 //printf("one_neg i: %d, l: %d, u: %d\n", i, l, u);
4505 info
[i
] = (constraint
)(info
[i
] | ONE_NEG
);
4507 //puts("neg remainder");
4508 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
4510 } else if (!(info
[i
] & ROT_NEG
)) {
4511 if (parallel_constraints(P
->Constraint
[l
],
4513 row
->p
, nvar
, exist
)) {
4514 negative_test_constraint7(P
->Constraint
[l
],
4516 row
->p
, nvar
, exist
,
4518 oppose_constraint(row
->p
, len
, &f
);
4519 Polyhedron
*T
= AddConstraints(row
->p
, 1, P
, MaxRays
);
4521 // printf("rot_neg i: %d, l: %d, u: %d\n", i, l, u);
4522 info
[i
] = (constraint
)(info
[i
] | ROT_NEG
);
4525 //puts("neg remainder");
4526 //Polyhedron_Print(stdout, P_VALUE_FMT, T);
4531 if (!(info
[i
] & ALL_POS
) && (info
[i
] & (ONE_NEG
| ROT_NEG
)))
4535 if (info
[i
] & ALL_POS
)
4542 for (int i = 0; i < exist; ++i)
4543 printf("%i: %i\n", i, info[i]);
4545 for (int i
= 0; i
< exist
; ++i
)
4546 if (info
[i
] & ALL_POS
) {
4548 fprintf(stderr
, "\nER: Positive\n");
4549 #endif /* DEBUG_ER */
4551 // Maybe we should chew off some of the fat here
4552 Matrix
*M
= Matrix_Alloc(P
->Dimension
, P
->Dimension
+1);
4553 for (int j
= 0; j
< P
->Dimension
; ++j
)
4554 value_set_si(M
->p
[j
][j
+ (j
>= i
+nvar
)], 1);
4555 Polyhedron
*T
= Polyhedron_Image(P
, M
, MaxRays
);
4557 evalue
*EP
= barvinok_enumerate_e(T
, exist
-1, nparam
, MaxRays
);
4564 for (int i
= 0; i
< exist
; ++i
)
4565 if (info
[i
] & ONE_NEG
) {
4567 fprintf(stderr
, "\nER: Negative\n");
4568 #endif /* DEBUG_ER */
4573 return barvinok_enumerate_e(P
, exist
-1, nparam
, MaxRays
);
4575 Polyhedron
*T
= Polyhedron_Copy(P
);
4576 SwapColumns(T
, nvar
+1, nvar
+1+i
);
4577 evalue
*EP
= barvinok_enumerate_e(T
, exist
-1, nparam
, MaxRays
);
4582 for (int i
= 0; i
< exist
; ++i
)
4583 if (info
[i
] & ROT_NEG
) {
4585 fprintf(stderr
, "\nER: Rotate\n");
4586 #endif /* DEBUG_ER */
4590 Polyhedron
*T
= rotate_along(P
, r
, nvar
, exist
, MaxRays
);
4591 evalue
*EP
= barvinok_enumerate_e(T
, exist
-1, nparam
, MaxRays
);
4595 for (int i
= 0; i
< exist
; ++i
)
4596 if (info
[i
] & INDEPENDENT
) {
4597 Polyhedron
*pos
, *neg
;
4599 /* Find constraint again and split off negative part */
4601 if (SplitOnVar(P
, i
, nvar
, exist
, MaxRays
,
4602 row
, f
, true, &pos
, &neg
)) {
4604 fprintf(stderr
, "\nER: Split\n");
4605 #endif /* DEBUG_ER */
4608 barvinok_enumerate_e(neg
, exist
-1, nparam
, MaxRays
);
4610 barvinok_enumerate_e(pos
, exist
, nparam
, MaxRays
);
4612 free_evalue_refs(E
);
4614 Polyhedron_Free(neg
);
4615 Polyhedron_Free(pos
);
4629 EP
= enumerate_line(P
, exist
, nparam
, MaxRays
);
4633 EP
= barvinok_enumerate_pip(P
, exist
, nparam
, MaxRays
);
4637 EP
= enumerate_redundant_ray(P
, exist
, nparam
, MaxRays
);
4641 EP
= enumerate_sure(P
, exist
, nparam
, MaxRays
);
4645 EP
= enumerate_ray(P
, exist
, nparam
, MaxRays
);
4649 EP
= enumerate_sure2(P
, exist
, nparam
, MaxRays
);
4653 F
= unfringe(P
, MaxRays
);
4654 if (!PolyhedronIncludes(F
, P
)) {
4656 fprintf(stderr
, "\nER: Fringed\n");
4657 #endif /* DEBUG_ER */
4658 EP
= barvinok_enumerate_e(F
, exist
, nparam
, MaxRays
);
4665 EP
= enumerate_vd(&P
, exist
, nparam
, MaxRays
);
4670 EP
= enumerate_sum(P
, exist
, nparam
, MaxRays
);
4677 Polyhedron
*pos
, *neg
;
4678 for (i
= 0; i
< exist
; ++i
)
4679 if (SplitOnVar(P
, i
, nvar
, exist
, MaxRays
,
4680 row
, f
, false, &pos
, &neg
))
4686 EP
= enumerate_or(pos
, exist
, nparam
, MaxRays
);
4698 static void split_param_compression(Matrix
*CP
, mat_ZZ
& map
, vec_ZZ
& offset
)
4700 Matrix
*T
= Transpose(CP
);
4701 matrix2zz(T
, map
, T
->NbRows
-1, T
->NbColumns
-1);
4702 values2zz(T
->p
[T
->NbRows
-1], offset
, T
->NbColumns
-1);
4707 * remove equalities that require a "compression" of the parameters
4709 #ifndef HAVE_COMPRESS_PARMS
4710 static Polyhedron
*remove_more_equalities(Polyhedron
*P
, unsigned nparam
,
4711 Matrix
**CP
, unsigned MaxRays
)
4716 static Polyhedron
*remove_more_equalities(Polyhedron
*P
, unsigned nparam
,
4717 Matrix
**CP
, unsigned MaxRays
)
4722 /* compress_parms doesn't like equalities that only involve parameters */
4723 for (int i
= 0; i
< P
->NbEq
; ++i
)
4724 assert(First_Non_Zero(P
->Constraint
[i
]+1, P
->Dimension
-nparam
) != -1);
4726 M
= Matrix_Alloc(P
->NbEq
, P
->Dimension
+2);
4727 Vector_Copy(P
->Constraint
[0], M
->p
[0], P
->NbEq
* (P
->Dimension
+2));
4728 *CP
= compress_parms(M
, nparam
);
4729 T
= align_matrix(*CP
, P
->Dimension
+1);
4730 Q
= Polyhedron_Preimage(P
, T
, MaxRays
);
4733 P
= remove_equalities_p(P
, P
->Dimension
-nparam
, NULL
);
4740 gen_fun
* barvinok_series(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
4744 unsigned nparam
= C
->Dimension
;
4746 CA
= align_context(C
, P
->Dimension
, MaxRays
);
4747 P
= DomainIntersection(P
, CA
, MaxRays
);
4748 Polyhedron_Free(CA
);
4755 assert(!Polyhedron_is_infinite_param(P
, nparam
));
4756 assert(P
->NbBid
== 0);
4757 assert(Polyhedron_has_positive_rays(P
, nparam
));
4759 P
= remove_equalities_p(P
, P
->Dimension
-nparam
, NULL
);
4761 P
= remove_more_equalities(P
, nparam
, &CP
, MaxRays
);
4762 assert(P
->NbEq
== 0);
4764 #ifdef USE_INCREMENTAL_BF
4765 partial_bfcounter
red(P
, nparam
);
4766 #elif defined USE_INCREMENTAL_DF
4767 partial_ireducer
red(P
, nparam
);
4769 partial_reducer
red(P
, nparam
);
4776 split_param_compression(CP
, map
, offset
);
4777 red
.gf
->substitute(CP
, map
, offset
);
4783 static Polyhedron
*skew_into_positive_orthant(Polyhedron
*D
, unsigned nparam
,
4789 for (Polyhedron
*P
= D
; P
; P
= P
->next
) {
4790 POL_ENSURE_VERTICES(P
);
4791 assert(!Polyhedron_is_infinite_param(P
, nparam
));
4792 assert(P
->NbBid
== 0);
4793 assert(Polyhedron_has_positive_rays(P
, nparam
));
4795 for (int r
= 0; r
< P
->NbRays
; ++r
) {
4796 if (value_notzero_p(P
->Ray
[r
][P
->Dimension
+1]))
4798 for (int i
= 0; i
< nparam
; ++i
) {
4800 if (value_posz_p(P
->Ray
[r
][i
+1]))
4803 M
= Matrix_Alloc(D
->Dimension
+1, D
->Dimension
+1);
4804 for (int i
= 0; i
< D
->Dimension
+1; ++i
)
4805 value_set_si(M
->p
[i
][i
], 1);
4807 Inner_Product(P
->Ray
[r
]+1, M
->p
[i
], D
->Dimension
+1, &tmp
);
4808 if (value_posz_p(tmp
))
4811 for (j
= P
->Dimension
- nparam
; j
< P
->Dimension
; ++j
)
4812 if (value_pos_p(P
->Ray
[r
][j
+1]))
4814 assert(j
< P
->Dimension
);
4815 value_pdivision(tmp
, P
->Ray
[r
][j
+1], P
->Ray
[r
][i
+1]);
4816 value_subtract(M
->p
[i
][j
], M
->p
[i
][j
], tmp
);
4822 D
= DomainImage(D
, M
, MaxRays
);
4828 gen_fun
* barvinok_enumerate_union_series(Polyhedron
*D
, Polyhedron
* C
,
4831 Polyhedron
*conv
, *D2
;
4833 unsigned nparam
= C
->Dimension
;
4837 D2
= skew_into_positive_orthant(D
, nparam
, MaxRays
);
4838 for (Polyhedron
*P
= D2
; P
; P
= P
->next
) {
4839 assert(P
->Dimension
== D2
->Dimension
);
4840 POL_ENSURE_VERTICES(P
);
4841 /* it doesn't matter which reducer we use, since we don't actually
4842 * reduce anything here
4844 partial_reducer
red(P
, P
->Dimension
);
4849 gf
->add_union(red
.gf
, MaxRays
);
4853 /* we actually only need the convex union of the parameter space
4854 * but the reducer classes currently expect a polyhedron in
4855 * the combined space
4857 conv
= DomainConvex(D2
, MaxRays
);
4858 #ifdef USE_INCREMENTAL_DF
4859 partial_ireducer
red(conv
, nparam
);
4861 partial_reducer
red(conv
, nparam
);
4863 for (int i
= 0; i
< gf
->term
.size(); ++i
) {
4864 for (int j
= 0; j
< gf
->term
[i
]->n
.power
.NumRows(); ++j
) {
4865 red
.reduce(gf
->term
[i
]->n
.coeff
[j
][0], gf
->term
[i
]->n
.coeff
[j
][1],
4866 gf
->term
[i
]->n
.power
[j
], gf
->term
[i
]->d
.power
);
4872 Polyhedron_Free(conv
);
4876 evalue
* barvinok_enumerate_union(Polyhedron
*D
, Polyhedron
* C
, unsigned MaxRays
)
4879 gen_fun
*gf
= barvinok_enumerate_union_series(D
, C
, MaxRays
);