1 #include <barvinok/util.h>
2 #include <barvinok/options.h>
5 /* If a vertex is described by A x + B p + c = 0, then
6 * M = [A B] and we want to compute a linear transformation L such
7 * that H L = A and H \Z contains both A \Z and B \Z.
9 * [ A B ] = [ H 0 ] [ U_11 U_12 ]
12 * U_11 is the required linear transformation.
13 * Note that this also works if M has more rows than there are variables,
14 * i.e., if some rows in M are linear combinations of other rows.
15 * These extra rows only affect and H and not U.
17 static Lattice
*extract_lattice(Matrix
*M
, unsigned nvar
)
20 Matrix
*H
, *Q
, *U
, *Li
;
24 left_hermite(M
, &H
, &Q
, &U
);
27 Li
= Matrix_Alloc(nvar
+1, nvar
+1);
28 L
= Matrix_Alloc(nvar
+1, nvar
+1);
29 value_set_si(Li
->p
[nvar
][nvar
], 1);
31 for (row
= 0; row
< nvar
; ++row
)
32 Vector_Copy(Q
->p
[row
], Li
->p
[row
], nvar
);
36 ok
= Matrix_Inverse(Li
, L
);
43 /* Returns the smallest (wrt inclusion) lattice that contains both X and Y */
44 static Lattice
*LatticeJoin(Lattice
*X
, Lattice
*Y
)
47 int dim
= X
->NbRows
-1;
51 Matrix
*M
, *H
, *U
, *Q
;
53 assert(X
->NbColumns
-1 == dim
);
54 assert(Y
->NbRows
-1 == dim
);
55 assert(Y
->NbColumns
-1 == dim
);
60 M
= Matrix_Alloc(dim
, 2*dim
);
61 value_lcm(X
->p
[dim
][dim
], Y
->p
[dim
][dim
], &lcm
);
63 value_division(tmp
, lcm
, X
->p
[dim
][dim
]);
64 for (i
= 0; i
< dim
; ++i
)
65 Vector_Scale(X
->p
[i
], M
->p
[i
], tmp
, dim
);
66 value_division(tmp
, lcm
, Y
->p
[dim
][dim
]);
67 for (i
= 0; i
< dim
; ++i
)
68 Vector_Scale(Y
->p
[i
], M
->p
[i
]+dim
, tmp
, dim
);
70 left_hermite(M
, &H
, &Q
, &U
);
75 L
= Matrix_Alloc(dim
+1, dim
+1);
76 value_assign(L
->p
[dim
][dim
], lcm
);
77 for (i
= 0; i
< dim
; ++i
)
78 Vector_Copy(H
->p
[i
], L
->p
[i
], dim
);
86 static void Param_Vertex_Common_Denominator(Param_Vertices
*V
)
92 assert(V
->Vertex
->NbRows
> 0);
93 dim
= V
->Vertex
->NbColumns
-2;
97 value_assign(lcm
, V
->Vertex
->p
[0][dim
+1]);
98 for (i
= 1; i
< V
->Vertex
->NbRows
; ++i
)
99 value_lcm(V
->Vertex
->p
[i
][dim
+1], lcm
, &lcm
);
101 for (i
= 0; i
< V
->Vertex
->NbRows
; ++i
) {
102 if (value_eq(V
->Vertex
->p
[i
][dim
+1], lcm
))
104 value_division(V
->Vertex
->p
[i
][dim
+1], lcm
, V
->Vertex
->p
[i
][dim
+1]);
105 Vector_Scale(V
->Vertex
->p
[i
], V
->Vertex
->p
[i
],
106 V
->Vertex
->p
[i
][dim
+1], dim
+1);
107 value_assign(V
->Vertex
->p
[i
][dim
+1], lcm
);
113 static void Param_Vertex_Image(Param_Vertices
*V
, Matrix
*T
)
115 unsigned nvar
= V
->Vertex
->NbRows
;
116 unsigned nparam
= V
->Vertex
->NbColumns
- 2;
120 Param_Vertex_Common_Denominator(V
);
121 Vertex
= Matrix_Alloc(V
->Vertex
->NbRows
, V
->Vertex
->NbColumns
);
122 Matrix_Product(T
, V
->Vertex
, Vertex
);
123 for (i
= 0; i
< nvar
; ++i
) {
124 value_assign(Vertex
->p
[i
][nparam
+1], V
->Vertex
->p
[i
][nparam
+1]);
125 Vector_Normalize(Vertex
->p
[i
], nparam
+2);
127 Matrix_Free(V
->Vertex
);
131 /* Scales the parametric polyhedron with constraints *P and vertices PP
132 * such that the number of integer points can be represented by a polynomial.
133 * Both *P and P->Vertex are adapted according to the scaling.
134 * The scaling factor is returned in *det.
135 * The transformation that maps the new coordinates to the original
136 * coordinates is returned in *Lat (if Lat != NULL).
137 * The enumerator of the scaled parametric polyhedron should be divided
138 * by this number to obtain an approximation of the enumerator of the
139 * original parametric polyhedron.
141 * The algorithm is described in "Approximating Ehrhart Polynomials using
142 * affine transformations" by B. Meister.
144 void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron
*PP
, Polyhedron
**P
,
146 Value
*det
, unsigned MaxRays
)
149 unsigned dim
= (*P
)->Dimension
;
152 Lattice
*L
= NULL
, *Li
;
161 nparam
= PP
->V
->Vertex
->NbColumns
- 2;
164 for (V
= PP
->V
; V
; V
= V
->next
) {
168 unsigned char *supporting
;
170 supporting
= supporting_constraints(*P
, V
, &n
);
171 M
= Matrix_Alloc(n
, (*P
)->Dimension
);
172 for (i
= 0, j
= 0; i
< (*P
)->NbConstraints
; ++i
)
174 Vector_Copy((*P
)->Constraint
[i
]+1, M
->p
[j
++], (*P
)->Dimension
);
176 L2
= extract_lattice(M
, nvar
);
182 Lattice
*L3
= LatticeJoin(L
, L2
);
190 *Lat
= Matrix_Copy(L
);
192 /* apply the variable expansion to the polyhedron (constraints) */
193 expansion
= Matrix_Alloc(nvar
+ nparam
+ 1, nvar
+ nparam
+ 1);
194 for (i
= 0; i
< nvar
; ++i
)
195 Vector_Copy(L
->p
[i
], expansion
->p
[i
], nvar
);
196 for (i
= nvar
; i
< nvar
+nparam
+1; ++i
)
197 value_assign(expansion
->p
[i
][i
], L
->p
[nvar
][nvar
]);
199 *P
= Polyhedron_Preimage(*P
, expansion
, MaxRays
);
200 Matrix_Free(expansion
);
202 /* apply the variable expansion to the parametric vertices */
203 Li
= Matrix_Alloc(nvar
+1, nvar
+1);
204 ok
= Matrix_Inverse(L
, Li
);
207 assert(value_one_p(Li
->p
[nvar
][nvar
]));
208 T
= Matrix_Alloc(nvar
, nvar
);
209 value_set_si(*det
, 1);
210 for (i
= 0; i
< nvar
; ++i
) {
211 value_multiply(*det
, *det
, Li
->p
[i
][i
]);
212 Vector_Copy(Li
->p
[i
], T
->p
[i
], nvar
);
215 for (V
= PP
->V
; V
; V
= V
->next
)
216 Param_Vertex_Image(V
, T
);
220 /* Scales the parametric polyhedron with constraints *P and vertices PP
221 * such that the number of integer points can be represented by a polynomial.
222 * Both *P and P->Vertex are adapted according to the scaling.
223 * The scaling factor is returned in *det.
224 * The transformation that maps the new coordinates to the original
225 * coordinates is returned in *Lat (if Lat != NULL).
226 * The enumerator of the scaled parametric polyhedron should be divided
227 * by this number to obtain an approximation of the enumerator of the
228 * original parametric polyhedron.
230 * The algorithm is described in "Approximating Ehrhart Polynomials using
231 * affine transformations" by B. Meister.
233 void Param_Polyhedron_Scale_Integer_Fast(Param_Polyhedron
*PP
, Polyhedron
**P
,
235 Value
*det
, unsigned MaxRays
)
238 int nb_param
, nb_vars
;
241 Value global_var_lcm
;
245 value_set_si(*det
, 1);
249 nb_param
= PP
->D
->Domain
->Dimension
;
250 nb_vars
= PP
->V
->Vertex
->NbRows
;
252 /* Scan the vertices and make an orthogonal expansion of the variable
254 /* a- prepare the array of common denominators */
255 denoms
= Vector_Alloc(nb_vars
);
256 value_init(global_var_lcm
);
259 /* b- scan the vertices and compute the variables' global lcms */
260 for (V
= PP
->V
; V
; V
= V
->next
) {
261 for (i
= 0; i
< nb_vars
; i
++) {
262 Vector_Gcd(V
->Vertex
->p
[i
], nb_param
, &tmp
);
263 Gcd(tmp
, V
->Vertex
->p
[i
][nb_param
+1], &tmp
);
264 value_division(tmp
, V
->Vertex
->p
[i
][nb_param
+1], tmp
);
265 Lcm3(denoms
->p
[i
], tmp
, &denoms
->p
[i
]);
270 value_set_si(global_var_lcm
, 1);
271 for (i
= 0; i
< nb_vars
; i
++) {
272 value_multiply(*det
, *det
, denoms
->p
[i
]);
273 Lcm3(global_var_lcm
, denoms
->p
[i
], &global_var_lcm
);
277 for (V
= PP
->V
; V
; V
= V
->next
)
278 for (i
= 0; i
< nb_vars
; i
++) {
279 Vector_Scale(V
->Vertex
->p
[i
], V
->Vertex
->p
[i
], denoms
->p
[i
], nb_param
+1);
280 Vector_Normalize(V
->Vertex
->p
[i
], nb_param
+2);
283 /* the expansion can be actually writen as global_var_lcm.L^{-1} */
284 /* this is equivalent to multiply the rows of P by denoms_det */
285 for (i
= 0; i
< nb_vars
; i
++)
286 value_division(denoms
->p
[i
], global_var_lcm
, denoms
->p
[i
]);
288 /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
289 /* c- make the quick expansion matrix */
290 expansion
= Matrix_Alloc(nb_vars
+nb_param
+1, nb_vars
+nb_param
+1);
291 for (i
= 0; i
< nb_vars
; i
++)
292 value_assign(expansion
->p
[i
][i
], denoms
->p
[i
]);
293 for (i
= nb_vars
; i
< nb_vars
+nb_param
+1; i
++)
294 value_assign(expansion
->p
[i
][i
], global_var_lcm
);
296 /* d- apply the variable expansion to the polyhedron */
298 *P
= Polyhedron_Preimage(*P
, expansion
, MaxRays
);
301 Lattice
*L
= Matrix_Alloc(nb_vars
+1, nb_vars
+1);
302 for (i
= 0; i
< nb_vars
; ++i
)
303 value_assign(L
->p
[i
][i
], denoms
->p
[i
]);
304 value_assign(L
->p
[nb_vars
][nb_vars
], global_var_lcm
);
308 Matrix_Free(expansion
);
309 value_clear(global_var_lcm
);
313 /* Compute negated sum of all positive or negative coefficients in a row */
314 static void negated_sum(Value
*v
, int len
, int negative
, Value
*sum
)
318 value_set_si(*sum
, 0);
319 for (j
= 0; j
< len
; ++j
)
320 if (negative
? value_neg_p(v
[j
]) : value_pos_p(v
[j
]))
321 value_subtract(*sum
, *sum
, v
[j
]);
324 /* adapted from mpolyhedron_inflate in PolyLib */
325 Polyhedron
*Polyhedron_Flate(Polyhedron
*P
, unsigned nparam
, int inflate
,
329 int nvar
= P
->Dimension
- nparam
;
330 Matrix
*C
= Polyhedron2Constraints(P
);
335 /* subtract the sum of the negative coefficients of each inequality */
336 for (i
= 0; i
< C
->NbRows
; ++i
) {
337 negated_sum(C
->p
[i
]+1, nvar
, inflate
, &sum
);
338 value_addto(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
], sum
);
341 P2
= Constraints2Polyhedron(C
, MaxRays
);
346 static Polyhedron
*flate_narrow2(Polyhedron
*P
, Lattice
*L
,
347 unsigned nparam
, int inflate
,
351 unsigned nvar
= P
->Dimension
- nparam
;
357 expansion
= Matrix_Alloc(nvar
+ nparam
+ 1, nvar
+ nparam
+ 1);
358 for (i
= 0; i
< nvar
; ++i
)
359 Vector_Copy(L
->p
[i
], expansion
->p
[i
], nvar
);
360 for (i
= nvar
; i
< nvar
+nparam
+1; ++i
)
361 value_assign(expansion
->p
[i
][i
], L
->p
[nvar
][nvar
]);
363 C
= Matrix_Alloc(P
->NbConstraints
+1, 1+P
->Dimension
+1);
365 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
366 negated_sum(P
->Constraint
[i
]+1, nvar
, inflate
, &sum
);
367 value_assign(C
->p
[i
][0], P
->Constraint
[i
][0]);
368 Vector_Matrix_Product(P
->Constraint
[i
]+1, expansion
, C
->p
[i
]+1);
369 if (value_zero_p(sum
))
371 Vector_Copy(C
->p
[i
]+1, C
->p
[i
+1]+1, P
->Dimension
+1);
372 value_addmul(C
->p
[i
][1+P
->Dimension
], sum
, L
->p
[nvar
][nvar
]);
373 ConstraintSimplify(C
->p
[i
], C
->p
[i
], P
->Dimension
+2, &sum
);
374 ConstraintSimplify(C
->p
[i
+1], C
->p
[i
+1], P
->Dimension
+2, &sum
);
375 if (value_ne(C
->p
[i
][1+P
->Dimension
], C
->p
[i
+1][1+P
->Dimension
])) {
377 value_decrement(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
]);
379 value_increment(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
]);
384 P2
= Constraints2Polyhedron(C
, MaxRays
);
387 Matrix_Free(expansion
);
391 static void linear_min(Polyhedron
*D
, Value
*obj
, Value
*min
)
396 POL_ENSURE_VERTICES(D
);
397 for (i
= 0; i
< D
->NbRays
; ++i
) {
398 Inner_Product(obj
, D
->Ray
[i
]+1, D
->Dimension
, &tmp
);
399 mpz_cdiv_q(tmp
, tmp
, D
->Ray
[i
][1+D
->Dimension
]);
400 if (!i
|| value_lt(tmp
, *min
))
401 value_assign(*min
, tmp
);
406 static void Vector_Oppose(Value
*p1
, Value
*p2
, unsigned len
)
410 for (i
= 0; i
< len
; ++i
)
411 value_oppose(p2
[i
], p1
[i
]);
414 static Polyhedron
*inflate_deflate_domain(Lattice
*L
, unsigned MaxRays
)
416 unsigned nvar
= L
->NbRows
-1;
421 M
= Matrix_Alloc(2*nvar
, 1+nvar
+1);
422 for (i
= 0; i
< nvar
; ++i
) {
423 value_set_si(M
->p
[2*i
][0], 1);
424 Vector_Copy(L
->p
[i
], M
->p
[2*i
]+1, nvar
);
425 Vector_Normalize(M
->p
[2*i
]+1, nvar
);
427 value_set_si(M
->p
[2*i
+1][0], 1);
428 Vector_Oppose(L
->p
[i
], M
->p
[2*i
+1]+1, nvar
);
429 value_assign(M
->p
[2*i
+1][1+nvar
], L
->p
[nvar
][nvar
]);
430 Vector_Normalize(M
->p
[2*i
+1]+1, nvar
+1);
431 value_decrement(M
->p
[2*i
+1][1+nvar
], M
->p
[2*i
+1][1+nvar
]);
433 D
= Constraints2Polyhedron(M
, MaxRays
);
439 static Polyhedron
*flate_narrow(Polyhedron
*P
, Lattice
*L
,
440 unsigned nparam
, int inflate
, unsigned MaxRays
)
443 unsigned nvar
= P
->Dimension
- nparam
;
450 D
= inflate_deflate_domain(L
, MaxRays
);
452 obj
= Vector_Alloc(nvar
);
453 C
= Polyhedron2Constraints(P
);
455 for (i
= 0; i
< C
->NbRows
; ++i
) {
457 Vector_Copy(C
->p
[i
]+1, obj
->p
, nvar
);
459 Vector_Oppose(C
->p
[i
]+1, obj
->p
, nvar
);
460 linear_min(D
, obj
->p
, &min
);
462 value_subtract(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
], min
);
464 value_addto(C
->p
[i
][1+P
->Dimension
], C
->p
[i
][1+P
->Dimension
], min
);
468 P2
= Constraints2Polyhedron(C
, MaxRays
);
475 static Polyhedron
*flate(Polyhedron
*P
, Lattice
*L
,
476 unsigned nparam
, int inflate
,
477 struct barvinok_options
*options
)
479 if (options
->scale_flags
& BV_APPROX_SCALE_NARROW2
)
480 return flate_narrow2(P
, L
, nparam
, inflate
, options
->MaxRays
);
481 else if (options
->scale_flags
& BV_APPROX_SCALE_NARROW
)
482 return flate_narrow(P
, L
, nparam
, inflate
, options
->MaxRays
);
484 return Polyhedron_Flate(P
, nparam
, inflate
, options
->MaxRays
);
487 static void Param_Polyhedron_Scale(Param_Polyhedron
*PP
, Polyhedron
**P
,
489 Value
*det
, struct barvinok_options
*options
)
491 if (options
->scale_flags
& BV_APPROX_SCALE_FAST
)
492 Param_Polyhedron_Scale_Integer_Fast(PP
, P
, L
, det
, options
->MaxRays
);
494 Param_Polyhedron_Scale_Integer_Slow(PP
, P
, L
, det
, options
->MaxRays
);
497 Polyhedron
*scale_init(Polyhedron
*P
, Polyhedron
*C
, struct scale_data
*scaling
,
498 struct barvinok_options
*options
)
500 unsigned nparam
= C
->Dimension
;
501 Polyhedron
*Porig
= P
;
503 int scale_narrow
= options
->scale_flags
& BV_APPROX_SCALE_NARROW
;
504 int scale_narrow2
= options
->scale_flags
& BV_APPROX_SCALE_NARROW2
;
507 value_init(scaling
->det
);
508 value_set_si(scaling
->det
, 1);
509 scaling
->save_approximation
= options
->polynomial_approximation
;
511 if (options
->polynomial_approximation
== BV_APPROX_SIGN_NONE
||
512 options
->polynomial_approximation
== BV_APPROX_SIGN_APPROX
)
515 if (scale_narrow
|| scale_narrow2
) {
516 Param_Polyhedron
*PP
;
517 unsigned PP_MaxRays
= options
->MaxRays
;
518 if (PP_MaxRays
& POL_NO_DUAL
)
520 PP
= Polyhedron2Param_Domain(P
, C
, PP_MaxRays
);
521 Param_Polyhedron_Scale(PP
, &P
, &L
, &scaling
->det
, options
);
522 Param_Polyhedron_Free(PP
);
527 /* Don't scale again (on this polytope) */
528 options
->approximation_method
= BV_APPROX_NONE
;
531 if (options
->polynomial_approximation
== BV_APPROX_SIGN_UPPER
)
532 P
= flate(P
, L
, nparam
, 1, options
);
533 if (options
->polynomial_approximation
== BV_APPROX_SIGN_LOWER
)
534 P
= flate(P
, L
, nparam
, 0, options
);
536 /* Don't deflate/inflate again (on this polytope) */
537 options
->polynomial_approximation
= BV_APPROX_SIGN_NONE
;
546 Polyhedron
*scale(Param_Polyhedron
*PP
, Polyhedron
*P
,
547 struct scale_data
*scaling
, int free_P
,
548 struct barvinok_options
*options
)
553 MaxRays
= options
->MaxRays
;
554 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
556 Param_Polyhedron_Scale(PP
, &T
, NULL
, &scaling
->det
, options
);
560 options
->MaxRays
= MaxRays
;
565 void scale_finish(evalue
*e
, struct scale_data
*scaling
,
566 struct barvinok_options
*options
)
568 if (value_notone_p(scaling
->det
))
569 evalue_div(e
, scaling
->det
);
570 value_clear(scaling
->det
);
571 /* reset options that may have been changed */
572 options
->approximation_method
= BV_APPROX_SCALE
;
573 options
->polynomial_approximation
= scaling
->save_approximation
;