update isl for extract_interface
[barvinok.git] / scale.c
blobfcded17661e29f17cbd898fadb1337194c0acdbf
1 #include <assert.h>
2 #include <barvinok/barvinok.h>
3 #include <barvinok/util.h>
4 #include <barvinok/options.h>
5 #include "scale.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.
15 * We compute
16 * [ A B ] = [ H 0 ] [ U_11 U_12 ]
17 * [ U_21 U_22 ]
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)
26 int row;
27 Matrix *H, *Q, *U, *Li;
28 Lattice *L;
29 int ok;
31 left_hermite(M, &H, &Q, &U);
32 Matrix_Free(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);
40 Matrix_Free(H);
41 Matrix_Free(Q);
43 ok = Matrix_Inverse(Li, L);
44 assert(ok);
45 Matrix_Free(Li);
47 return L;
50 /* Returns the smallest (wrt inclusion) lattice that contains both X and Y */
51 static Lattice *LatticeJoin(Lattice *X, Lattice *Y)
53 int i;
54 int dim = X->NbRows-1;
55 Value lcm;
56 Value tmp;
57 Lattice *L;
58 Matrix *M, *H, *U, *Q;
60 assert(X->NbColumns-1 == dim);
61 assert(Y->NbRows-1 == dim);
62 assert(Y->NbColumns-1 == dim);
64 value_init(lcm);
65 value_init(tmp);
67 M = Matrix_Alloc(dim, 2*dim);
68 value_lcm(lcm, X->p[dim][dim], Y->p[dim][dim]);
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);
78 Matrix_Free(M);
79 Matrix_Free(Q);
80 Matrix_Free(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);
86 Matrix_Free(H);
88 value_clear(tmp);
89 value_clear(lcm);
90 return L;
93 static void Param_Vertex_Image(Param_Vertices *V, Matrix *T)
95 unsigned nvar = V->Vertex->NbRows;
96 unsigned nparam = V->Vertex->NbColumns - 2;
97 Matrix *Vertex;
98 int i;
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);
108 V->Vertex = Vertex;
111 static void apply_expansion(Param_Polyhedron *PP, Polyhedron **P,
112 Matrix *expansion, unsigned MaxRays)
114 int i;
115 unsigned nparam = PP->V->Vertex->NbColumns - 2;
116 unsigned nvar = PP->V->Vertex->NbRows;
117 Vector *constraint;
119 constraint = Vector_Alloc(nvar+nparam+1);
120 for (i = 0; i < PP->Constraints->NbRows; ++i) {
121 Vector_Matrix_Product(PP->Constraints->p[i]+1, expansion, constraint->p);
122 Vector_Copy(constraint->p, PP->Constraints->p[i]+1, nvar+nparam+1);
123 Vector_Normalize(PP->Constraints->p[i]+1, nvar+nparam+1);
125 Vector_Free(constraint);
126 if (P)
127 *P = Polyhedron_Preimage(*P, expansion, MaxRays);
130 /* Scales the parametric polyhedron with constraints *P and vertices PP
131 * such that the number of integer points can be represented by a polynomial.
132 * Both *P and P->Vertex are adapted according to the scaling.
133 * The scaling factor is returned in *det.
134 * The transformation that maps the new coordinates to the original
135 * coordinates is returned in *Lat (if Lat != NULL).
136 * The enumerator of the scaled parametric polyhedron should be divided
137 * by this number to obtain an approximation of the enumerator of the
138 * original parametric polyhedron.
140 * The algorithm is described in "Approximating Ehrhart Polynomials using
141 * affine transformations" by B. Meister.
143 void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron *PP, Polyhedron **P,
144 Lattice **Lat,
145 Value *det, unsigned MaxRays)
147 Param_Vertices *V;
148 unsigned dim = (*P)->Dimension;
149 unsigned nparam;
150 unsigned nvar;
151 Lattice *L = NULL, *Li;
152 Matrix *T;
153 Matrix *expansion;
154 int i;
155 int ok;
157 if (!PP->nbV)
158 return;
160 nparam = PP->V->Vertex->NbColumns - 2;
161 nvar = dim - nparam;
163 for (V = PP->V; V; V = V->next) {
164 Lattice *L2;
165 Matrix *M;
166 int i, j, n;
167 unsigned *supporting;
168 int ix;
169 unsigned bx;
171 supporting = supporting_constraints(PP->Constraints, V, &n);
172 M = Matrix_Alloc(n, PP->Constraints->NbColumns-2);
173 for (i = 0, j = 0, ix = 0, bx = MSB; i < PP->Constraints->NbRows; ++i) {
174 if (supporting[ix] & bx)
175 Vector_Copy(PP->Constraints->p[i]+1, M->p[j++],
176 PP->Constraints->NbColumns-2);
177 NEXT(ix, bx);
179 free(supporting);
180 L2 = extract_lattice(M, nvar);
181 Matrix_Free(M);
183 if (!L)
184 L = L2;
185 else {
186 Lattice *L3 = LatticeJoin(L, L2);
187 Matrix_Free(L);
188 Matrix_Free(L2);
189 L = L3;
193 if (Lat)
194 *Lat = Matrix_Copy(L);
196 /* apply the variable expansion to the polyhedron (constraints) */
197 expansion = Matrix_Alloc(nvar + nparam + 1, nvar + nparam + 1);
198 for (i = 0; i < nvar; ++i)
199 Vector_Copy(L->p[i], expansion->p[i], nvar);
200 for (i = nvar; i < nvar+nparam+1; ++i)
201 value_assign(expansion->p[i][i], L->p[nvar][nvar]);
203 apply_expansion(PP, P, expansion, MaxRays);
204 Matrix_Free(expansion);
206 /* apply the variable expansion to the parametric vertices */
207 Li = Matrix_Alloc(nvar+1, nvar+1);
208 ok = Matrix_Inverse(L, Li);
209 assert(ok);
210 Matrix_Free(L);
211 assert(value_one_p(Li->p[nvar][nvar]));
212 T = Matrix_Alloc(nvar, nvar);
213 value_set_si(*det, 1);
214 for (i = 0; i < nvar; ++i) {
215 value_multiply(*det, *det, Li->p[i][i]);
216 Vector_Copy(Li->p[i], T->p[i], nvar);
218 Matrix_Free(Li);
219 for (V = PP->V; V; V = V->next)
220 Param_Vertex_Image(V, T);
221 Matrix_Free(T);
224 /* Scales the parametric polyhedron with constraints *P and vertices PP
225 * such that the number of integer points can be represented by a polynomial.
226 * Both *P and P->Vertex are adapted according to the scaling.
227 * The scaling factor is returned in *det.
228 * The transformation that maps the new coordinates to the original
229 * coordinates is returned in *Lat (if Lat != NULL).
230 * The enumerator of the scaled parametric polyhedron should be divided
231 * by this number to obtain an approximation of the enumerator of the
232 * original parametric polyhedron.
234 * The algorithm is described in "Approximating Ehrhart Polynomials using
235 * affine transformations" by B. Meister.
237 void Param_Polyhedron_Scale_Integer_Fast(Param_Polyhedron *PP, Polyhedron **P,
238 Lattice **Lat,
239 Value *det, unsigned MaxRays)
241 int i;
242 int nb_param, nb_vars;
243 Vector *denoms;
244 Param_Vertices *V;
245 Value global_var_lcm;
246 Value tmp;
247 Matrix *expansion;
249 value_set_si(*det, 1);
250 if (!PP->nbV)
251 return;
253 nb_param = PP->D->Domain->Dimension;
254 nb_vars = PP->V->Vertex->NbRows;
256 /* Scan the vertices and make an orthogonal expansion of the variable
257 space */
258 /* a- prepare the array of common denominators */
259 denoms = Vector_Alloc(nb_vars);
260 value_init(global_var_lcm);
262 value_init(tmp);
263 /* b- scan the vertices and compute the variables' global lcms */
264 for (V = PP->V; V; V = V->next) {
265 for (i = 0; i < nb_vars; i++) {
266 Vector_Gcd(V->Vertex->p[i], nb_param, &tmp);
267 value_gcd(tmp, tmp, V->Vertex->p[i][nb_param+1]);
268 value_division(tmp, V->Vertex->p[i][nb_param+1], tmp);
269 Lcm3(denoms->p[i], tmp, &denoms->p[i]);
272 value_clear(tmp);
274 value_set_si(global_var_lcm, 1);
275 for (i = 0; i < nb_vars; i++) {
276 value_multiply(*det, *det, denoms->p[i]);
277 Lcm3(global_var_lcm, denoms->p[i], &global_var_lcm);
280 /* scale vertices */
281 for (V = PP->V; V; V = V->next)
282 for (i = 0; i < nb_vars; i++) {
283 Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
284 Vector_Normalize(V->Vertex->p[i], nb_param+2);
287 /* the expansion can be actually writen as global_var_lcm.L^{-1} */
288 /* this is equivalent to multiply the rows of P by denoms_det */
289 for (i = 0; i < nb_vars; i++)
290 value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
292 /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
293 /* c- make the quick expansion matrix */
294 expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
295 for (i = 0; i < nb_vars; i++)
296 value_assign(expansion->p[i][i], denoms->p[i]);
297 for (i = nb_vars; i < nb_vars+nb_param+1; i++)
298 value_assign(expansion->p[i][i], global_var_lcm);
300 /* d- apply the variable expansion to the polyhedron */
301 apply_expansion(PP, P, expansion, MaxRays);
303 if (Lat) {
304 Lattice *L = Matrix_Alloc(nb_vars+1, nb_vars+1);
305 for (i = 0; i < nb_vars; ++i)
306 value_assign(L->p[i][i], denoms->p[i]);
307 value_assign(L->p[nb_vars][nb_vars], global_var_lcm);
308 *Lat = L;
311 Matrix_Free(expansion);
312 value_clear(global_var_lcm);
313 Vector_Free(denoms);
316 /* Compute negated sum of all positive or negative coefficients in a row */
317 static void negated_sum(Value *v, int len, int negative, Value *sum)
319 int j;
321 value_set_si(*sum, 0);
322 for (j = 0; j < len; ++j)
323 if (negative ? value_neg_p(v[j]) : value_pos_p(v[j]))
324 value_subtract(*sum, *sum, v[j]);
327 /* adapted from mpolyhedron_inflate in PolyLib */
328 Polyhedron *Polyhedron_Flate(Polyhedron *P, unsigned nparam, int inflate,
329 unsigned MaxRays)
331 Value sum;
332 int nvar = P->Dimension - nparam;
333 Matrix *C = Polyhedron2Constraints(P);
334 Polyhedron *P2;
335 int i;
337 value_init(sum);
338 /* subtract the sum of the negative coefficients of each inequality */
339 for (i = 0; i < C->NbRows; ++i) {
340 negated_sum(C->p[i]+1, nvar, inflate, &sum);
341 value_addto(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], sum);
343 value_clear(sum);
344 P2 = Constraints2Polyhedron(C, MaxRays);
345 Matrix_Free(C);
347 if (inflate) {
348 Polyhedron *C, *CA;
349 C = Polyhedron_Project(P, nparam);
350 CA = align_context(C, P->Dimension, MaxRays);
351 P = P2;
352 P2 = DomainIntersection(P, CA, MaxRays);
353 Polyhedron_Free(C);
354 Polyhedron_Free(CA);
355 Polyhedron_Free(P);
358 return P2;
361 static Polyhedron *flate_narrow2(Polyhedron *P, Lattice *L,
362 unsigned nparam, int inflate,
363 unsigned MaxRays)
365 Value sum;
366 unsigned nvar = P->Dimension - nparam;
367 Matrix *expansion;
368 Matrix *C;
369 int i, j;
370 Polyhedron *P2;
372 expansion = Matrix_Alloc(nvar + nparam + 1, nvar + nparam + 1);
373 for (i = 0; i < nvar; ++i)
374 Vector_Copy(L->p[i], expansion->p[i], nvar);
375 for (i = nvar; i < nvar+nparam+1; ++i)
376 value_assign(expansion->p[i][i], L->p[nvar][nvar]);
378 C = Matrix_Alloc(P->NbConstraints+1, 1+P->Dimension+1);
379 value_init(sum);
380 for (i = 0; i < P->NbConstraints; ++i) {
381 negated_sum(P->Constraint[i]+1, nvar, inflate, &sum);
382 value_assign(C->p[i][0], P->Constraint[i][0]);
383 Vector_Matrix_Product(P->Constraint[i]+1, expansion, C->p[i]+1);
384 if (value_zero_p(sum))
385 continue;
386 Vector_Copy(C->p[i]+1, C->p[i+1]+1, P->Dimension+1);
387 value_addmul(C->p[i][1+P->Dimension], sum, L->p[nvar][nvar]);
388 ConstraintSimplify(C->p[i], C->p[i], P->Dimension+2, &sum);
389 ConstraintSimplify(C->p[i+1], C->p[i+1], P->Dimension+2, &sum);
390 if (value_ne(C->p[i][1+P->Dimension], C->p[i+1][1+P->Dimension])) {
391 if (inflate)
392 value_decrement(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension]);
393 else
394 value_increment(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension]);
397 value_clear(sum);
398 C->NbRows--;
399 P2 = Constraints2Polyhedron(C, MaxRays);
400 Matrix_Free(C);
402 Matrix_Free(expansion);
404 if (inflate) {
405 Polyhedron *C, *CA;
406 C = Polyhedron_Project(P, nparam);
407 CA = align_context(C, P->Dimension, MaxRays);
408 P = P2;
409 P2 = DomainIntersection(P, CA, MaxRays);
410 Polyhedron_Free(C);
411 Polyhedron_Free(CA);
412 Polyhedron_Free(P);
415 return P2;
418 static void linear_min(Polyhedron *D, Value *obj, Value *min)
420 int i;
421 Value tmp;
422 value_init(tmp);
423 POL_ENSURE_VERTICES(D);
424 for (i = 0; i < D->NbRays; ++i) {
425 Inner_Product(obj, D->Ray[i]+1, D->Dimension, &tmp);
426 mpz_cdiv_q(tmp, tmp, D->Ray[i][1+D->Dimension]);
427 if (!i || value_lt(tmp, *min))
428 value_assign(*min, tmp);
430 value_clear(tmp);
433 static Polyhedron *inflate_deflate_domain(Lattice *L, unsigned MaxRays)
435 unsigned nvar = L->NbRows-1;
436 int i;
437 Matrix *M;
438 Polyhedron *D;
440 M = Matrix_Alloc(2*nvar, 1+nvar+1);
441 for (i = 0; i < nvar; ++i) {
442 value_set_si(M->p[2*i][0], 1);
443 Vector_Copy(L->p[i], M->p[2*i]+1, nvar);
444 Vector_Normalize(M->p[2*i]+1, nvar);
446 value_set_si(M->p[2*i+1][0], 1);
447 Vector_Oppose(L->p[i], M->p[2*i+1]+1, nvar);
448 value_assign(M->p[2*i+1][1+nvar], L->p[nvar][nvar]);
449 Vector_Normalize(M->p[2*i+1]+1, nvar+1);
450 value_decrement(M->p[2*i+1][1+nvar], M->p[2*i+1][1+nvar]);
452 D = Constraints2Polyhedron(M, MaxRays);
453 Matrix_Free(M);
455 return D;
458 static Polyhedron *flate_narrow(Polyhedron *P, Lattice *L,
459 unsigned nparam, int inflate, unsigned MaxRays)
461 int i;
462 unsigned nvar = P->Dimension - nparam;
463 Vector *obj;
464 Value min;
465 Matrix *C;
466 Polyhedron *D;
467 Polyhedron *P2;
469 D = inflate_deflate_domain(L, MaxRays);
470 value_init(min);
471 obj = Vector_Alloc(nvar);
472 C = Polyhedron2Constraints(P);
474 for (i = 0; i < C->NbRows; ++i) {
475 if (inflate)
476 Vector_Copy(C->p[i]+1, obj->p, nvar);
477 else
478 Vector_Oppose(C->p[i]+1, obj->p, nvar);
479 linear_min(D, obj->p, &min);
480 if (inflate)
481 value_subtract(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], min);
482 else
483 value_addto(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], min);
486 Polyhedron_Free(D);
487 P2 = Constraints2Polyhedron(C, MaxRays);
488 Matrix_Free(C);
489 Vector_Free(obj);
490 value_clear(min);
492 if (inflate) {
493 Polyhedron *C, *CA;
494 C = Polyhedron_Project(P, nparam);
495 CA = align_context(C, P->Dimension, MaxRays);
496 P = P2;
497 P2 = DomainIntersection(P, CA, MaxRays);
498 Polyhedron_Free(C);
499 Polyhedron_Free(CA);
500 Polyhedron_Free(P);
503 return P2;
506 static Polyhedron *flate(Polyhedron *P, Lattice *L,
507 unsigned nparam, int inflate,
508 struct barvinok_options *options)
510 if (options->approx->scale_flags & BV_APPROX_SCALE_NARROW2)
511 return flate_narrow2(P, L, nparam, inflate, options->MaxRays);
512 else if (options->approx->scale_flags & BV_APPROX_SCALE_NARROW)
513 return flate_narrow(P, L, nparam, inflate, options->MaxRays);
514 else
515 return Polyhedron_Flate(P, nparam, inflate, options->MaxRays);
518 static void Param_Polyhedron_Scale(Param_Polyhedron *PP, Polyhedron **P,
519 Lattice **L,
520 Value *det, struct barvinok_options *options)
522 if (options->approx->scale_flags & BV_APPROX_SCALE_FAST)
523 Param_Polyhedron_Scale_Integer_Fast(PP, P, L, det, options->MaxRays);
524 else
525 Param_Polyhedron_Scale_Integer_Slow(PP, P, L, det, options->MaxRays);
528 static evalue *enumerate_flated(Polyhedron *P, Polyhedron *C, Lattice *L,
529 struct barvinok_options *options)
531 unsigned nparam = C->Dimension;
532 evalue *eres;
533 int save_approximation = options->approx->approximation;
535 if (options->approx->approximation == BV_APPROX_SIGN_UPPER)
536 P = flate(P, L, nparam, 1, options);
537 if (options->approx->approximation == BV_APPROX_SIGN_LOWER)
538 P = flate(P, L, nparam, 0, options);
540 /* Don't deflate/inflate again (on this polytope) */
541 options->approx->approximation = BV_APPROX_SIGN_NONE;
542 eres = barvinok_enumerate_with_options(P, C, options);
543 options->approx->approximation = save_approximation;
544 Polyhedron_Free(P);
546 return eres;
549 static evalue *PP_enumerate_narrow_flated(Param_Polyhedron *PP,
550 Polyhedron *P, Polyhedron *C,
551 struct barvinok_options *options)
553 Polyhedron *Porig = P;
554 int scale_narrow2 = options->approx->scale_flags & BV_APPROX_SCALE_NARROW2;
555 Lattice *L = NULL;
556 Value det;
557 evalue *eres;
559 value_init(det);
560 value_set_si(det, 1);
562 Param_Polyhedron_Scale(PP, &P, &L, &det, options);
563 Param_Polyhedron_Free(PP);
564 if (scale_narrow2) {
565 Polyhedron_Free(P);
566 P = Porig;
568 /* Don't scale again (on this polytope) */
569 options->approx->method = BV_APPROX_NONE;
570 eres = enumerate_flated(P, C, L, options);
571 options->approx->method = BV_APPROX_SCALE;
572 Matrix_Free(L);
573 if (P != Porig)
574 Polyhedron_Free(P);
575 if (value_notone_p(det))
576 evalue_div(eres, det);
577 value_clear(det);
578 return eres;
581 #define INT_BITS (sizeof(unsigned) * 8)
583 static Param_Polyhedron *Param_Polyhedron_Domain(Param_Polyhedron *PP,
584 Param_Domain *D,
585 Polyhedron *rVD)
587 int nv;
588 Param_Polyhedron *PP_D;
589 int i, ix;
590 unsigned bx;
591 Param_Vertices **next, *V;
592 int facet_len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
594 PP_D = ALLOC(Param_Polyhedron);
595 PP_D->D = ALLOC(Param_Domain);
596 PP_D->D->next = NULL;
597 PP_D->D->Domain = Domain_Copy(rVD);
598 PP_D->V = NULL;
599 PP_D->Constraints = Matrix_Copy(PP->Constraints);
600 PP_D->Rays = NULL;
602 nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
603 PP_D->D->F = ALLOCN(unsigned, nv);
604 memset(PP_D->D->F, 0, nv * sizeof(unsigned));
606 next = &PP_D->V;
607 i = 0;
608 ix = 0;
609 bx = MSB;
610 FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
611 Param_Vertices *V2 = ALLOC(Param_Vertices);
612 V2->Vertex = Matrix_Copy(V->Vertex);
613 V2->Domain = NULL;
614 V2->next = NULL;
615 V2->Facets = ALLOCN(unsigned, facet_len);
616 memcpy(V2->Facets, V->Facets, facet_len * sizeof(unsigned));
617 *next = V2;
618 next = &V2->next;
619 PP_D->D->F[ix] |= bx;
620 NEXT(ix, bx);
621 ++i;
622 END_FORALL_PVertex_in_ParamPolyhedron;
623 PP_D->nbV = i;
625 return PP_D;
628 static evalue *enumerate_narrow_flated(Polyhedron *P, Polyhedron *C,
629 struct barvinok_options *options)
631 unsigned Rat_MaxRays = options->MaxRays;
632 Param_Polyhedron *PP;
633 PP = Polyhedron2Param_Polyhedron(P, C, options);
634 POL_UNSET(Rat_MaxRays, POL_INTEGER);
636 if ((options->approx->scale_flags & BV_APPROX_SCALE_CHAMBER) && PP->D->next) {
637 int nd = -1;
638 evalue *tmp, *eres = NULL;
639 Polyhedron *TC = true_context(P, C, options->MaxRays);
641 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
642 Polyhedron *P2, *CA;
643 Param_Polyhedron *PP_D;
644 /* Intersect with D->Domain, so we only have the relevant constraints
645 * left. Don't use rVD, though, since we still want to recognize
646 * the defining constraints of the parametric vertices.
648 CA = align_context(D->Domain, P->Dimension, options->MaxRays);
649 P2 = DomainIntersection(P, CA, Rat_MaxRays);
650 POL_ENSURE_VERTICES(P2);
651 Polyhedron_Free(CA);
652 /* Use rVD for context, to avoid overlapping domains in
653 * results of computations in different chambers.
655 PP_D = Param_Polyhedron_Domain(PP, D, rVD);
656 tmp = PP_enumerate_narrow_flated(PP_D, P2, rVD, options);
657 Polyhedron_Free(P2);
658 if (!eres)
659 eres = tmp;
660 else {
661 eadd(tmp, eres);
662 free_evalue_refs(tmp);
663 free(tmp);
665 Polyhedron_Free(rVD);
666 END_FORALL_REDUCED_DOMAIN
667 Param_Polyhedron_Free(PP);
668 if (!eres)
669 eres = evalue_zero();
670 Polyhedron_Free(TC);
671 return eres;
672 } else
673 return PP_enumerate_narrow_flated(PP, P, C, options);
676 /* If scaling is to be performed in combination with deflation/inflation,
677 * do both and return the result.
678 * Otherwise return NULL.
680 evalue *scale_bound(Polyhedron *P, Polyhedron *C,
681 struct barvinok_options *options)
683 int scale_narrow = options->approx->scale_flags & BV_APPROX_SCALE_NARROW;
684 int scale_narrow2 = options->approx->scale_flags & BV_APPROX_SCALE_NARROW2;
686 if (options->approx->approximation == BV_APPROX_SIGN_NONE ||
687 options->approx->approximation == BV_APPROX_SIGN_APPROX)
688 return NULL;
690 if (scale_narrow || scale_narrow2)
691 return enumerate_narrow_flated(P, C, options);
692 else
693 return enumerate_flated(P, C, NULL, options);
696 evalue *scale(Param_Polyhedron *PP, Polyhedron *P, Polyhedron *C,
697 struct barvinok_options *options)
699 Polyhedron *T = P;
700 unsigned MaxRays;
701 evalue *eres = NULL;
702 Value det;
704 if ((options->approx->scale_flags & BV_APPROX_SCALE_CHAMBER) && PP->D->next) {
705 int nd = -1;
706 evalue *tmp;
707 Polyhedron *TC = true_context(P, C, options->MaxRays);
709 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
710 Param_Polyhedron *PP_D = Param_Polyhedron_Domain(PP, D, rVD);
711 tmp = scale(PP_D, P, rVD, options);
712 if (!eres)
713 eres = tmp;
714 else {
715 eadd(tmp, eres);
716 free_evalue_refs(tmp);
717 free(tmp);
719 Param_Polyhedron_Free(PP_D);
720 Polyhedron_Free(rVD);
721 END_FORALL_REDUCED_DOMAIN
722 if (!eres)
723 eres = evalue_zero();
724 Polyhedron_Free(TC);
725 return eres;
728 value_init(det);
729 value_set_si(det, 1);
731 MaxRays = options->MaxRays;
732 POL_UNSET(options->MaxRays, POL_INTEGER);
733 Param_Polyhedron_Scale(PP, &T, NULL, &det, options);
734 options->MaxRays = MaxRays;
736 eres = Param_Polyhedron_Enumerate(PP, T, C, options);
737 if (P != T)
738 Polyhedron_Free(T);
740 if (value_notone_p(det))
741 evalue_div(eres, det);
742 value_clear(det);
744 return eres;