2 #include <barvinok/barvinok.h>
3 #include <barvinok/util.h>
4 #include <barvinok/options.h>
6 #include "reduce_domain.h"
7 #include "param_util.h"
9 #define ALLOC(type) (type*)malloc(sizeof(type))
10 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
12 /* If a vertex is described by A x + B p + c = 0, then
13 * M = [A B] and we want to compute a linear transformation L such
14 * that H L = A and H \Z contains both A \Z and B \Z.
16 * [ A B ] = [ H 0 ] [ U_11 U_12 ]
19 * U_11 is the required linear transformation.
20 * Note that this also works if M has more rows than there are variables,
21 * i.e., if some rows in M are linear combinations of other rows.
22 * These extra rows only affect and H and not U.
24 static Lattice
*extract_lattice(Matrix
*M
, unsigned nvar
)
27 Matrix
*H
, *Q
, *U
, *Li
;
31 left_hermite(M
, &H
, &Q
, &U
);
34 Li
= Matrix_Alloc(nvar
+1, nvar
+1);
35 L
= Matrix_Alloc(nvar
+1, nvar
+1);
36 value_set_si(Li
->p
[nvar
][nvar
], 1);
38 for (row
= 0; row
< nvar
; ++row
)
39 Vector_Copy(Q
->p
[row
], Li
->p
[row
], nvar
);
43 ok
= Matrix_Inverse(Li
, L
);
50 /* Returns the smallest (wrt inclusion) lattice that contains both X and Y */
51 static Lattice
*LatticeJoin(Lattice
*X
, Lattice
*Y
)
54 int dim
= X
->NbRows
-1;
58 Matrix
*M
, *H
, *U
, *Q
;
60 assert(X
->NbColumns
-1 == dim
);
61 assert(Y
->NbRows
-1 == dim
);
62 assert(Y
->NbColumns
-1 == dim
);
67 M
= Matrix_Alloc(dim
, 2*dim
);
68 value_lcm(X
->p
[dim
][dim
], Y
->p
[dim
][dim
], &lcm
);
70 value_division(tmp
, lcm
, X
->p
[dim
][dim
]);
71 for (i
= 0; i
< dim
; ++i
)
72 Vector_Scale(X
->p
[i
], M
->p
[i
], tmp
, dim
);
73 value_division(tmp
, lcm
, Y
->p
[dim
][dim
]);
74 for (i
= 0; i
< dim
; ++i
)
75 Vector_Scale(Y
->p
[i
], M
->p
[i
]+dim
, tmp
, dim
);
77 left_hermite(M
, &H
, &Q
, &U
);
82 L
= Matrix_Alloc(dim
+1, dim
+1);
83 value_assign(L
->p
[dim
][dim
], lcm
);
84 for (i
= 0; i
< dim
; ++i
)
85 Vector_Copy(H
->p
[i
], L
->p
[i
], dim
);
93 static void Param_Vertex_Image(Param_Vertices
*V
, Matrix
*T
)
95 unsigned nvar
= V
->Vertex
->NbRows
;
96 unsigned nparam
= V
->Vertex
->NbColumns
- 2;
100 Param_Vertex_Common_Denominator(V
);
101 Vertex
= Matrix_Alloc(V
->Vertex
->NbRows
, V
->Vertex
->NbColumns
);
102 Matrix_Product(T
, V
->Vertex
, Vertex
);
103 for (i
= 0; i
< nvar
; ++i
) {
104 value_assign(Vertex
->p
[i
][nparam
+1], V
->Vertex
->p
[i
][nparam
+1]);
105 Vector_Normalize(Vertex
->p
[i
], nparam
+2);
107 Matrix_Free(V
->Vertex
);
111 /* Scales the parametric polyhedron with constraints *P and vertices PP
112 * such that the number of integer points can be represented by a polynomial.
113 * Both *P and P->Vertex are adapted according to the scaling.
114 * The scaling factor is returned in *det.
115 * The transformation that maps the new coordinates to the original
116 * coordinates is returned in *Lat (if Lat != NULL).
117 * The enumerator of the scaled parametric polyhedron should be divided
118 * by this number to obtain an approximation of the enumerator of the
119 * original parametric polyhedron.
121 * The algorithm is described in "Approximating Ehrhart Polynomials using
122 * affine transformations" by B. Meister.
124 void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron
*PP
, Polyhedron
**P
,
126 Value
*det
, unsigned MaxRays
)
129 unsigned dim
= (*P
)->Dimension
;
132 Lattice
*L
= NULL
, *Li
;
141 nparam
= PP
->V
->Vertex
->NbColumns
- 2;
144 for (V
= PP
->V
; V
; V
= V
->next
) {
148 unsigned char *supporting
;
150 supporting
= supporting_constraints(*P
, V
, &n
);
151 M
= Matrix_Alloc(n
, (*P
)->Dimension
);
152 for (i
= 0, j
= 0; i
< (*P
)->NbConstraints
; ++i
)
154 Vector_Copy((*P
)->Constraint
[i
]+1, M
->p
[j
++], (*P
)->Dimension
);
156 L2
= extract_lattice(M
, nvar
);
162 Lattice
*L3
= LatticeJoin(L
, L2
);
170 *Lat
= Matrix_Copy(L
);
172 /* apply the variable expansion to the polyhedron (constraints) */
173 expansion
= Matrix_Alloc(nvar
+ nparam
+ 1, nvar
+ nparam
+ 1);
174 for (i
= 0; i
< nvar
; ++i
)
175 Vector_Copy(L
->p
[i
], expansion
->p
[i
], nvar
);
176 for (i
= nvar
; i
< nvar
+nparam
+1; ++i
)
177 value_assign(expansion
->p
[i
][i
], L
->p
[nvar
][nvar
]);
179 *P
= Polyhedron_Preimage(*P
, expansion
, MaxRays
);
180 Matrix_Free(expansion
);
182 /* apply the variable expansion to the parametric vertices */
183 Li
= Matrix_Alloc(nvar
+1, nvar
+1);
184 ok
= Matrix_Inverse(L
, Li
);
187 assert(value_one_p(Li
->p
[nvar
][nvar
]));
188 T
= Matrix_Alloc(nvar
, nvar
);
189 value_set_si(*det
, 1);
190 for (i
= 0; i
< nvar
; ++i
) {
191 value_multiply(*det
, *det
, Li
->p
[i
][i
]);
192 Vector_Copy(Li
->p
[i
], T
->p
[i
], nvar
);
195 for (V
= PP
->V
; V
; V
= V
->next
)
196 Param_Vertex_Image(V
, T
);
200 /* Scales the parametric polyhedron with constraints *P and vertices PP
201 * such that the number of integer points can be represented by a polynomial.
202 * Both *P and P->Vertex are adapted according to the scaling.
203 * The scaling factor is returned in *det.
204 * The transformation that maps the new coordinates to the original
205 * coordinates is returned in *Lat (if Lat != NULL).
206 * The enumerator of the scaled parametric polyhedron should be divided
207 * by this number to obtain an approximation of the enumerator of the
208 * original parametric polyhedron.
210 * The algorithm is described in "Approximating Ehrhart Polynomials using
211 * affine transformations" by B. Meister.
213 void Param_Polyhedron_Scale_Integer_Fast(Param_Polyhedron
*PP
, Polyhedron
**P
,
215 Value
*det
, unsigned MaxRays
)
218 int nb_param
, nb_vars
;
221 Value global_var_lcm
;
225 value_set_si(*det
, 1);
229 nb_param
= PP
->D
->Domain
->Dimension
;
230 nb_vars
= PP
->V
->Vertex
->NbRows
;
232 /* Scan the vertices and make an orthogonal expansion of the variable
234 /* a- prepare the array of common denominators */
235 denoms
= Vector_Alloc(nb_vars
);
236 value_init(global_var_lcm
);
239 /* b- scan the vertices and compute the variables' global lcms */
240 for (V
= PP
->V
; V
; V
= V
->next
) {
241 for (i
= 0; i
< nb_vars
; i
++) {
242 Vector_Gcd(V
->Vertex
->p
[i
], nb_param
, &tmp
);
243 Gcd(tmp
, V
->Vertex
->p
[i
][nb_param
+1], &tmp
);
244 value_division(tmp
, V
->Vertex
->p
[i
][nb_param
+1], tmp
);
245 Lcm3(denoms
->p
[i
], tmp
, &denoms
->p
[i
]);
250 value_set_si(global_var_lcm
, 1);
251 for (i
= 0; i
< nb_vars
; i
++) {
252 value_multiply(*det
, *det
, denoms
->p
[i
]);
253 Lcm3(global_var_lcm
, denoms
->p
[i
], &global_var_lcm
);
257 for (V
= PP
->V
; V
; V
= V
->next
)
258 for (i
= 0; i
< nb_vars
; i
++) {
259 Vector_Scale(V
->Vertex
->p
[i
], V
->Vertex
->p
[i
], denoms
->p
[i
], nb_param
+1);
260 Vector_Normalize(V
->Vertex
->p
[i
], nb_param
+2);
263 /* the expansion can be actually writen as global_var_lcm.L^{-1} */
264 /* this is equivalent to multiply the rows of P by denoms_det */
265 for (i
= 0; i
< nb_vars
; i
++)
266 value_division(denoms
->p
[i
], global_var_lcm
, denoms
->p
[i
]);
268 /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
269 /* c- make the quick expansion matrix */
270 expansion
= Matrix_Alloc(nb_vars
+nb_param
+1, nb_vars
+nb_param
+1);
271 for (i
= 0; i
< nb_vars
; i
++)
272 value_assign(expansion
->p
[i
][i
], denoms
->p
[i
]);
273 for (i
= nb_vars
; i
< nb_vars
+nb_param
+1; i
++)
274 value_assign(expansion
->p
[i
][i
], global_var_lcm
);
276 /* d- apply the variable expansion to the polyhedron */
278 *P
= Polyhedron_Preimage(*P
, expansion
, MaxRays
);
281 Lattice
*L
= Matrix_Alloc(nb_vars
+1, nb_vars
+1);
282 for (i
= 0; i
< nb_vars
; ++i
)
283 value_assign(L
->p
[i
][i
], denoms
->p
[i
]);
284 value_assign(L
->p
[nb_vars
][nb_vars
], global_var_lcm
);
288 Matrix_Free(expansion
);
289 value_clear(global_var_lcm
);
293 /* Compute negated sum of all positive or negative coefficients in a row */
294 static void negated_sum(Value
*v
, int len
, int negative
, Value
*sum
)
298 value_set_si(*sum
, 0);
299 for (j
= 0; j
< len
; ++j
)
300 if (negative
? value_neg_p(v
[j
]) : value_pos_p(v
[j
]))
301 value_subtract(*sum
, *sum
, v
[j
]);
304 /* adapted from mpolyhedron_inflate in PolyLib */
305 Polyhedron
*Polyhedron_Flate(Polyhedron
*P
, unsigned nparam
, int inflate
,
309 int nvar
= P
->Dimension
- nparam
;
310 Matrix
*C
= Polyhedron2Constraints(P
);
315 /* subtract the sum of the negative coefficients of each inequality */
316 for (i
= 0; i
< C
->NbRows
; ++i
) {
317 negated_sum(C
->p
[i
]+1, nvar
, inflate
, &sum
);
318 value_addto(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
], sum
);
321 P2
= Constraints2Polyhedron(C
, MaxRays
);
326 C
= Polyhedron_Project(P
, nparam
);
327 CA
= align_context(C
, P
->Dimension
, MaxRays
);
329 P2
= DomainIntersection(P
, CA
, MaxRays
);
338 static Polyhedron
*flate_narrow2(Polyhedron
*P
, Lattice
*L
,
339 unsigned nparam
, int inflate
,
343 unsigned nvar
= P
->Dimension
- nparam
;
349 expansion
= Matrix_Alloc(nvar
+ nparam
+ 1, nvar
+ nparam
+ 1);
350 for (i
= 0; i
< nvar
; ++i
)
351 Vector_Copy(L
->p
[i
], expansion
->p
[i
], nvar
);
352 for (i
= nvar
; i
< nvar
+nparam
+1; ++i
)
353 value_assign(expansion
->p
[i
][i
], L
->p
[nvar
][nvar
]);
355 C
= Matrix_Alloc(P
->NbConstraints
+1, 1+P
->Dimension
+1);
357 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
358 negated_sum(P
->Constraint
[i
]+1, nvar
, inflate
, &sum
);
359 value_assign(C
->p
[i
][0], P
->Constraint
[i
][0]);
360 Vector_Matrix_Product(P
->Constraint
[i
]+1, expansion
, C
->p
[i
]+1);
361 if (value_zero_p(sum
))
363 Vector_Copy(C
->p
[i
]+1, C
->p
[i
+1]+1, P
->Dimension
+1);
364 value_addmul(C
->p
[i
][1+P
->Dimension
], sum
, L
->p
[nvar
][nvar
]);
365 ConstraintSimplify(C
->p
[i
], C
->p
[i
], P
->Dimension
+2, &sum
);
366 ConstraintSimplify(C
->p
[i
+1], C
->p
[i
+1], P
->Dimension
+2, &sum
);
367 if (value_ne(C
->p
[i
][1+P
->Dimension
], C
->p
[i
+1][1+P
->Dimension
])) {
369 value_decrement(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
]);
371 value_increment(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
]);
376 P2
= Constraints2Polyhedron(C
, MaxRays
);
379 Matrix_Free(expansion
);
383 C
= Polyhedron_Project(P
, nparam
);
384 CA
= align_context(C
, P
->Dimension
, MaxRays
);
386 P2
= DomainIntersection(P
, CA
, MaxRays
);
395 static void linear_min(Polyhedron
*D
, Value
*obj
, Value
*min
)
400 POL_ENSURE_VERTICES(D
);
401 for (i
= 0; i
< D
->NbRays
; ++i
) {
402 Inner_Product(obj
, D
->Ray
[i
]+1, D
->Dimension
, &tmp
);
403 mpz_cdiv_q(tmp
, tmp
, D
->Ray
[i
][1+D
->Dimension
]);
404 if (!i
|| value_lt(tmp
, *min
))
405 value_assign(*min
, tmp
);
410 static Polyhedron
*inflate_deflate_domain(Lattice
*L
, unsigned MaxRays
)
412 unsigned nvar
= L
->NbRows
-1;
417 M
= Matrix_Alloc(2*nvar
, 1+nvar
+1);
418 for (i
= 0; i
< nvar
; ++i
) {
419 value_set_si(M
->p
[2*i
][0], 1);
420 Vector_Copy(L
->p
[i
], M
->p
[2*i
]+1, nvar
);
421 Vector_Normalize(M
->p
[2*i
]+1, nvar
);
423 value_set_si(M
->p
[2*i
+1][0], 1);
424 Vector_Oppose(L
->p
[i
], M
->p
[2*i
+1]+1, nvar
);
425 value_assign(M
->p
[2*i
+1][1+nvar
], L
->p
[nvar
][nvar
]);
426 Vector_Normalize(M
->p
[2*i
+1]+1, nvar
+1);
427 value_decrement(M
->p
[2*i
+1][1+nvar
], M
->p
[2*i
+1][1+nvar
]);
429 D
= Constraints2Polyhedron(M
, MaxRays
);
435 static Polyhedron
*flate_narrow(Polyhedron
*P
, Lattice
*L
,
436 unsigned nparam
, int inflate
, unsigned MaxRays
)
439 unsigned nvar
= P
->Dimension
- nparam
;
446 D
= inflate_deflate_domain(L
, MaxRays
);
448 obj
= Vector_Alloc(nvar
);
449 C
= Polyhedron2Constraints(P
);
451 for (i
= 0; i
< C
->NbRows
; ++i
) {
453 Vector_Copy(C
->p
[i
]+1, obj
->p
, nvar
);
455 Vector_Oppose(C
->p
[i
]+1, obj
->p
, nvar
);
456 linear_min(D
, obj
->p
, &min
);
458 value_subtract(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
], min
);
460 value_addto(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
], min
);
464 P2
= Constraints2Polyhedron(C
, MaxRays
);
471 C
= Polyhedron_Project(P
, nparam
);
472 CA
= align_context(C
, P
->Dimension
, MaxRays
);
474 P2
= DomainIntersection(P
, CA
, MaxRays
);
483 static Polyhedron
*flate(Polyhedron
*P
, Lattice
*L
,
484 unsigned nparam
, int inflate
,
485 struct barvinok_options
*options
)
487 if (options
->scale_flags
& BV_APPROX_SCALE_NARROW2
)
488 return flate_narrow2(P
, L
, nparam
, inflate
, options
->MaxRays
);
489 else if (options
->scale_flags
& BV_APPROX_SCALE_NARROW
)
490 return flate_narrow(P
, L
, nparam
, inflate
, options
->MaxRays
);
492 return Polyhedron_Flate(P
, nparam
, inflate
, options
->MaxRays
);
495 static void Param_Polyhedron_Scale(Param_Polyhedron
*PP
, Polyhedron
**P
,
497 Value
*det
, struct barvinok_options
*options
)
499 if (options
->scale_flags
& BV_APPROX_SCALE_FAST
)
500 Param_Polyhedron_Scale_Integer_Fast(PP
, P
, L
, det
, options
->MaxRays
);
502 Param_Polyhedron_Scale_Integer_Slow(PP
, P
, L
, det
, options
->MaxRays
);
505 static evalue
*enumerate_flated(Polyhedron
*P
, Polyhedron
*C
, Lattice
*L
,
506 struct barvinok_options
*options
)
508 unsigned nparam
= C
->Dimension
;
510 int save_approximation
= options
->polynomial_approximation
;
512 if (options
->polynomial_approximation
== BV_APPROX_SIGN_UPPER
)
513 P
= flate(P
, L
, nparam
, 1, options
);
514 if (options
->polynomial_approximation
== BV_APPROX_SIGN_LOWER
)
515 P
= flate(P
, L
, nparam
, 0, options
);
517 /* Don't deflate/inflate again (on this polytope) */
518 options
->polynomial_approximation
= BV_APPROX_SIGN_NONE
;
519 eres
= barvinok_enumerate_with_options(P
, C
, options
);
520 options
->polynomial_approximation
= save_approximation
;
526 static evalue
*PP_enumerate_narrow_flated(Param_Polyhedron
*PP
,
527 Polyhedron
*P
, Polyhedron
*C
,
528 struct barvinok_options
*options
)
530 Polyhedron
*Porig
= P
;
531 int scale_narrow2
= options
->scale_flags
& BV_APPROX_SCALE_NARROW2
;
537 value_set_si(det
, 1);
539 Param_Polyhedron_Scale(PP
, &P
, &L
, &det
, options
);
540 Param_Polyhedron_Free(PP
);
545 /* Don't scale again (on this polytope) */
546 options
->approximation_method
= BV_APPROX_NONE
;
547 eres
= enumerate_flated(P
, C
, L
, options
);
548 options
->approximation_method
= BV_APPROX_SCALE
;
552 if (value_notone_p(det
))
553 evalue_div(eres
, det
);
558 static Param_Polyhedron
*Param_Polyhedron_Domain(Param_Polyhedron
*PP
,
563 Param_Polyhedron
*PP_D
;
566 Param_Vertices
**next
, *V
;
568 PP_D
= ALLOC(Param_Polyhedron
);
569 PP_D
->D
= ALLOC(Param_Domain
);
570 PP_D
->D
->next
= NULL
;
571 PP_D
->D
->Domain
= Domain_Copy(rVD
);
574 nv
= (PP
->nbV
- 1)/(8*sizeof(int)) + 1;
575 PP_D
->D
->F
= ALLOCN(unsigned, nv
);
576 memset(PP_D
->D
->F
, 0, nv
* sizeof(unsigned));
582 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
583 Param_Vertices
*V2
= ALLOC(Param_Vertices
);
584 V2
->Vertex
= Matrix_Copy(V
->Vertex
);
589 PP_D
->D
->F
[ix
] |= bx
;
592 END_FORALL_PVertex_in_ParamPolyhedron
;
598 static evalue
*enumerate_narrow_flated(Polyhedron
*P
, Polyhedron
*C
,
599 struct barvinok_options
*options
)
601 unsigned Rat_MaxRays
= options
->MaxRays
;
602 Param_Polyhedron
*PP
;
603 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
604 POL_UNSET(Rat_MaxRays
, POL_INTEGER
);
606 if ((options
->scale_flags
& BV_APPROX_SCALE_CHAMBER
) && PP
->D
->next
) {
608 evalue
*tmp
, *eres
= NULL
;
609 Polyhedron
*TC
= true_context(P
, C
, options
->MaxRays
);
611 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
613 Param_Polyhedron
*PP_D
;
614 /* Intersect with D->Domain, so we only have the relevant constraints
615 * left. Don't use rVD, though, since we still want to recognize
616 * the defining constraints of the parametric vertices.
618 CA
= align_context(D
->Domain
, P
->Dimension
, options
->MaxRays
);
619 P2
= DomainIntersection(P
, CA
, Rat_MaxRays
);
620 POL_ENSURE_VERTICES(P2
);
622 /* Use rVD for context, to avoid overlapping domains in
623 * results of computations in different chambers.
625 PP_D
= Param_Polyhedron_Domain(PP
, D
, rVD
);
626 tmp
= PP_enumerate_narrow_flated(PP_D
, P2
, rVD
, options
);
632 free_evalue_refs(tmp
);
635 Polyhedron_Free(rVD
);
636 END_FORALL_REDUCED_DOMAIN
637 Param_Polyhedron_Free(PP
);
639 eres
= evalue_zero();
643 return PP_enumerate_narrow_flated(PP
, P
, C
, options
);
646 /* If scaling is to be performed in combination with deflation/inflation,
647 * do both and return the result.
648 * Otherwise return NULL.
650 evalue
*scale_bound(Polyhedron
*P
, Polyhedron
*C
,
651 struct barvinok_options
*options
)
653 int scale_narrow
= options
->scale_flags
& BV_APPROX_SCALE_NARROW
;
654 int scale_narrow2
= options
->scale_flags
& BV_APPROX_SCALE_NARROW2
;
656 if (options
->polynomial_approximation
== BV_APPROX_SIGN_NONE
||
657 options
->polynomial_approximation
== BV_APPROX_SIGN_APPROX
)
660 if (scale_narrow
|| scale_narrow2
)
661 return enumerate_narrow_flated(P
, C
, options
);
663 return enumerate_flated(P
, C
, NULL
, options
);
666 evalue
*scale(Param_Polyhedron
*PP
, Polyhedron
*P
, Polyhedron
*C
,
667 struct barvinok_options
*options
)
674 if ((options
->scale_flags
& BV_APPROX_SCALE_CHAMBER
) && PP
->D
->next
) {
677 Polyhedron
*TC
= true_context(P
, C
, options
->MaxRays
);
679 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
680 Param_Polyhedron
*PP_D
= Param_Polyhedron_Domain(PP
, D
, rVD
);
681 tmp
= scale(PP_D
, P
, rVD
, options
);
686 free_evalue_refs(tmp
);
689 Param_Polyhedron_Free(PP_D
);
690 Polyhedron_Free(rVD
);
691 END_FORALL_REDUCED_DOMAIN
693 eres
= evalue_zero();
699 value_set_si(det
, 1);
701 MaxRays
= options
->MaxRays
;
702 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
703 Param_Polyhedron_Scale(PP
, &T
, NULL
, &det
, options
);
704 options
->MaxRays
= MaxRays
;
706 eres
= Param_Polyhedron_Enumerate(PP
, T
, C
, options
);
710 if (value_notone_p(det
))
711 evalue_div(eres
, det
);