4 #include <NTL/mat_ZZ.h>
6 #include <barvinok/barvinok.h>
7 #include <barvinok/util.h>
8 #include "conversion.h"
9 #include "decomposer.h"
19 * Returns the largest absolute value in the vector
21 static ZZ
max(vec_ZZ
& v
)
24 for (int i
= 1; i
< v
.length(); ++i
)
30 /* Remove common divisor of elements of cols of B */
31 static void normalize_cols(mat_ZZ
& B
)
34 for (int i
= 0; i
< B
.NumCols(); ++i
) {
36 for (int j
= 1 ; gcd
!= 1 && j
< B
.NumRows(); ++j
)
37 GCD(gcd
, gcd
, B
[j
][i
]);
39 for (int j
= 0; j
< B
.NumRows(); ++j
)
46 cone(const mat_ZZ
& r
, int row
, const vec_ZZ
& w
) {
52 cone(const signed_cone
& sc
) {
55 set_closed(sc
.closed
);
58 det
= determinant(rays
);
60 void set_closed(int *cl
) {
63 closed
= new int[rays
.NumRows()];
64 for (int i
= 0; i
< rays
.NumRows(); ++i
)
68 bool needs_split(barvinok_options
*options
) {
72 if (options
->primal
&& index
<= options
->max_index
)
75 Matrix
*M
= rays2matrix(rays
);
76 Matrix
*inv
= Matrix_Alloc(M
->NbRows
, M
->NbColumns
);
77 int ok
= Matrix_Inverse(M
, inv
);
81 matrix2zz(inv
, B
, inv
->NbRows
- 1, inv
->NbColumns
- 1);
84 if (!options
->primal
&& options
->max_index
> 1) {
87 index
= abs(determinant(B2
));
88 if (index
<= options
->max_index
)
95 void short_vector(vec_ZZ
& v
, vec_ZZ
& lambda
, barvinok_options
*options
) {
98 long r
= LLL(det2
, B
, U
, options
->LLL_a
, options
->LLL_b
);
102 for (int i
= 1; i
< B
.NumRows(); ++i
) {
115 for (i
= 0; i
< lambda
.length(); ++i
)
118 if (i
== lambda
.length()) {
136 static void decompose(const signed_cone
& sc
, signed_cone_consumer
& scc
,
137 barvinok_options
*options
)
139 vector
<cone
*> nonuni
;
140 cone
* c
= new cone(sc
);
144 if (c
->needs_split(options
)) {
148 options
->stats
->base_cones
++;
149 scc
.handle(signed_cone(sc
.C
, sc
.rays
, sc
.sign
, to_ulong(c
->index
),
150 sc
.closed
), options
);
160 int closed
[c
->rays
.NumRows()];
161 while (!nonuni
.empty()) {
164 c
->short_vector(v
, lambda
, options
);
165 for (int i
= 0; i
< c
->rays
.NumRows(); ++i
) {
168 cone
*pc
= new cone(c
->rays
, i
, v
);
170 for (int j
= 0; j
< c
->rays
.NumRows(); ++j
) {
172 closed
[j
] = c
->closed
[j
];
175 closed
[j
] = c
->closed
[j
];
177 closed
[j
] = !c
->closed
[j
];
178 } else if (sign(lambda
[i
]) == sign(lambda
[j
])) {
179 if (c
->closed
[i
] == c
->closed
[j
])
182 closed
[j
] = c
->closed
[j
];
184 closed
[j
] = c
->closed
[i
] && c
->closed
[j
];
186 pc
->set_closed(closed
);
188 assert (pc
->det
!= 0);
189 if (pc
->needs_split(options
)) {
190 assert(abs(pc
->det
) < abs(c
->det
));
191 nonuni
.push_back(pc
);
194 options
->stats
->base_cones
++;
195 scc
.handle(signed_cone(pc
->rays
, sign(pc
->det
) * s
,
197 pc
->closed
), options
);
202 while (!nonuni
.empty()) {
215 struct polar_signed_cone_consumer
: public signed_cone_consumer
{
216 signed_cone_consumer
& scc
;
218 polar_signed_cone_consumer(signed_cone_consumer
& scc
) : scc(scc
) {}
219 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
) {
220 Polyhedron
*C
= sc
.C
;
222 Matrix
*M
= Matrix_Alloc(sc
.rays
.NumRows()+1, sc
.rays
.NumCols()+2);
223 for (int i
= 0; i
< sc
.rays
.NumRows(); ++i
) {
224 zz2values(sc
.rays
[i
], M
->p
[i
]+1);
225 value_set_si(M
->p
[i
][0], 1);
227 value_set_si(M
->p
[sc
.rays
.NumRows()][0], 1);
228 value_set_si(M
->p
[sc
.rays
.NumRows()][1+sc
.rays
.NumCols()], 1);
229 C
= Rays2Polyhedron(M
, M
->NbRows
+1);
230 assert(C
->NbConstraints
== C
->NbRays
);
233 Polyhedron_Polarize(C
);
235 scc
.handle(signed_cone(C
, r
, sc
.sign
, sc
.det
), options
);
241 /* Remove common divisor of elements of rows of B */
242 static void normalize_rows(mat_ZZ
& B
)
245 for (int i
= 0; i
< B
.NumRows(); ++i
) {
247 for (int j
= 1 ; gcd
!= 1 && j
< B
.NumCols(); ++j
)
248 GCD(gcd
, gcd
, B
[i
][j
]);
250 for (int j
= 0; j
< B
.NumCols(); ++j
)
255 static void polar_decompose(Polyhedron
*cone
, signed_cone_consumer
& scc
,
256 barvinok_options
*options
)
258 POL_ENSURE_VERTICES(cone
);
259 Polyhedron_Polarize(cone
);
260 if (cone
->NbRays
- 1 != cone
->Dimension
) {
261 Polyhedron
*tmp
= cone
;
262 cone
= triangulate_cone(cone
, options
->MaxRays
);
263 Polyhedron_Free(tmp
);
265 polar_signed_cone_consumer
pssc(scc
);
268 for (Polyhedron
*Polar
= cone
; Polar
; Polar
= Polar
->next
) {
271 decompose(signed_cone(Polar
, r
, 1), pssc
, options
);
280 static void primal_decompose(Polyhedron
*cone
, signed_cone_consumer
& scc
,
281 barvinok_options
*options
)
283 POL_ENSURE_VERTICES(cone
);
285 if (cone
->NbRays
- 1 == cone
->Dimension
)
288 parts
= triangulate_cone(cone
, options
->MaxRays
);
289 int closed
[cone
->Dimension
];
290 Vector
*average
= NULL
;
294 average
= Vector_Alloc(cone
->Dimension
);
295 for (int i
= 0; i
< cone
->NbRays
; ++i
) {
296 if (value_notzero_p(cone
->Ray
[i
][1+cone
->Dimension
]))
298 Vector_Add(average
->p
, cone
->Ray
[i
]+1, average
->p
, cone
->Dimension
);
303 for (Polyhedron
*simple
= parts
; simple
; simple
= simple
->next
) {
304 for (int i
= 0, r
= 0; r
< simple
->NbRays
; ++r
) {
305 if (value_notzero_p(simple
->Ray
[r
][1+simple
->Dimension
]))
307 if (simple
== cone
) {
311 for (f
= 0; f
< simple
->NbConstraints
; ++f
) {
312 Inner_Product(simple
->Ray
[r
]+1, simple
->Constraint
[f
]+1,
313 simple
->Dimension
, &tmp
);
314 if (value_notzero_p(tmp
))
317 assert(f
< simple
->NbConstraints
);
318 Inner_Product(simple
->Constraint
[f
]+1, average
->p
,
319 simple
->Dimension
, &tmp
);
320 if (value_notzero_p(tmp
))
321 closed
[i
] = value_pos_p(tmp
);
323 int p
= First_Non_Zero(simple
->Constraint
[f
]+1,
325 closed
[i
] = value_pos_p(simple
->Constraint
[f
][1+p
]);
331 decompose(signed_cone(simple
, ray
, 1, 0, closed
), scc
, options
);
337 Vector_Free(average
);
344 Vector_Free(average
);
350 void barvinok_decompose(Polyhedron
*C
, signed_cone_consumer
& scc
,
351 barvinok_options
*options
)
354 primal_decompose(C
, scc
, options
);
356 polar_decompose(C
, scc
, options
);
359 void vertex_decomposer::decompose_at_vertex(Param_Vertices
*V
, int _i
,
360 barvinok_options
*options
)
362 Polyhedron
*C
= supporting_cone_p(P
, V
);
366 barvinok_decompose(C
, scc
, options
);
369 struct posneg_collector
: public signed_cone_consumer
{
370 posneg_collector(Polyhedron
*pos
, Polyhedron
*neg
) : pos(pos
), neg(neg
) {}
371 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
) {
372 Polyhedron
*p
= Polyhedron_Copy(sc
.C
);
386 * Barvinok's Decomposition of a simplicial cone
388 * Returns two lists of polyhedra
390 void barvinok_decompose(Polyhedron
*C
, Polyhedron
**ppos
, Polyhedron
**pneg
)
392 barvinok_options
*options
= barvinok_options_new_with_defaults();
393 posneg_collector
pc(*ppos
, *pneg
);
394 POL_ENSURE_VERTICES(C
);
397 decompose(signed_cone(C
, r
, 1), pc
, options
);
400 barvinok_options_free(options
);