barvinok_maximize: add input conversion
[barvinok.git] / bernstein.cc
blob7c90c1a82a1204f977395e9ca055c5886aa24996
1 #include <vector>
2 #include <bernstein/bernstein.h>
3 #include <bernstein/piecewise_lst.h>
4 #include <barvinok/barvinok.h>
5 #include <barvinok/util.h>
6 #include <barvinok/bernstein.h>
7 #include <barvinok/options.h>
9 using namespace GiNaC;
10 using namespace bernstein;
12 using std::pair;
13 using std::vector;
14 using std::cerr;
15 using std::endl;
17 namespace barvinok {
19 ex evalue2ex(evalue *e, const exvector& vars)
21 if (value_notzero_p(e->d))
22 return value2numeric(e->x.n)/value2numeric(e->d);
23 if (e->x.p->type != polynomial)
24 return fail();
25 ex poly = 0;
26 for (int i = e->x.p->size-1; i >= 0; --i) {
27 poly *= vars[e->x.p->pos-1];
28 ex t = evalue2ex(&e->x.p->arr[i], vars);
29 if (is_exactly_a<fail>(t))
30 return t;
31 poly += t;
33 return poly;
36 static int type_offset(enode *p)
38 return p->type == fractional ? 1 :
39 p->type == flooring ? 1 : 0;
42 typedef pair<bool, const evalue *> typed_evalue;
44 static ex evalue2ex_add_var(evalue *e, exvector& extravar,
45 vector<typed_evalue>& expr, bool is_fract)
47 ex base_var = 0;
49 for (int i = 0; i < expr.size(); ++i) {
50 if (is_fract == expr[i].first && eequal(e, expr[i].second)) {
51 base_var = extravar[i];
52 break;
55 if (base_var != 0)
56 return base_var;
58 char name[20];
59 snprintf(name, sizeof(name), "f%c%d", is_fract ? 'r' : 'l', expr.size());
60 extravar.push_back(base_var = symbol(name));
61 expr.push_back(typed_evalue(is_fract, e));
63 return base_var;
66 static ex evalue2ex_r(const evalue *e, const exvector& vars,
67 exvector& extravar, vector<typed_evalue>& expr,
68 Vector *coset)
70 if (value_notzero_p(e->d))
71 return value2numeric(e->x.n)/value2numeric(e->d);
72 ex base_var = 0;
73 ex poly = 0;
75 switch (e->x.p->type) {
76 case polynomial:
77 base_var = vars[e->x.p->pos-1];
78 break;
79 case flooring:
80 base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, false);
81 break;
82 case fractional:
83 base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, true);
84 break;
85 case periodic:
86 assert(coset);
87 return evalue2ex_r(&e->x.p->arr[VALUE_TO_INT(coset->p[e->x.p->pos-1])],
88 vars, extravar, expr, coset);
89 default:
90 return fail();
93 int offset = type_offset(e->x.p);
94 for (int i = e->x.p->size-1; i >= offset; --i) {
95 poly *= base_var;
96 ex t = evalue2ex_r(&e->x.p->arr[i], vars, extravar, expr, coset);
97 if (is_exactly_a<fail>(t))
98 return t;
99 poly += t;
101 return poly;
104 /* For each t = floor(e/d), set up two constraints
106 * e - d t >= 0
107 * -e + d t + d-1 >= 0
109 * e is assumed to be an affine expression.
111 * For each t = fract(e/d), set up two constraints
113 * -d t + d-1 >= 0
114 * t >= 0
116 static Matrix *setup_constraints(const vector<typed_evalue> expr, int nvar)
118 int extra = expr.size();
119 if (!extra)
120 return NULL;
121 Matrix *M = Matrix_Alloc(2*extra, 1+extra+nvar+1);
122 for (int i = 0; i < extra; ++i) {
123 Value *d = &M->p[2*i][1+i];
124 value_set_si(*d, 1);
125 evalue_denom(expr[i].second, d);
126 if (expr[i].first) {
127 value_set_si(M->p[2*i][0], 1);
128 value_decrement(M->p[2*i][1+extra+nvar], *d);
129 value_oppose(*d, *d);
130 value_set_si(M->p[2*i+1][0], 1);
131 value_set_si(M->p[2*i+1][1+i], 1);
132 } else {
133 const evalue *e;
134 for (e = expr[i].second; value_zero_p(e->d); e = &e->x.p->arr[0]) {
135 assert(e->x.p->type == polynomial);
136 assert(e->x.p->size == 2);
137 evalue *c = &e->x.p->arr[1];
138 value_multiply(M->p[2*i][1+extra+e->x.p->pos-1], *d, c->x.n);
139 value_division(M->p[2*i][1+extra+e->x.p->pos-1],
140 M->p[2*i][1+extra+e->x.p->pos-1], c->d);
142 value_multiply(M->p[2*i][1+extra+nvar], *d, e->x.n);
143 value_division(M->p[2*i][1+extra+nvar], M->p[2*i][1+extra+nvar], e->d);
144 value_oppose(*d, *d);
145 value_set_si(M->p[2*i][0], -1);
146 Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
147 value_set_si(M->p[2*i][0], 1);
148 value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
149 value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
152 return M;
155 static bool evalue_is_periodic(const evalue *e, Vector *periods)
157 int i, offset;
158 bool is_periodic = false;
160 if (value_notzero_p(e->d))
161 return false;
163 assert(e->x.p->type != partition);
164 if (e->x.p->type == periodic) {
165 Value size;
166 value_init(size);
167 value_set_si(size, e->x.p->size);
168 value_lcm(periods->p[e->x.p->pos-1], size, &periods->p[e->x.p->pos-1]);
169 value_clear(size);
170 is_periodic = true;
172 offset = type_offset(e->x.p);
173 for (i = e->x.p->size-1; i >= offset; --i)
174 is_periodic = evalue_is_periodic(&e->x.p->arr[i], periods) || is_periodic;
175 return is_periodic;
178 static ex evalue2lst(const evalue *e, const exvector& vars,
179 exvector& extravar, vector<typed_evalue>& expr,
180 Vector *periods)
182 Vector *coset = Vector_Alloc(periods->Size);
183 lst list;
184 for (;;) {
185 int i;
186 list.append(evalue2ex_r(e, vars, extravar, expr, coset));
187 for (i = coset->Size-1; i >= 0; --i) {
188 value_increment(coset->p[i], coset->p[i]);
189 if (value_lt(coset->p[i], periods->p[i]))
190 break;
191 value_set_si(coset->p[i], 0);
193 if (i < 0)
194 break;
196 Vector_Free(coset);
197 return list;
200 ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar,
201 Matrix **C, Vector **p)
203 vector<typed_evalue> expr;
204 Vector *periods = Vector_Alloc(vars.size());
205 assert(p);
206 assert(C);
207 for (int i = 0; i < periods->Size; ++i)
208 value_set_si(periods->p[i], 1);
209 if (evalue_is_periodic(e, periods)) {
210 *p = periods;
211 *C = NULL;
212 lst list;
213 return list;
214 } else {
215 Vector_Free(periods);
216 *p = NULL;
217 ex poly = evalue2ex_r(e, vars, floorvar, expr, NULL);
218 Matrix *M = setup_constraints(expr, vars.size());
219 *C = M;
220 return poly;
224 /* if the evalue is a relation, we use the relation to cut off the
225 * the edges of the domain
227 static void evalue_extract_poly(evalue *e, int i, Polyhedron **D, evalue **poly,
228 unsigned MaxRays)
230 *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
231 *poly = e = &e->x.p->arr[2*i+1];
232 if (value_notzero_p(e->d))
233 return;
234 if (e->x.p->type != relation)
235 return;
236 if (e->x.p->size > 2)
237 return;
238 evalue *fr = &e->x.p->arr[0];
239 assert(value_zero_p(fr->d));
240 assert(fr->x.p->type == fractional);
241 assert(fr->x.p->size == 3);
242 Matrix *T = Matrix_Alloc(2, (*D)->Dimension+1);
243 value_set_si(T->p[1][(*D)->Dimension], 1);
245 /* convert argument of fractional to polylib */
246 /* the argument is assumed to be linear */
247 evalue *p = &fr->x.p->arr[0];
248 evalue_denom(p, &T->p[1][(*D)->Dimension]);
249 for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
250 assert(p->x.p->type == polynomial);
251 assert(p->x.p->size == 2);
252 assert(value_notzero_p(p->x.p->arr[1].d));
253 int pos = p->x.p->pos - 1;
254 value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
255 value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
256 value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
258 int pos = (*D)->Dimension;
259 value_assign(T->p[0][pos], p->x.n);
260 value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
261 value_division(T->p[0][pos], T->p[0][pos], p->d);
263 Polyhedron *E = NULL;
264 for (Polyhedron *P = *D; P; P = P->next) {
265 Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
266 I = DomainConstraintSimplify(I, MaxRays);
267 Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
268 Polyhedron_Free(I);
269 Polyhedron *next = P->next;
270 P->next = NULL;
271 Polyhedron *S = DomainIntersection(P, R, MaxRays);
272 Polyhedron_Free(R);
273 P->next = next;
274 if (emptyQ2(S))
275 Polyhedron_Free(S);
276 else
277 E = DomainConcat(S, E);
279 Matrix_Free(T);
281 *D = E;
282 *poly = &e->x.p->arr[1];
285 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
286 Polyhedron *ctx, const exvector& params)
288 piecewise_lst *pl;
289 barvinok_options *options = barvinok_options_new_with_defaults();
290 pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
291 barvinok_options_free(options);
292 return pl;
295 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
296 Polyhedron *D, const ex& poly,
297 Polyhedron *ctx,
298 const exvector& params, const exvector& floorvar,
299 barvinok_options *options)
301 if (emptyQ2(D))
302 return pl_all;
304 unsigned PP_MaxRays = options->MaxRays;
305 if (PP_MaxRays & POL_NO_DUAL)
306 PP_MaxRays = 0;
308 for (Polyhedron *P = D; P; P = P->next) {
309 Param_Polyhedron *PP;
310 Polyhedron *next = P->next;
311 Polyhedron *P1 = P;
312 P->next = NULL;
313 if (Polyhedron_is_unbounded(P, ctx->Dimension, options->MaxRays)) {
314 fprintf(stderr, "warning: unbounded domain skipped\n");
315 Polyhedron_Print(stderr, P_VALUE_FMT, P);
316 P->next = next;
317 continue;
319 PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
320 piecewise_lst *pl = new piecewise_lst(params);
321 for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
322 matrix VM = domainVertices(PP, Q, params);
323 lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
324 pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
326 Param_Polyhedron_Free(PP);
327 if (!pl_all)
328 pl_all = pl;
329 else {
330 pl_all->combine(*pl);
331 delete pl;
333 P->next = next;
335 return pl_all;
338 /* Compute the coefficients of the polynomial corresponding to each coset
339 * on its own domain. This allows us to cut the domain on multiples of
340 * the period.
341 * To perform the cutting for a coset "i mod n = c" we map the domain
342 * to the quotient space trough "i = i' n + c", simplify the constraints
343 * (implicitly) and then map back to the original space.
345 static piecewise_lst *bernstein_coefficients_periodic(piecewise_lst *pl_all,
346 Polyhedron *D, const evalue *e,
347 Polyhedron *ctx, const exvector& vars,
348 const exvector& params, Vector *periods,
349 barvinok_options *options)
351 assert(D->Dimension == periods->Size);
352 Matrix *T = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
353 Matrix *T2 = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
354 Vector *coset = Vector_Alloc(periods->Size);
355 exvector extravar;
356 vector<typed_evalue> expr;
357 exvector allvars = vars;
358 allvars.insert(allvars.end(), params.begin(), params.end());
360 value_set_si(T2->p[D->Dimension][D->Dimension], 1);
361 for (int i = 0; i < D->Dimension; ++i) {
362 value_assign(T->p[i][i], periods->p[i]);
363 value_lcm(T2->p[D->Dimension][D->Dimension], periods->p[i],
364 &T2->p[D->Dimension][D->Dimension]);
366 value_set_si(T->p[D->Dimension][D->Dimension], 1);
367 for (int i = 0; i < D->Dimension; ++i)
368 value_division(T2->p[i][i], T2->p[D->Dimension][D->Dimension],
369 periods->p[i]);
370 for (;;) {
371 int i;
372 ex poly = evalue2ex_r(e, allvars, extravar, expr, coset);
373 assert(extravar.size() == 0);
374 assert(expr.size() == 0);
375 Polyhedron *E = DomainPreimage(D, T, options->MaxRays);
376 Polyhedron *F = DomainPreimage(E, T2, options->MaxRays);
377 Polyhedron_Free(E);
378 if (!emptyQ2(F))
379 pl_all = bernstein_coefficients(pl_all, F, poly, ctx, params,
380 vars, options);
381 Polyhedron_Free(F);
382 for (i = D->Dimension-1; i >= 0; --i) {
383 value_increment(coset->p[i], coset->p[i]);
384 value_increment(T->p[i][D->Dimension], T->p[i][D->Dimension]);
385 value_subtract(T2->p[i][D->Dimension], T2->p[i][D->Dimension],
386 T2->p[i][i]);
387 if (value_lt(coset->p[i], periods->p[i]))
388 break;
389 value_set_si(coset->p[i], 0);
390 value_set_si(T->p[i][D->Dimension], 0);
391 value_set_si(T2->p[i][D->Dimension], 0);
393 if (i < 0)
394 break;
396 Vector_Free(coset);
397 Matrix_Free(T);
398 Matrix_Free(T2);
399 return pl_all;
402 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
403 Polyhedron *ctx, const exvector& params,
404 barvinok_options *options)
406 unsigned nparam = ctx->Dimension;
407 if (EVALUE_IS_ZERO(*e))
408 return pl_all;
409 assert(value_zero_p(e->d));
410 assert(e->x.p->type == partition);
411 assert(e->x.p->size >= 2);
412 unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
414 exvector vars = constructVariableVector(nvars, "v");
415 exvector allvars = vars;
416 allvars.insert(allvars.end(), params.begin(), params.end());
418 for (int i = 0; i < e->x.p->size/2; ++i) {
419 Polyhedron *E;
420 evalue *EP;
421 Matrix *M;
422 Vector *periods;
423 exvector floorvar;
425 evalue_extract_poly(e, i, &E, &EP, options->MaxRays);
426 ex poly = evalue2ex(EP, allvars, floorvar, &M, &periods);
427 floorvar.insert(floorvar.end(), vars.begin(), vars.end());
428 if (M) {
429 Polyhedron *AE = align_context(E, M->NbColumns-2, options->MaxRays);
430 if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
431 Domain_Free(E);
432 E = DomainAddConstraints(AE, M, options->MaxRays);
433 Matrix_Free(M);
434 Domain_Free(AE);
436 if (is_exactly_a<fail>(poly)) {
437 if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
438 Domain_Free(E);
439 delete pl_all;
440 return NULL;
442 if (periods)
443 pl_all = bernstein_coefficients_periodic(pl_all, E, EP, ctx, vars,
444 params, periods, options);
445 else
446 pl_all = bernstein_coefficients(pl_all, E, poly, ctx, params,
447 floorvar, options);
448 if (periods)
449 Vector_Free(periods);
450 if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
451 Domain_Free(E);
453 return pl_all;