From adf5021392c6f641217bc3d6e4b3e36dc77b5e47 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Tue, 26 Feb 2008 14:44:28 +0100 Subject: [PATCH] Bernoulli_sum_evalue: cut off some redundant parts where sum should be zero --- bernoulli.c | 11 ++++++++--- doc/implementation.tex | 38 +++++++++++++++++++++++++++++++++++--- 2 files changed, 43 insertions(+), 6 deletions(-) diff --git a/bernoulli.c b/bernoulli.c index 1293a2f..f86ceeb 100644 --- a/bernoulli.c +++ b/bernoulli.c @@ -336,9 +336,10 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) * upper is the constraint * -a i + u' >= 0 i <= u'/a = u */ - M2 = Matrix_Alloc(2, 2+T->Dimension); + M2 = Matrix_Alloc(3, 2+T->Dimension); value_set_si(M2->p[0][0], 1); value_set_si(M2->p[1][0], 1); + value_set_si(M2->p[2][0], 1); /* Case 1: * 1 <= l l' - b >= 0 */ @@ -410,10 +411,12 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) /* Case 4: * l < 1 -l' + b - 1 >= 0 * 0 < l l' - 1 >= 0 + * u >= 1 u' - a >= 0 */ bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1); bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1); - D = AddConstraints(M2->p_Init, 2, T, data->MaxRays); + bound_constraint(M2->p[2]+1, T->Dimension, upper+1, 1, 1, 0); + D = AddConstraints(M2->p_Init, 3, T, data->MaxRays); POL_ENSURE_VERTICES(D); if (emptyQ2(D)) Polyhedron_Free(D); @@ -431,10 +434,12 @@ static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data) /* Case 5: * -u < 1 u' + a - 1 >= 0 * 0 < -u -u' - 1 >= 0 + * l <= -1 -l' - b >= 0 */ bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1); bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1); - D = AddConstraints(M2->p_Init, 2, T, data->MaxRays); + bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, -1, 0); + D = AddConstraints(M2->p_Init, 3, T, data->MaxRays); POL_ENSURE_VERTICES(D); if (emptyQ2(D)) Polyhedron_Free(D); diff --git a/doc/implementation.tex b/doc/implementation.tex index 2809291..68628c9 100644 --- a/doc/implementation.tex +++ b/doc/implementation.tex @@ -1977,6 +1977,7 @@ cases. \end{align*} \item $l(\vec {\hat x}) \le 0$ and $u(\vec{\hat x}) \ge 0$ +\label{i:l:u} \begin{align*} \sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n @@ -2007,14 +2008,13 @@ integer, then the above 3 cases partition (the integer points in) the projection of $P$ onto the remaining variables. However, if some of the coefficients are rational, then the lower and upper bound can lie in the open interval $(0,1)$ for some -values of $\vec{\hat x}$. We therefore also need to consider +values of $\vec{\hat x}$. We may therefore also want to consider the following two cases. -Note that the results will not be correct in these cases, but -not taking them into account would lead to a greater loss in accuracy. \begin{enumerate} \setcounter{enumi}{\value{saveenumi}} \item $0 < l(\vec {\hat x}) < 1$ +\label{i:l:fractional} $$ \sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n @@ -2024,6 +2024,7 @@ S_n(u(\vec{\hat x})+1) $$ \item $0 < -u(\vec {\hat x}) < 1$ +\label{i:u:fractional} $$ \sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n @@ -2034,6 +2035,37 @@ S_n(-l(\vec{\hat x})+1) $$ \end{enumerate} +Note that we may add the constraint $u \ge 1$ to +case~\ref{i:l:fractional} and the constraint $l \le -1$ +to case~\ref{i:u:fractional}, since the correct value for +these two cases would be zero if these extra constraints do not hold. + +An alternative to adding the above two cases would be +to simply ignore them, i.e., assume a value of $0$. +Another alternative would be to reduce case~\ref{i:l:u} +to +$$ +l(\vec {\hat x}) \le -1\quad\hbox{and}\quad u(\vec{\hat x}) \ge 1 +, +$$ +while extending cases~\ref{i:l:fractional} and~\ref{i:u:fractional} +to +$$ +-1 < l(\vec {\hat x}) < 1\quad\hbox{and}\quad u \ge 1 +$$ +and +$$ +-1 < u(\vec {\hat x}) < 1\quad\hbox{and}\quad l \le -1 +, +$$ +respectively, with the remaining cases +($-1 < l \le u < 1$) having value $0$. +There does not appear to be a consistently better choice +here, as each of these three approaches seems to yield better +results on some examples. +The last approach has the additional drawback that we +would also have to deal with 5 cases, even if the bounds +are integers. \subsection{Summation using local Euler-Maclaurin formula} \label{s:euler} -- 2.11.4.GIT