From 01a8851f08646d301a4f9ab59c4cdef2dcd3cacc Mon Sep 17 00:00:00 2001
From: Sven Verdoolaege <skimo@kotnet.org>
Date: Mon, 19 Feb 2007 10:39:11 +0100
Subject: [PATCH] evalue_bernstein_coefficients: handle floorings in evalue

---
 barvinok/bernstein.h |   2 +
 bernstein.cc         | 116 ++++++++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 116 insertions(+), 2 deletions(-)

diff --git a/barvinok/bernstein.h b/barvinok/bernstein.h
index a780401..fb48e93 100644
--- a/barvinok/bernstein.h
+++ b/barvinok/bernstein.h
@@ -14,5 +14,7 @@ bernstein::piecewise_lst *evalue_bernstein_coefficients(
 	    Polyhedron *ctx, const GiNaC::exvector& params,
 	    barvinok_options *options);
 GiNaC::ex evalue2ex(evalue *e, const GiNaC::exvector& vars);
+GiNaC::ex evalue2ex(const evalue *e, const GiNaC::exvector& vars,
+		     GiNaC::exvector& floorvar, Matrix **C);
 
 }
diff --git a/bernstein.cc b/bernstein.cc
index 0f72699..c35e99c 100644
--- a/bernstein.cc
+++ b/bernstein.cc
@@ -1,3 +1,4 @@
+#include <vector>
 #include <bernstein/bernstein.h>
 #include <bernstein/piecewise_lst.h>
 #include <barvinok/barvinok.h>
@@ -8,6 +9,10 @@
 using namespace GiNaC;
 using namespace bernstein;
 
+using std::vector;
+using std::cerr;
+using std::endl;
+
 namespace barvinok {
 
 ex evalue2ex(evalue *e, const exvector& vars)
@@ -27,6 +32,101 @@ ex evalue2ex(evalue *e, const exvector& vars)
     return poly;
 }
 
+static int type_offset(enode *p)
+{
+   return p->type == fractional ? 1 : 
+	  p->type == flooring ? 1 : 0;
+}
+
+static ex evalue2ex_r(const evalue *e, const exvector& vars,
+		      exvector& floorvar, vector<const evalue *>& floors)
+{
+    if (value_notzero_p(e->d))
+	return value2numeric(e->x.n)/value2numeric(e->d);
+    ex base_var = 0;
+    ex poly = 0;
+
+    switch (e->x.p->type) {
+    case polynomial:
+	base_var = vars[e->x.p->pos-1];
+	break;
+    case flooring:
+	for (int i = 0; i < floors.size(); ++i) {
+	    if (eequal(&e->x.p->arr[0], floors[i])) {
+		base_var = floorvar[i];
+		break;
+	    }
+	}
+	if (base_var != 0)
+	    break;
+	char name[20];
+	snprintf(name, sizeof(name), "f%d", floors.size());
+	floorvar.push_back(base_var = symbol(name));
+	floors.push_back(&e->x.p->arr[0]);
+	break;
+    default:
+	return fail();
+    }
+
+    int offset = type_offset(e->x.p);
+    for (int i = e->x.p->size-1; i >= offset; --i) {
+	poly *= base_var;
+	ex t = evalue2ex_r(&e->x.p->arr[i], vars, floorvar, floors);
+	if (is_exactly_a<fail>(t))
+	    return t;
+	poly += t;
+    }
+    return poly;
+}
+
+/* For each t = floor(e/d), set up two constraints
+ *
+ *	 e - d t       >= 0
+ *	-e + d t + d-1 >= 0
+ *
+ * e is assumed to be an affine expression.
+ */
+static Matrix *setup_floor_constraints(const vector<const evalue*> floors, int nvar)
+{
+    int extra = floors.size();
+    if (!extra)
+	return NULL;
+    Matrix *M = Matrix_Alloc(2*extra, 1+extra+nvar+1);
+    for (int i = 0; i < extra; ++i) {
+	Value *d = &M->p[2*i][1+i];
+	value_set_si(*d, 1);
+	evalue_denom(floors[i], d);
+	const evalue *e;
+	for (e = floors[i]; value_zero_p(e->d); e = &e->x.p->arr[0]) {
+	    assert(e->x.p->type == polynomial);
+	    assert(e->x.p->size == 2);
+	    evalue *c = &e->x.p->arr[1];
+	    value_multiply(M->p[2*i][1+extra+e->x.p->pos-1], *d, c->x.n);
+	    value_division(M->p[2*i][1+extra+e->x.p->pos-1],
+			   M->p[2*i][1+extra+e->x.p->pos-1], c->d);
+	}
+	value_multiply(M->p[2*i][1+extra+nvar], *d, e->x.n);
+	value_division(M->p[2*i][1+extra+nvar], M->p[2*i][1+extra+nvar], e->d);
+	value_oppose(*d, *d);
+	value_set_si(M->p[2*i][0], -1);
+	Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
+	value_set_si(M->p[2*i][0], 1);
+	value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
+	value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
+    }
+    return M;
+}
+
+ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar, Matrix **C)
+{
+    vector<const evalue *> floors;
+    ex poly = evalue2ex_r(e, vars, floorvar, floors);
+    assert(C);
+    Matrix *M = setup_floor_constraints(floors, vars.size());
+    *C = M;
+    return poly;
+}
+
 /* if the evalue is a relation, we use the relation to cut off the 
  * the edges of the domain
  */
@@ -121,8 +221,20 @@ piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 	Param_Polyhedron *PP;
 	Polyhedron *E;
 	evalue *EP;
+	Matrix *M;
+	exvector floorvar;
+
 	evalue_extract_poly(e, i, &E, &EP, options->MaxRays);
-	ex poly = evalue2ex(EP, allvars);
+	ex poly = evalue2ex(EP, allvars, floorvar, &M);
+	floorvar.insert(floorvar.end(), vars.begin(), vars.end());
+	if (M) {
+	    Polyhedron *AE = align_context(E, M->NbColumns-2, options->MaxRays);
+	    if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
+		Domain_Free(E);
+	    E = DomainAddConstraints(AE, M, options->MaxRays);
+	    Matrix_Free(M);
+	    Domain_Free(AE);
+	}
 	if (is_exactly_a<fail>(poly)) {
 	    if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 		Domain_Free(E);
@@ -137,7 +249,7 @@ piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 	    PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 	    for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
 		matrix VM = domainVertices(PP, Q, params);
-		lst coeffs = bernsteinExpansion(VM, poly, vars, params);
+		lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
 		pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
 	    }
 	    Param_Polyhedron_Free(PP);
-- 
2.11.4.GIT