3 #include <bernstein/bernstein.h>
4 #include <bernstein/piecewise_lst.h>
5 #include <barvinok/barvinok.h>
6 #include <barvinok/util.h>
7 #include <barvinok/bernstein.h>
8 #include <barvinok/options.h>
9 #include "reduce_domain.h"
11 using namespace GiNaC
;
12 using namespace bernstein
;
21 ex
evalue2ex(evalue
*e
, const exvector
& vars
)
23 if (value_pos_p(e
->d
))
24 return value2numeric(e
->x
.n
)/value2numeric(e
->d
);
25 if (EVALUE_IS_NAN(*e
))
27 if (e
->x
.p
->type
!= polynomial
)
30 for (int i
= e
->x
.p
->size
-1; i
>= 0; --i
) {
31 poly
*= vars
[e
->x
.p
->pos
-1];
32 ex t
= evalue2ex(&e
->x
.p
->arr
[i
], vars
);
33 if (is_exactly_a
<fail
>(t
))
40 static int type_offset(enode
*p
)
42 return p
->type
== fractional
? 1 :
43 p
->type
== flooring
? 1 : 0;
46 typedef pair
<bool, const evalue
*> typed_evalue
;
48 static ex
evalue2ex_add_var(evalue
*e
, exvector
& extravar
,
49 vector
<typed_evalue
>& expr
, bool is_fract
)
53 for (int i
= 0; i
< expr
.size(); ++i
) {
54 if (is_fract
== expr
[i
].first
&& eequal(e
, expr
[i
].second
)) {
55 base_var
= extravar
[i
];
63 snprintf(name
, sizeof(name
), "f%c%d", is_fract
? 'r' : 'l', expr
.size());
64 extravar
.push_back(base_var
= symbol(name
));
65 expr
.push_back(typed_evalue(is_fract
, e
));
70 /* For the argument e=(f/d) of a fractional, return (d-1)/d times
71 * a variable in [0,1] (see setup_constraints).
73 static ex
evalue2ex_get_fract(evalue
*e
, exvector
& extravar
,
74 vector
<typed_evalue
>& expr
)
82 den
= value2numeric(d
);
86 ex base_var
= evalue2ex_add_var(e
, extravar
, expr
, true);
91 static ex
evalue2ex_r(const evalue
*e
, const exvector
& vars
,
92 exvector
& extravar
, vector
<typed_evalue
>& expr
,
95 if (value_notzero_p(e
->d
))
96 return value2numeric(e
->x
.n
)/value2numeric(e
->d
);
100 switch (e
->x
.p
->type
) {
102 base_var
= vars
[e
->x
.p
->pos
-1];
105 base_var
= evalue2ex_add_var(&e
->x
.p
->arr
[0], extravar
, expr
, false);
108 base_var
= evalue2ex_get_fract(&e
->x
.p
->arr
[0], extravar
, expr
);
112 return evalue2ex_r(&e
->x
.p
->arr
[VALUE_TO_INT(coset
->p
[e
->x
.p
->pos
-1])],
113 vars
, extravar
, expr
, coset
);
118 int offset
= type_offset(e
->x
.p
);
119 for (int i
= e
->x
.p
->size
-1; i
>= offset
; --i
) {
121 ex t
= evalue2ex_r(&e
->x
.p
->arr
[i
], vars
, extravar
, expr
, coset
);
122 if (is_exactly_a
<fail
>(t
))
129 /* For each t = floor(e/d), set up two constraints
132 * -e + d t + d-1 >= 0
134 * e is assumed to be an affine expression.
136 * For each t = fract(e/d), set up two constraints
141 static Matrix
*setup_constraints(const vector
<typed_evalue
> expr
, int nvar
)
143 int extra
= expr
.size();
146 Matrix
*M
= Matrix_Alloc(2*extra
, 1+extra
+nvar
+1);
147 for (int i
= 0; i
< extra
; ++i
) {
149 value_set_si(M
->p
[2*i
][0], 1);
150 value_set_si(M
->p
[2*i
][1+i
], -1);
151 value_set_si(M
->p
[2*i
][1+extra
+nvar
], 1);
152 value_set_si(M
->p
[2*i
+1][0], 1);
153 value_set_si(M
->p
[2*i
+1][1+i
], 1);
155 Value
*d
= &M
->p
[2*i
][1+i
];
156 evalue_extract_affine(expr
[i
].second
, M
->p
[2*i
]+1+extra
,
157 M
->p
[2*i
]+1+extra
+nvar
, d
);
158 value_oppose(*d
, *d
);
159 value_set_si(M
->p
[2*i
][0], -1);
160 Vector_Scale(M
->p
[2*i
], M
->p
[2*i
+1], M
->p
[2*i
][0], 1+extra
+nvar
+1);
161 value_set_si(M
->p
[2*i
][0], 1);
162 value_subtract(M
->p
[2*i
+1][1+extra
+nvar
], M
->p
[2*i
+1][1+extra
+nvar
], *d
);
163 value_decrement(M
->p
[2*i
+1][1+extra
+nvar
], M
->p
[2*i
+1][1+extra
+nvar
]);
169 static bool evalue_is_periodic(const evalue
*e
, Vector
*periods
)
172 bool is_periodic
= false;
174 if (value_notzero_p(e
->d
))
177 assert(e
->x
.p
->type
!= partition
);
178 if (e
->x
.p
->type
== periodic
) {
181 value_set_si(size
, e
->x
.p
->size
);
182 value_lcm(periods
->p
[e
->x
.p
->pos
-1], periods
->p
[e
->x
.p
->pos
-1], size
);
186 offset
= type_offset(e
->x
.p
);
187 for (i
= e
->x
.p
->size
-1; i
>= offset
; --i
)
188 is_periodic
= evalue_is_periodic(&e
->x
.p
->arr
[i
], periods
) || is_periodic
;
192 static ex
evalue2lst(const evalue
*e
, const exvector
& vars
,
193 exvector
& extravar
, vector
<typed_evalue
>& expr
,
196 Vector
*coset
= Vector_Alloc(periods
->Size
);
200 list
.append(evalue2ex_r(e
, vars
, extravar
, expr
, coset
));
201 for (i
= coset
->Size
-1; i
>= 0; --i
) {
202 value_increment(coset
->p
[i
], coset
->p
[i
]);
203 if (value_lt(coset
->p
[i
], periods
->p
[i
]))
205 value_set_si(coset
->p
[i
], 0);
214 ex
evalue2ex(const evalue
*e
, const exvector
& vars
, exvector
& floorvar
,
215 Matrix
**C
, Vector
**p
)
217 vector
<typed_evalue
> expr
;
218 Vector
*periods
= Vector_Alloc(vars
.size());
221 for (int i
= 0; i
< periods
->Size
; ++i
)
222 value_set_si(periods
->p
[i
], 1);
223 if (evalue_is_periodic(e
, periods
)) {
229 Vector_Free(periods
);
231 ex poly
= evalue2ex_r(e
, vars
, floorvar
, expr
, NULL
);
232 Matrix
*M
= setup_constraints(expr
, vars
.size());
238 /* if the evalue is a relation, we use the relation to cut off the
239 * the edges of the domain
241 static Polyhedron
*relation_domain(Polyhedron
*D
, evalue
*fr
, unsigned MaxRays
)
243 assert(value_zero_p(fr
->d
));
244 assert(fr
->x
.p
->type
== fractional
);
245 assert(fr
->x
.p
->size
== 3);
246 Matrix
*T
= Matrix_Alloc(2, D
->Dimension
+1);
247 value_set_si(T
->p
[1][D
->Dimension
], 1);
249 /* convert argument of fractional to polylib */
250 /* the argument is assumed to be linear */
251 evalue
*p
= &fr
->x
.p
->arr
[0];
252 evalue_denom(p
, &T
->p
[1][D
->Dimension
]);
253 for (;value_zero_p(p
->d
); p
= &p
->x
.p
->arr
[0]) {
254 assert(p
->x
.p
->type
== polynomial
);
255 assert(p
->x
.p
->size
== 2);
256 assert(value_notzero_p(p
->x
.p
->arr
[1].d
));
257 int pos
= p
->x
.p
->pos
- 1;
258 value_assign(T
->p
[0][pos
], p
->x
.p
->arr
[1].x
.n
);
259 value_multiply(T
->p
[0][pos
], T
->p
[0][pos
], T
->p
[1][D
->Dimension
]);
260 value_division(T
->p
[0][pos
], T
->p
[0][pos
], p
->x
.p
->arr
[1].d
);
262 int pos
= D
->Dimension
;
263 value_assign(T
->p
[0][pos
], p
->x
.n
);
264 value_multiply(T
->p
[0][pos
], T
->p
[0][pos
], T
->p
[1][D
->Dimension
]);
265 value_division(T
->p
[0][pos
], T
->p
[0][pos
], p
->d
);
267 Polyhedron
*E
= NULL
;
268 for (Polyhedron
*P
= D
; P
; P
= P
->next
) {
269 Polyhedron
*I
= Polyhedron_Image(P
, T
, MaxRays
);
270 I
= DomainConstraintSimplify(I
, MaxRays
);
271 Polyhedron
*R
= Polyhedron_Preimage(I
, T
, MaxRays
);
273 Polyhedron
*next
= P
->next
;
275 Polyhedron
*S
= DomainIntersection(P
, R
, MaxRays
);
281 E
= DomainConcat(S
, E
);
288 piecewise_lst
*evalue_bernstein_coefficients(piecewise_lst
*pl_all
, evalue
*e
,
289 Polyhedron
*ctx
, const exvector
& params
)
292 barvinok_options
*options
= barvinok_options_new_with_defaults();
293 pl
= evalue_bernstein_coefficients(pl_all
, e
, ctx
, params
, options
);
294 barvinok_options_free(options
);
298 static piecewise_lst
*bernstein_coefficients(piecewise_lst
*pl_all
,
299 Polyhedron
*D
, const ex
& poly
,
301 const exvector
& params
, const exvector
& floorvar
,
302 barvinok_options
*options
);
304 /* Recursively apply Bernstein expansion on P, optimizing over dims[i]
305 * variables in each level. The context ctx is assumed to have been adapted
306 * to the first level in the recursion.
308 static piecewise_lst
*bernstein_coefficients_recursive(piecewise_lst
*pl_all
,
309 Polyhedron
*P
, const vector
<int>& dims
, const ex
& poly
,
311 const exvector
& params
, const exvector
& vars
,
312 barvinok_options
*options
)
314 assert(dims
.size() > 0);
315 assert(ctx
->Dimension
== P
->Dimension
- dims
[0]);
318 for (int j
= 0; j
< dims
.size(); ++j
) {
320 pl_vars
.insert(pl_vars
.end(), vars
.begin()+done
, vars
.begin()+done
+dims
[j
]);
322 pl_params
.insert(pl_params
.end(), vars
.begin()+done
+dims
[j
], vars
.end());
323 pl_params
.insert(pl_params
.end(), params
.begin(), params
.end());
326 pl
= bernstein_coefficients(NULL
, P
, poly
, ctx
,
327 pl_params
, pl_vars
, options
);
329 piecewise_lst
*new_pl
= NULL
;
330 Polyhedron
*U
= Universe_Polyhedron(pl_params
.size());
332 for (int i
= 0; i
< pl
->list
.size(); ++i
) {
333 Polyhedron
*D
= pl
->list
[i
].first
;
334 lst polys
= pl
->list
[i
].second
;
335 new_pl
= bernstein_coefficients(new_pl
, D
, polys
, U
, pl_params
,
354 pl_all
->combine(*pl
);
361 static piecewise_lst
*bernstein_coefficients_full_recurse(piecewise_lst
*pl_all
,
362 Polyhedron
*P
, const ex
& poly
,
364 const exvector
& params
, const exvector
& vars
,
365 barvinok_options
*options
)
367 Polyhedron
*CR
= align_context(ctx
, P
->Dimension
-1, options
->MaxRays
);
368 vector
<int> dims(vars
.size());
369 for (int i
= 0; i
< dims
.size(); ++i
)
371 pl_all
= bernstein_coefficients_recursive(pl_all
, P
, dims
, poly
, CR
,
372 params
, vars
, options
);
378 static piecewise_lst
*bernstein_coefficients_product(piecewise_lst
*pl_all
,
379 Polyhedron
*F
, Matrix
*T
, const ex
& poly
,
381 const exvector
& params
, const exvector
& vars
,
382 barvinok_options
*options
)
386 for (Polyhedron
*G
= F
; G
; G
= G
->next
)
390 unsigned nparam
= params
.size();
391 unsigned nvar
= vars
.size();
392 unsigned constraints
;
394 Polyhedron
*C
= NULL
;
396 /* More context constraints */
397 if (F
->Dimension
== ctx
->Dimension
) {
407 M
= Matrix_Alloc(F
->NbConstraints
, 1+nvar
+nparam
+1);
408 for (int i
= 0; i
< F
->NbConstraints
; ++i
) {
409 Vector_Copy(F
->Constraint
[i
], M
->p
[i
], 1+F
->Dimension
-nparam
);
410 Vector_Copy(F
->Constraint
[i
]+1+F
->Dimension
-nparam
,
411 M
->p
[i
]+1+nvar
, nparam
+1);
413 P
= Constraints2Polyhedron(M
, options
->MaxRays
);
417 constraints
= C
? C
->NbConstraints
: 0;
418 constraints
+= ctx
->NbConstraints
;
419 for (Polyhedron
*G
= F
->next
; G
; G
= G
->next
) {
420 constraints
+= G
->NbConstraints
;
424 unsigned total_var
= nvar
-(F
->Dimension
-nparam
);
427 M
= Matrix_Alloc(constraints
, 1+total_var
+nparam
+1);
428 for (Polyhedron
*G
= F
->next
; G
; G
= G
->next
) {
429 unsigned this_var
= G
->Dimension
- nparam
;
430 for (int i
= 0; i
< G
->NbConstraints
; ++i
) {
431 value_assign(M
->p
[c
+i
][0], G
->Constraint
[i
][0]);
432 Vector_Copy(G
->Constraint
[i
]+1, M
->p
[c
+i
]+1+skip
, this_var
);
433 Vector_Copy(G
->Constraint
[i
]+1+this_var
, M
->p
[c
+i
]+1+total_var
,
436 c
+= G
->NbConstraints
;
439 assert(skip
== total_var
);
441 for (int i
= 0; i
< C
->NbConstraints
; ++i
) {
442 value_assign(M
->p
[c
+i
][0], C
->Constraint
[i
][0]);
443 Vector_Copy(C
->Constraint
[i
]+1, M
->p
[c
+i
]+1+total_var
,
446 c
+= C
->NbConstraints
;
448 for (int i
= 0; i
< ctx
->NbConstraints
; ++i
) {
449 value_assign(M
->p
[c
+i
][0], ctx
->Constraint
[i
][0]);
450 Vector_Copy(ctx
->Constraint
[i
]+1, M
->p
[c
+i
]+1+total_var
, nparam
+1);
452 PC
= Constraints2Polyhedron(M
, options
->MaxRays
);
455 exvector newvars
= constructVariableVector(nvar
, "t");
456 matrix
subs(1, nvar
);
457 for (int i
= 0; i
< nvar
; ++i
)
458 for (int j
= 0; j
< nvar
; ++j
)
459 subs(0,i
) += value2numeric(T
->p
[i
][j
]) * newvars
[j
];
461 ex newpoly
= replaceVariablesInPolynomial(poly
, vars
, subs
);
463 vector
<int> dims(factors
);
464 for (int i
= 0; F
; ++i
, F
= F
->next
)
465 dims
[i
] = F
->Dimension
-nparam
;
467 pl_all
= bernstein_coefficients_recursive(pl_all
, P
, dims
, newpoly
, PC
,
468 params
, newvars
, options
);
476 static piecewise_lst
*bernstein_coefficients_polyhedron(piecewise_lst
*pl_all
,
477 Polyhedron
*P
, const ex
& poly
,
479 const exvector
& params
, const exvector
& floorvar
,
480 barvinok_options
*options
)
482 if (Polyhedron_is_unbounded(P
, ctx
->Dimension
, options
->MaxRays
)) {
483 fprintf(stderr
, "warning: unbounded domain skipped\n");
484 Polyhedron_Print(stderr
, P_VALUE_FMT
, P
);
488 if (options
->bernstein_recurse
& BV_BERNSTEIN_FACTORS
) {
490 Polyhedron
*F
= Polyhedron_Factor(P
, ctx
->Dimension
, &T
, options
->MaxRays
);
492 pl_all
= bernstein_coefficients_product(pl_all
, F
, T
, poly
, ctx
, params
,
499 if (floorvar
.size() > 1 &&
500 options
->bernstein_recurse
& BV_BERNSTEIN_INTERVALS
)
501 return bernstein_coefficients_full_recurse(pl_all
, P
, poly
, ctx
,
502 params
, floorvar
, options
);
504 unsigned PP_MaxRays
= options
->MaxRays
;
505 if (PP_MaxRays
& POL_NO_DUAL
)
508 Param_Polyhedron
*PP
= Polyhedron2Param_Domain(P
, ctx
, PP_MaxRays
);
511 piecewise_lst
*pl
= new piecewise_lst(params
, options
->bernstein_optimize
);
514 Polyhedron
*TC
= true_context(P
, ctx
, options
->MaxRays
);
515 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, PD
, rVD
)
516 matrix VM
= domainVertices(PP
, PD
, params
);
517 lst coeffs
= bernsteinExpansion(VM
, poly
, floorvar
, params
);
518 pl
->add_guarded_lst(rVD
, coeffs
);
519 END_FORALL_REDUCED_DOMAIN
522 Param_Polyhedron_Free(PP
);
526 pl_all
->combine(*pl
);
533 static piecewise_lst
*bernstein_coefficients(piecewise_lst
*pl_all
,
534 Polyhedron
*D
, const ex
& poly
,
536 const exvector
& params
, const exvector
& floorvar
,
537 barvinok_options
*options
)
539 if (!D
->next
&& emptyQ2(D
))
542 for (Polyhedron
*P
= D
; P
; P
= P
->next
) {
543 /* This shouldn't happen */
546 Polyhedron
*next
= P
->next
;
548 pl_all
= bernstein_coefficients_polyhedron(pl_all
, P
, poly
, ctx
,
549 params
, floorvar
, options
);
555 /* Compute the coefficients of the polynomial corresponding to each coset
556 * on its own domain. This allows us to cut the domain on multiples of
558 * To perform the cutting for a coset "i mod n = c" we map the domain
559 * to the quotient space trough "i = i' n + c", simplify the constraints
560 * (implicitly) and then map back to the original space.
562 static piecewise_lst
*bernstein_coefficients_periodic(piecewise_lst
*pl_all
,
563 Polyhedron
*D
, const evalue
*e
,
564 Polyhedron
*ctx
, const exvector
& vars
,
565 const exvector
& params
, Vector
*periods
,
566 barvinok_options
*options
)
568 assert(D
->Dimension
== periods
->Size
);
569 Matrix
*T
= Matrix_Alloc(D
->Dimension
+1, D
->Dimension
+1);
570 Matrix
*T2
= Matrix_Alloc(D
->Dimension
+1, D
->Dimension
+1);
571 Vector
*coset
= Vector_Alloc(periods
->Size
);
573 vector
<typed_evalue
> expr
;
574 exvector allvars
= vars
;
575 allvars
.insert(allvars
.end(), params
.begin(), params
.end());
577 value_set_si(T2
->p
[D
->Dimension
][D
->Dimension
], 1);
578 for (int i
= 0; i
< D
->Dimension
; ++i
) {
579 value_assign(T
->p
[i
][i
], periods
->p
[i
]);
580 value_lcm(T2
->p
[D
->Dimension
][D
->Dimension
],
581 T2
->p
[D
->Dimension
][D
->Dimension
], periods
->p
[i
]);
583 value_set_si(T
->p
[D
->Dimension
][D
->Dimension
], 1);
584 for (int i
= 0; i
< D
->Dimension
; ++i
)
585 value_division(T2
->p
[i
][i
], T2
->p
[D
->Dimension
][D
->Dimension
],
589 ex poly
= evalue2ex_r(e
, allvars
, extravar
, expr
, coset
);
590 assert(extravar
.size() == 0);
591 assert(expr
.size() == 0);
592 Polyhedron
*E
= DomainPreimage(D
, T
, options
->MaxRays
);
593 Polyhedron
*F
= DomainPreimage(E
, T2
, options
->MaxRays
);
596 pl_all
= bernstein_coefficients(pl_all
, F
, poly
, ctx
, params
,
599 for (i
= D
->Dimension
-1; i
>= 0; --i
) {
600 value_increment(coset
->p
[i
], coset
->p
[i
]);
601 value_increment(T
->p
[i
][D
->Dimension
], T
->p
[i
][D
->Dimension
]);
602 value_subtract(T2
->p
[i
][D
->Dimension
], T2
->p
[i
][D
->Dimension
],
604 if (value_lt(coset
->p
[i
], periods
->p
[i
]))
606 value_set_si(coset
->p
[i
], 0);
607 value_set_si(T
->p
[i
][D
->Dimension
], 0);
608 value_set_si(T2
->p
[i
][D
->Dimension
], 0);
619 piecewise_lst
*bernstein_coefficients_relation(piecewise_lst
*pl_all
,
620 Polyhedron
*D
, evalue
*EP
, Polyhedron
*ctx
,
621 const exvector
& allvars
, const exvector
& vars
,
622 const exvector
& params
, barvinok_options
*options
)
624 if (value_zero_p(EP
->d
) && EP
->x
.p
->type
== relation
) {
625 Polyhedron
*E
= relation_domain(D
, &EP
->x
.p
->arr
[0], options
->MaxRays
);
627 pl_all
= bernstein_coefficients_relation(pl_all
, E
, &EP
->x
.p
->arr
[1],
628 ctx
, allvars
, vars
, params
,
632 /* In principle, we could cut off the edges of this domain too */
633 if (EP
->x
.p
->size
> 2)
634 pl_all
= bernstein_coefficients_relation(pl_all
, D
, &EP
->x
.p
->arr
[2],
635 ctx
, allvars
, vars
, params
,
643 ex poly
= evalue2ex(EP
, allvars
, floorvar
, &M
, &periods
);
644 floorvar
.insert(floorvar
.end(), vars
.begin(), vars
.end());
647 Polyhedron
*AE
= align_context(D
, M
->NbColumns
-2, options
->MaxRays
);
648 E
= DomainAddConstraints(AE
, M
, options
->MaxRays
);
652 if (is_exactly_a
<fail
>(poly
)) {
657 pl_all
= bernstein_coefficients_periodic(pl_all
, E
, EP
, ctx
, vars
,
658 params
, periods
, options
);
660 pl_all
= bernstein_coefficients(pl_all
, E
, poly
, ctx
, params
,
663 Vector_Free(periods
);
670 piecewise_lst
*evalue_bernstein_coefficients(piecewise_lst
*pl_all
, evalue
*e
,
671 Polyhedron
*ctx
, const exvector
& params
,
672 barvinok_options
*options
)
674 unsigned nparam
= ctx
->Dimension
;
675 if (EVALUE_IS_ZERO(*e
))
677 assert(value_zero_p(e
->d
));
678 assert(e
->x
.p
->type
== partition
);
679 assert(e
->x
.p
->size
>= 2);
680 unsigned nvars
= EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
- nparam
;
682 exvector vars
= constructVariableVector(nvars
, "v");
683 exvector allvars
= vars
;
684 allvars
.insert(allvars
.end(), params
.begin(), params
.end());
686 for (int i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
687 pl_all
= bernstein_coefficients_relation(pl_all
,
688 EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]), &e
->x
.p
->arr
[2*i
+1],
689 ctx
, allvars
, vars
, params
, options
);