new constraints poly-eq-k, poly-ne, and poly-ne-k

Vasco Pedro
1 parent acfea957
Showing 6 changed files with 1022 additions and 0 deletions Show diff stats
README
src/constraints.c
src/constraints.h
src/constraints/poly-eq-k.c
src/constraints/poly-ne-k.c
src/constraints/poly-ne.c
@@ -173,6 +173,9 @@ Global constraints
   fd_sum_prod(fd_int X[], fd_int Y[], int n, int k)	X . Y = k
  
   fd_poly_eq(int C[], fd_int Y[], int n, fd_int z)	z <= C . Y <= z
+  fd_poly_eq_k(int C[], fd_int Y[], int n, int k)	C . Y = k
+  fd_poly_ne(int C[], fd_int Y[], int n, fd_int z)	C . Y != z
+  fd_poly_ne_k(int C[], fd_int Y[], int n, int k)	C . Y != k
   fd_knapsack2(fd_int X[], fd_int Y[], int n, fd_int z)	z <= X . Y <= z
  
   fd_exactly(fd_int X[], int n, int c, int k)		 #{ i | X[i] = k } = c
@@ -149,6 +149,9 @@ void fd__constraint_free_memory()
 #include <constraints/poly-eq.c>
 #include <constraints/element-var.c>
 #include <constraints/exactly-vars.c>
+#include <constraints/poly-eq-k.c>
+#include <constraints/poly-ne.c>
+#include <constraints/poly-ne-k.c>
  
 #ifdef CONSTRAINT_CLASS
  
@@ -178,6 +181,9 @@ _FD_CONSTRAINT_INITIALISATION(max, FD_CONSTR_MAX)
 _FD_CONSTRAINT_INITIALISATION(poly_eq, FD_CONSTR_POLY_EQ)
 _FD_CONSTRAINT_INITIALISATION(element_var, FD_CONSTR_ELEMENT_VAR)
 _FD_CONSTRAINT_INITIALISATION(exactly_vars, FD_CONSTR_EXACTLY_VARS)
+_FD_CONSTRAINT_INITIALISATION(poly_eq_k, FD_CONSTR_POLY_EQ_K)
+_FD_CONSTRAINT_INITIALISATION(poly_ne, FD_CONSTR_POLY_NE)
+_FD_CONSTRAINT_INITIALISATION(poly_ne_k, FD_CONSTR_POLY_NE_K)
  
 void _fd_init_constraints()
 {
@@ -212,6 +218,9 @@ void _fd_init_constraints()
   _fd_init_constraint(poly_eq);
   _fd_init_constraint(element_var);
   _fd_init_constraint(exactly_vars);
+  _fd_init_constraint(poly_eq_k);
+  _fd_init_constraint(poly_ne);
+  _fd_init_constraint(poly_ne_k);
  
   done = 1;
 }
@@ -33,6 +33,9 @@ typedef enum {
   FD_CONSTR_POLY_EQ,
   FD_CONSTR_ELEMENT_VAR,
   FD_CONSTR_EXACTLY_VARS,
+  FD_CONSTR_POLY_EQ_K,
+  FD_CONSTR_POLY_NE,
+  FD_CONSTR_POLY_NE_K,
   FD_CONSTR_KINDS		// number of kinds of constraints
 } _fd_constraint_kind;
  
@@ -0,0 +1,469 @@
+/* poly-eq-k(C, X, k) == sum(C . X) = k */
+
+static int fd_poly_eq_k_filter(fd_constraint this)
+{
+  int k;
+  int min, max;
+  int terms = this->nvariables;
+  int *mins, *maxs;
+  int i;
+#ifdef CONSTRAINT_TEMPS
+  int *base;
+
+  assert(!fd__constraint_data_valid(this));
+
+  if (!constraint_memory[this->index])
+    constraint_memory[this->index] = malloc((2 * terms + 2) * sizeof(int));
+
+  base = constraint_memory[this->index];
+
+  mins = base + 2;
+  maxs = mins + terms;
+#else
+  mins = alloca(terms * sizeof(*mins));
+  maxs = alloca(terms * sizeof(*maxs));
+#endif
+
+  k = this->constants[terms];
+
+  // sum the minima and the maxima of the terms
+  min = max = 0;
+
+  for (i = 0; i < terms; ++i)
+    {
+      int vl, vh;
+
+      if (this->constants[i] > 0)
+	{
+	  vl = mins[i] = _fd_var_min(VAR(this, i));
+	  vh = maxs[i] = _fd_var_max(VAR(this, i));
+	}
+      else
+	{
+	  vl = maxs[i] = _fd_var_max(VAR(this, i));
+	  vh = mins[i] = _fd_var_min(VAR(this, i));
+	}
+
+      min += this->constants[i] * vl;
+      max += this->constants[i] * vh;
+    }
+
+  if (min > k || max < k)
+    return FD_NOSOLUTION;
+
+  if (min == max)
+    return FD_OK;
+
+  // XXX: poor man's propagation
+  if (min == k)
+    {
+      for (i = 0; i < terms; ++i)
+	{
+	  int c = this->constants[i];
+
+	  if (c && mins[i] != maxs[i])
+	    {
+	      if (c > 0)
+		_fd_var_del_gt(mins[i], VAR(this, i));
+	      else
+		_fd_var_del_lt(maxs[i], VAR(this, i));
+
+	      // check needed for when variables occur more than once
+	      // and with opposite signs
+	      if (fd_domain_empty(VAR(this, i)))
+		return FD_NOSOLUTION;
+
+	      _fd_revise_connected(this, VAR(this, i));
+	    }
+	}
+
+      fd__constraint_set_entailed(this);
+    }
+  else if (max == k)
+    {
+      for (i = 0; i < terms; ++i)
+	{
+	  int c = this->constants[i];
+
+	  if (c && mins[i] != maxs[i])
+	    {
+	      if (c > 0)
+		_fd_var_del_lt(maxs[i], VAR(this, i));
+	      else
+		_fd_var_del_gt(mins[i], VAR(this, i));
+
+	      if (fd_domain_empty(VAR(this, i)))
+		return FD_NOSOLUTION;
+
+	      _fd_revise_connected(this, VAR(this, i));
+	    }
+	}
+
+      fd__constraint_set_entailed(this);
+    }
+  else
+    {
+      if (max > k)
+	for (i = 0; i < terms; ++i)
+	  {
+	    int c = this->constants[i];
+	    int xmin = mins[i];
+	    int xmax = maxs[i];
+
+	    if (c > 0)
+	      {
+		// enforce sum{j!=i}(c[j] * min[j]) + c[i] * max[i] <= k
+
+		if ((xmax - xmin) * c > k - min)
+		  {
+		    // xmax = floor((k - min) / c) + xmin
+		    xmax = (k - min) / c + xmin;
+
+		    _fd_var_del_gt(xmax, VAR(this, i));
+
+		    if (fd_domain_empty(VAR(this, i)))
+		      return FD_NOSOLUTION;
+
+		    _fd_revise_connected(this, VAR(this, i));
+		  }
+	      }
+	    else if (c < 0)
+	      {
+		// enforce sum{j!=i}(c[j] * min[j]) + c[i] * min[i] <= k
+
+		if ((xmin - xmax) * c > k - min)
+		  {
+		    // xmin = ceil((k - min) / c) + xmax
+		    xmin = (k - min) / c + xmax;
+
+		    _fd_var_del_lt(xmin, VAR(this, i));
+
+		    if (fd_domain_empty(VAR(this, i)))
+		      return FD_NOSOLUTION;
+
+		    _fd_revise_connected(this, VAR(this, i));
+		  }
+	      }
+	  }
+
+      if (min < k)
+	for (i = 0; i < terms; ++i)
+	  {
+	    int c = this->constants[i];
+	    int xmin = mins[i];
+	    int xmax = maxs[i];
+
+	    if (c > 0)
+	      {
+		// enforce sum{j!=i}(c[j] * max[j]) + c[i] * min[i] >= k
+
+		if ((xmax - xmin) * c > max - k)
+		  {
+		    // xmin = ceil((k - max) / c) + xmax
+		    xmin = (k - max) / c + xmax;
+
+		    _fd_var_del_lt(xmin, VAR(this, i));
+
+		    if (fd_domain_empty(VAR(this, i))) // domains may have holes
+		      return FD_NOSOLUTION;
+
+		    _fd_revise_connected(this, VAR(this, i));
+		  }
+	      }
+	    else if (c < 0)
+	      {
+		// enforce sum{j!=i}(c[j] * max[j]) + c[i] * max[i] >= k
+
+		if ((xmax - xmin) * c < k - max)
+		  {
+		    // xmax = floor((k - max) / c) + xmin
+		    xmax = (k - max) / c + xmin;
+
+		    _fd_var_del_gt(xmax, VAR(this, i));
+
+		    if (fd_domain_empty(VAR(this, i))) // domains may have holes
+		      return FD_NOSOLUTION;
+
+		    _fd_revise_connected(this, VAR(this, i));
+		  }
+	      }
+	  }
+    }
+
+#ifdef CONSTRAINT_TEMPS
+  // save values
+  *base = min;
+  *(base + 1) = max;
+
+  fd__constraint_remember(this);
+#endif
+
+  return FD_OK;
+}
+
+static int fd_poly_eq_k_propagate2(fd_constraint this, fd_int culprit)
+{
+#ifdef CONSTRAINT_TEMPS
+  int k;
+  int min, max;
+  int terms = this->nvariables;
+  int *mins, *maxs;
+  int i;
+  int *base;
+  int x, c;
+  int nmin, nmin_x, nmax, nmax_x;
+
+  if (!fd__constraint_data_valid(this))
+    return fd_poly_eq_k_filter(this);	// ignores culprit
+
+  // bounds filtering
+
+  base = constraint_memory[this->index];
+
+  mins = base + 2;
+  maxs = mins + terms;
+
+  min = *base;
+  max = *(base + 1);
+
+  k = this->constants[terms];
+
+  // the culprit appears in one of the terms, find out which one(s)
+  for (x = 0; culprit != VAR(this, x); ++x)
+    ;
+
+  nmin_x = _fd_var_min(VAR(this, x));
+  nmax_x = _fd_var_max(VAR(this, x));
+
+  if (nmin_x == mins[x] && nmax_x == maxs[x])
+    return FD_OK;
+
+  nmin = min;
+  nmax = max;
+
+  do
+    {
+      c = this->constants[x];
+
+      if (c > 0)
+	{
+	  nmin = nmin + (nmin_x - mins[x]) * c;
+	  nmax = nmax - (maxs[x] - nmax_x) * c;
+	}
+      else if (c < 0)
+	{
+	  nmin = nmin - (maxs[x] - nmax_x) * c;
+	  nmax = nmax + (nmin_x - mins[x]) * c;
+	}
+
+      if (nmin > k || nmax < k)
+	return FD_NOSOLUTION;
+
+      mins[x] = nmin_x;
+      maxs[x] = nmax_x;
+
+      while (++x < terms && culprit != VAR(this, x))
+	;
+    }
+  while (x < terms);
+
+  // XXX: poor man's propagation
+  if (nmin == k)
+    {
+      for (i = 0; i < terms; ++i)
+	{
+	  int c = this->constants[i];
+
+	  if (c && mins[i] != maxs[i])
+	    {
+	      if (c > 0)
+		{
+		  _fd_var_del_gt(mins[i], VAR(this, i));
+
+		  maxs[i] = mins[i];
+		}
+	      else
+		{
+		  _fd_var_del_lt(maxs[i], VAR(this, i));
+
+		  mins[i] = maxs[i];
+		}
+
+	      if (fd_domain_empty(VAR(this, i)))
+		return FD_NOSOLUTION;
+
+	      _fd_revise_connected(this, VAR(this, i));
+	    }
+	}
+
+      fd__constraint_set_entailed(this);
+    }
+  else if (nmax == k)
+    {
+      for (i = 0; i < terms; ++i)
+	{
+	  int c = this->constants[i];
+
+	  if (c && mins[i] != maxs[i])
+	    {
+	      if (c > 0)
+		{
+		  _fd_var_del_lt(maxs[i], VAR(this, i));
+
+		  mins[i] = maxs[i];
+		}
+	      else
+		{
+		  _fd_var_del_gt(mins[i], VAR(this, i));
+
+		  maxs[i] = mins[i];
+		}
+
+	      if (fd_domain_empty(VAR(this, i)))
+		return FD_NOSOLUTION;
+
+	      _fd_revise_connected(this, VAR(this, i));
+	    }
+	}
+
+      fd__constraint_set_entailed(this);
+    }
+
+#if 0	// this is too expensive
+    {
+      if (nmax < max && nmax > k)
+	for (i = 0; i < terms; ++i)
+	  {
+	    int c = this->constants[i];
+	    int xmin = mins[i];
+	    int xmax = maxs[i];
+
+	    if (c > 0)
+	      {
+		// enforce sum{j!=i}(c[j] * min[j]) + c[i] * max[i] <= k
+
+		if ((xmax - xmin) * c > k - nmin)
+		  {
+		    // xmax = floor((k - nmin) / c) + xmin
+		    xmax = (k - nmin) / c + xmin;
+
+		    if (_fd_var_del_gt(xmax, VAR(this, i)))
+		      {
+			if (fd_domain_empty(VAR(this, i)))
+			  return FD_NOSOLUTION;
+
+			_fd_revise_connected(this, VAR(this, i));
+		      }
+
+		    maxs[i] = xmax;
+		  }
+	      }
+	    else if (c < 0)
+	      {
+		// enforce sum{j!=i}(c[j] * min[j]) + c[i] * min[i] <= k
+
+		if ((xmin - xmax) * c > k - nmin)
+		  {
+		    // xmin = ceil((k - nmin) / c) + xmax
+		    xmin = (k - nmin) / c + xmax;
+
+		    if (_fd_var_del_lt(xmin, VAR(this, i)))
+		      {
+			if (fd_domain_empty(VAR(this, i)))
+			  return FD_NOSOLUTION;
+
+			_fd_revise_connected(this, VAR(this, i));
+		      }
+
+		    mins[i] = xmin;
+		  }
+	      }
+	  }
+
+      if (nmin > min && nmin < k)
+	for (i = 0; i < terms; ++i)
+	  {
+	    int c = this->constants[i];
+	    int xmin = mins[i];
+	    int xmax = maxs[i];
+
+	    if (c > 0)
+	      {
+		// enforce sum{j!=i}(c[j] * max[j]) + c[i] * min[i] >= k
+
+		if ((xmax - xmin) * c > nmax - k)
+		  {
+		    // xmin = ceil((k - nmax) / c) + xmax
+		    xmin = (k - nmax) / c + xmax;
+
+		    if (_fd_var_del_lt(xmin, VAR(this, i)))
+		      {
+
+			if (fd_domain_empty(VAR(this, i))) // domains may have holes
+			  return FD_NOSOLUTION;
+
+			_fd_revise_connected(this, VAR(this, i));
+		      }
+
+		    mins[i] = xmin;
+		  }
+	      }
+	    else if (c < 0)
+	      {
+		// enforce sum{j!=i}(c[j] * max[j]) + c[i] * max[i] >= k
+
+		if ((xmax - xmin) * c < k - nmax)
+		  {
+		    // xmax = floor((k - nmax) / c) + xmin
+		    xmax = (k - nmax) / c + xmin;
+
+		    if (_fd_var_del_gt(xmax, VAR(this, i)))
+		      {
+			if (fd_domain_empty(VAR(this, i))) // domains may have holes
+			  return FD_NOSOLUTION;
+
+			_fd_revise_connected(this, VAR(this, i));
+		      }
+
+		    maxs[i] = xmax;
+		  }
+	      }
+	  }
+    }
+#endif
+
+  *base = nmin;
+  *(base + 1) = nmax;
+
+  return FD_OK;
+#else /* CONSTRAINT_TEMPS */
+  return fd_poly_eq_k_filter(this);	// ignores culprit
+#endif /* CONSTRAINT_TEMPS */
+}
+
+fd_constraint fd_poly_eq_k(int cs[], fd_int xs[], int nterms, int k)
+{
+  fd_constraint c = _fd_constraint_new(nterms, nterms + 1);
+  int i;
+
+  if (c)
+    {
+      for (i = 0; i < nterms; ++i)
+	c->variables[i] = FD_INT2C_VAR(xs[i]);
+      for (i = 0; i < nterms; ++i)
+	c->constants[i] = cs[i];
+      c->constants[nterms] = k;
+#ifdef CONSTRAINT_CLASS
+      c->kind = FD_CONSTR_POLY_EQ_K;
+#else /* CONSTRAINT_CLASS */
+      c->propagator2 = fd_poly_eq_k_propagate2;
+#endif /* CONSTRAINT_CLASS */
+
+      for (i = 0; i < c->nvariables; ++i)
+	_fd_var_add_constraint(VAR(c, i), c);
+
+      _fd_add_constraint(c);
+    }
+
+  return c;
+}
@@ -0,0 +1,250 @@
+/* poly-ne-k(C, X, k) == sum(C . X) != k */
+
+static int fd_poly_ne_k_filter(fd_constraint this)
+{
+  int k;
+  int min, max;
+  int terms = this->nvariables;
+  int *mins, *maxs;
+  int unset = 0;		// non-singleton variables
+  int i;
+#ifdef CONSTRAINT_TEMPS
+  int *base;
+
+  assert(!fd__constraint_data_valid(this));
+
+  if (!constraint_memory[this->index])
+    constraint_memory[this->index] = malloc((2 * terms + 3) * sizeof(int));
+
+  base = constraint_memory[this->index];
+
+  mins = base + 3;
+  maxs = mins + terms;
+#else
+  mins = alloca(terms * sizeof(*mins));
+  maxs = alloca(terms * sizeof(*maxs));
+#endif
+
+  k = this->constants[terms];
+
+  // sum the minima and the maxima of the terms
+  min = max = 0;
+
+  for (i = 0; i < terms; ++i)
+    {
+      int vl, vh;
+      int c = this->constants[i];
+
+      if (c > 0)
+	{
+	  vl = mins[i] = _fd_var_min(VAR(this, i));
+	  vh = maxs[i] = _fd_var_max(VAR(this, i));
+	}
+      else
+	{
+	  vl = maxs[i] = _fd_var_max(VAR(this, i));
+	  vh = mins[i] = _fd_var_min(VAR(this, i));
+	}
+
+      if (c && vl != vh)
+	unset++;
+
+      min += c * vl;
+      max += c * vh;
+    }
+
+  if (min > k || max < k)
+    {
+      fd__constraint_set_entailed(this);
+
+      return FD_OK;
+    }
+
+  if (min == max)
+    {
+      if (min == k)
+	return FD_NOSOLUTION;
+
+      fd__constraint_set_entailed(this);
+
+      return FD_OK;
+    }
+
+  if (unset == 1)
+    {
+      int m, v;
+
+      // find out which is the non-singleton variable
+      for (i = terms - 1; i >= 0; --i)
+	if (this->constants[i] && mins[i] != maxs[i])
+	  break;
+
+      if (this->constants[i] > 0)
+	m = min - this->constants[i] * mins[i];
+      else
+	m = min - this->constants[i] * maxs[i];
+
+      v = (k - m) / this->constants[i];
+
+      if (v * this->constants[i] + m == k)
+	fd_update_domain_and_check(del_val, v, VAR(this, i));
+
+      fd__constraint_set_entailed(this);
+    }
+
+#ifdef CONSTRAINT_TEMPS
+  // save values
+  *base = min;
+  *(base + 1) = max;
+  *(base + 2) = unset;
+
+  fd__constraint_remember(this);
+#endif
+
+  return FD_OK;
+}
+
+static int fd_poly_ne_k_propagate2(fd_constraint this, fd_int culprit)
+{
+#ifdef CONSTRAINT_TEMPS
+  int k;
+  int min, max;
+  int terms = this->nvariables;
+  int *mins, *maxs;
+  int unset;
+  int i;
+  int *base;
+  int x, c;
+  int nmin, nmin_x, nmax, nmax_x;
+
+  if (!fd__constraint_data_valid(this))
+    return fd_poly_ne_k_filter(this);	// ignores culprit
+
+  // bounds filtering
+
+  base = constraint_memory[this->index];
+
+  mins = base + 3;
+  maxs = mins + terms;
+
+  k = this->constants[terms];
+
+  min = *base;
+  max = *(base + 1);
+
+  unset = *(base + 2);
+
+  // the culprit appears in one of the terms, find out which one(s)
+  for (x = 0; culprit != VAR(this, x); ++x)
+    ;
+
+  nmin_x = _fd_var_min(VAR(this, x));
+  nmax_x = _fd_var_max(VAR(this, x));
+
+  if (nmin_x == mins[x] && nmax_x == maxs[x])
+    return FD_OK;
+
+  nmin = min;
+  nmax = max;
+
+  do
+    {
+      c = this->constants[x];
+
+      if (c && nmin_x == nmax_x)
+	unset--;
+
+      if (c > 0)
+	{
+	  nmin = nmin + (nmin_x - mins[x]) * c;
+	  nmax = nmax - (maxs[x] - nmax_x) * c;
+	}
+      else if (c < 0)
+	{
+	  nmin = nmin - (maxs[x] - nmax_x) * c;
+	  nmax = nmax + (nmin_x - mins[x]) * c;
+	}
+
+      mins[x] = nmin_x;
+      maxs[x] = nmax_x;
+
+      while (++x < terms && culprit != VAR(this, x))
+	;
+    }
+  while (x < terms);
+
+  if (nmin > k || nmax < k)
+    {
+      fd__constraint_set_entailed(this);
+
+      return FD_OK;
+    }
+
+  if (nmin == nmax)
+    {
+      if (nmin == k)
+	return FD_NOSOLUTION;
+
+      fd__constraint_set_entailed(this);
+
+      return FD_OK;
+    }
+
+  if (unset == 1)
+    {
+      int m, v;
+
+      // find out which is the non-singleton variable
+      for (i = terms - 1; i >= 0; --i)
+	if (this->constants[i] && mins[i] != maxs[i])
+	  break;
+
+      if (this->constants[i] > 0)
+	m = min - this->constants[i] * mins[i];
+      else
+	m = min - this->constants[i] * maxs[i];
+
+      v = (k - m) / this->constants[i];
+
+      if (v * this->constants[i] + m == k)
+	fd_update_domain_and_check(del_val, v, VAR(this, i));
+
+      fd__constraint_set_entailed(this);
+    }
+
+  *base = nmin;
+  *(base + 1) = nmax;
+  *(base + 2) = unset;
+
+  return FD_OK;
+#else /* CONSTRAINT_TEMPS */
+  return fd_poly_ne_k_filter(this);	// ignores culprit
+#endif /* CONSTRAINT_TEMPS */
+}
+
+fd_constraint fd_poly_ne_k(int cs[], fd_int xs[], int nterms, int k)
+{
+  fd_constraint c = _fd_constraint_new(nterms, nterms + 1);
+  int i;
+
+  if (c)
+    {
+      for (i = 0; i < nterms; ++i)
+	c->variables[i] = FD_INT2C_VAR(xs[i]);
+      for (i = 0; i < nterms; ++i)
+	c->constants[i] = cs[i];
+      c->constants[nterms] = k;
+#ifdef CONSTRAINT_CLASS
+      c->kind = FD_CONSTR_POLY_NE_K;
+#else /* CONSTRAINT_CLASS */
+      c->propagator2 = fd_poly_ne_k_propagate2;
+#endif /* CONSTRAINT_CLASS */
+
+      for (i = 0; i < c->nvariables; ++i)
+	_fd_var_add_constraint(VAR(c, i), c);
+
+      _fd_add_constraint(c);
+    }
+
+  return c;
+}
@@ -0,0 +1,288 @@
+/* poly-ne(C, X, y) == sum(C . X) != y */
+
+static int fd_poly_ne_filter(fd_constraint this)
+{
+  int ub, lb;
+  int min, max;
+  int terms = this->nconstants;
+  int *mins, *maxs;
+  int unset = 0;		// non-singleton variables
+  int i;
+#ifdef CONSTRAINT_TEMPS
+  int *base;
+
+  assert(!fd__constraint_data_valid(this));
+
+  if (!constraint_memory[this->index])
+    constraint_memory[this->index] = malloc((2 * terms + 5) * sizeof(int));
+
+  base = constraint_memory[this->index];
+
+  mins = base + 5;
+  maxs = mins + terms;
+#else
+  mins = alloca(terms * sizeof(*mins));
+  maxs = alloca(terms * sizeof(*maxs));
+#endif
+
+  lb = _fd_var_min(VAR(this, this->nvariables - 1));	// lower bound
+  ub = _fd_var_max(VAR(this, this->nvariables - 1));	// upper bound
+
+  // sum the minima and the maxima of the terms
+  min = max = 0;
+
+  for (i = 0; i < terms; ++i)
+    {
+      int vl, vh;
+      int c = this->constants[i];
+
+      if (c > 0)
+	{
+	  vl = mins[i] = _fd_var_min(VAR(this, i));
+	  vh = maxs[i] = _fd_var_max(VAR(this, i));
+	}
+      else
+	{
+	  vl = maxs[i] = _fd_var_max(VAR(this, i));
+	  vh = mins[i] = _fd_var_min(VAR(this, i));
+	}
+
+      if (c && vl != vh)
+	unset++;
+
+      min += c * vl;
+      max += c * vh;
+    }
+
+  if (min > ub || max < lb)
+    {
+      fd__constraint_set_entailed(this);
+
+      return FD_OK;
+    }
+
+  if (min == max)
+    {
+      fd_update_domain_and_check(del_val, min, VAR(this, this->nvariables - 1));
+
+      fd__constraint_set_entailed(this);
+    }
+  else if (unset == 1 && lb == ub)
+    {
+      int m, v;
+
+      // find out which is the non-singleton variable
+      for (i = terms - 1; i >= 0; --i)
+	if (this->constants[i] && mins[i] != maxs[i])
+	  break;
+
+      if (this->constants[i] > 0)
+	m = min - this->constants[i] * mins[i];
+      else
+	m = min - this->constants[i] * maxs[i];
+
+      v = (lb - m) / this->constants[i];
+
+      if (v * this->constants[i] + m == lb)
+	fd_update_domain_and_check(del_val, v, VAR(this, i));
+
+      fd__constraint_set_entailed(this);
+    }
+
+#ifdef CONSTRAINT_TEMPS
+  // save values
+  *base = lb;
+  *(base + 1) = ub;
+  *(base + 2) = min;
+  *(base + 3) = max;
+  *(base + 4) = unset;
+
+  fd__constraint_remember(this);
+#endif
+
+  return FD_OK;
+}
+
+static int fd_poly_ne_propagate2(fd_constraint this, fd_int culprit)
+{
+#ifdef CONSTRAINT_TEMPS
+  int ub, lb;
+  int min, max;
+  int terms = this->nconstants;
+  int *mins, *maxs;
+  int unset;
+  int i;
+  int *base;
+  int x, c;
+  int nmin, nmin_x, nmax, nmax_x;
+
+  if (!fd__constraint_data_valid(this))
+    return fd_poly_ne_filter(this);	// ignores culprit
+
+  // bounds filtering
+
+  base = constraint_memory[this->index];
+
+  mins = base + 5;
+  maxs = mins + terms;
+
+  lb = *base;
+  ub = *(base + 1);
+
+  min = *(base + 2);
+  max = *(base + 3);
+
+  unset = *(base + 4);
+
+  if (culprit == VAR(this, this->nvariables - 1))
+    {
+      // the culprit is the sum variable
+      int nlb, nub;
+
+      nlb = _fd_var_min(culprit);
+      nub = _fd_var_max(culprit);
+
+      if (nlb == lb && nub == ub)
+	return FD_OK;
+
+      if (min > nub || max < nlb)
+	fd__constraint_set_entailed(this);
+
+      if (unset == 1 && nlb == nub)
+	{
+	  int m, v;
+
+	  // find out which is the non-singleton variable
+	  for (i = this->nvariables - 2; i >= 0; --i)
+	    if (this->constants[i] && mins[i] != maxs[i])
+	      break;
+
+	  if (this->constants[i] > 0)
+	    m = min - this->constants[i] * mins[i];
+	  else
+	    m = min - this->constants[i] * maxs[i];
+
+	  v = (nlb - m) / this->constants[i];
+
+	  if (v * this->constants[i] + m == nlb)
+	    fd_update_domain_and_check(del_val, v, VAR(this, i));
+
+	  fd__constraint_set_entailed(this);
+	}
+
+      *base = nlb;
+      *(base + 1) = nub;
+
+      return FD_OK;
+    }
+
+  // the culprit appears in one of the terms, find out which one(s)
+  for (x = 0; culprit != VAR(this, x); ++x)
+    ;
+
+  nmin_x = _fd_var_min(VAR(this, x));
+  nmax_x = _fd_var_max(VAR(this, x));
+
+  if (nmin_x == mins[x] && nmax_x == maxs[x])
+    return FD_OK;
+
+  nmin = min;
+  nmax = max;
+
+  do
+    {
+      c = this->constants[x];
+
+      if (c && nmin_x == nmax_x)
+	unset--;
+
+      if (c > 0)
+	{
+	  nmin = nmin + (nmin_x - mins[x]) * c;
+	  nmax = nmax - (maxs[x] - nmax_x) * c;
+	}
+      else if (c < 0)
+	{
+	  nmin = nmin - (maxs[x] - nmax_x) * c;
+	  nmax = nmax + (nmin_x - mins[x]) * c;
+	}
+
+      mins[x] = nmin_x;
+      maxs[x] = nmax_x;
+
+      while (++x < terms && culprit != VAR(this, x))
+	;
+    }
+  while (x < terms);
+
+  if (nmin > ub || nmax < lb)
+    {
+      fd__constraint_set_entailed(this);
+
+      return FD_OK;
+    }
+
+  if (nmin == nmax)
+    {
+      fd_update_domain_and_check(del_val, nmin, VAR(this, this->nvariables - 1));
+
+      fd__constraint_set_entailed(this);
+    }
+  else if (unset == 1 && lb == ub)
+    {
+      int m, v;
+
+      // find out which is the non-singleton variable
+      for (i = terms - 1; i >= 0; --i)
+	if (this->constants[i] && mins[i] != maxs[i])
+	  break;
+
+      if (this->constants[i] > 0)
+	m = min - this->constants[i] * mins[i];
+      else
+	m = min - this->constants[i] * maxs[i];
+
+      v = (lb - m) / this->constants[i];
+
+      if (v * this->constants[i] + m == lb)
+	fd_update_domain_and_check(del_val, v, VAR(this, i));
+
+      fd__constraint_set_entailed(this);
+    }
+
+  *(base + 2) = nmin;
+  *(base + 3) = nmax;
+  *(base + 4) = unset;
+
+  return FD_OK;
+#else /* CONSTRAINT_TEMPS */
+  return fd_poly_ne_filter(this);	// ignores culprit
+#endif /* CONSTRAINT_TEMPS */
+}
+
+fd_constraint fd_poly_ne(int cs[], fd_int xs[], int nterms, fd_int y)
+{
+  fd_constraint c = _fd_constraint_new(nterms + 1, nterms);
+  int i;
+
+  if (c)
+    {
+      for (i = 0; i < nterms; ++i)
+	c->variables[i] = FD_INT2C_VAR(xs[i]);
+      c->variables[nterms] = FD_INT2C_VAR(y);
+      for (i = 0; i < nterms; ++i)
+	c->constants[i] = cs[i];
+#ifdef CONSTRAINT_CLASS
+      c->kind = FD_CONSTR_POLY_NE;
+#else /* CONSTRAINT_CLASS */
+      c->propagator2 = fd_poly_ne_propagate2;
+#endif /* CONSTRAINT_CLASS */
+
+      for (i = 0; i < c->nvariables; ++i)
+	_fd_var_add_constraint(VAR(c, i), c);
+
+      _fd_add_constraint(c);
+    }
+
+  return c;
+}
...	...	@@ -173,6 +173,9 @@ Global constraints
173	173	fd_sum_prod(fd_int X[], fd_int Y[], int n, int k) X . Y = k
174	174
175	175	fd_poly_eq(int C[], fd_int Y[], int n, fd_int z) z <= C . Y <= z
	176	+ fd_poly_eq_k(int C[], fd_int Y[], int n, int k) C . Y = k
	177	+ fd_poly_ne(int C[], fd_int Y[], int n, fd_int z) C . Y != z
	178	+ fd_poly_ne_k(int C[], fd_int Y[], int n, int k) C . Y != k
176	179	fd_knapsack2(fd_int X[], fd_int Y[], int n, fd_int z) z <= X . Y <= z
177	180
178	181	fd_exactly(fd_int X[], int n, int c, int k) #{ i \| X[i] = k } = c
...	...
...	...	@@ -149,6 +149,9 @@ void fd__constraint_free_memory()
149	149	#include <constraints/poly-eq.c>
150	150	#include <constraints/element-var.c>
151	151	#include <constraints/exactly-vars.c>
	152	+#include <constraints/poly-eq-k.c>
	153	+#include <constraints/poly-ne.c>
	154	+#include <constraints/poly-ne-k.c>
152	155
153	156	#ifdef CONSTRAINT_CLASS
154	157
...	...	@@ -178,6 +181,9 @@ _FD_CONSTRAINT_INITIALISATION(max, FD_CONSTR_MAX)
178	181	_FD_CONSTRAINT_INITIALISATION(poly_eq, FD_CONSTR_POLY_EQ)
179	182	_FD_CONSTRAINT_INITIALISATION(element_var, FD_CONSTR_ELEMENT_VAR)
180	183	_FD_CONSTRAINT_INITIALISATION(exactly_vars, FD_CONSTR_EXACTLY_VARS)
	184	+_FD_CONSTRAINT_INITIALISATION(poly_eq_k, FD_CONSTR_POLY_EQ_K)
	185	+_FD_CONSTRAINT_INITIALISATION(poly_ne, FD_CONSTR_POLY_NE)
	186	+_FD_CONSTRAINT_INITIALISATION(poly_ne_k, FD_CONSTR_POLY_NE_K)
181	187
182	188	void _fd_init_constraints()
183	189	{
...	...	@@ -212,6 +218,9 @@ void _fd_init_constraints()
212	218	_fd_init_constraint(poly_eq);
213	219	_fd_init_constraint(element_var);
214	220	_fd_init_constraint(exactly_vars);
	221	+ _fd_init_constraint(poly_eq_k);
	222	+ _fd_init_constraint(poly_ne);
	223	+ _fd_init_constraint(poly_ne_k);
215	224
216	225	done = 1;
217	226	}
...	...
...	...	@@ -33,6 +33,9 @@ typedef enum {
33	33	FD_CONSTR_POLY_EQ,
34	34	FD_CONSTR_ELEMENT_VAR,
35	35	FD_CONSTR_EXACTLY_VARS,
	36	+ FD_CONSTR_POLY_EQ_K,
	37	+ FD_CONSTR_POLY_NE,
	38	+ FD_CONSTR_POLY_NE_K,
36	39	FD_CONSTR_KINDS // number of kinds of constraints
37	40	} _fd_constraint_kind;
38	41
...	...
...	...	@@ -0,0 +1,469 @@
	1	+/* poly-eq-k(C, X, k) == sum(C . X) = k */
	2	+
	3	+static int fd_poly_eq_k_filter(fd_constraint this)
	4	+{
	5	+ int k;
	6	+ int min, max;
	7	+ int terms = this->nvariables;
	8	+ int mins, maxs;
	9	+ int i;
	10	+#ifdef CONSTRAINT_TEMPS
	11	+ int *base;
	12	+
	13	+ assert(!fd__constraint_data_valid(this));
	14	+
	15	+ if (!constraint_memory[this->index])
	16	+ constraint_memory[this->index] = malloc((2 * terms + 2) * sizeof(int));
	17	+
	18	+ base = constraint_memory[this->index];
	19	+
	20	+ mins = base + 2;
	21	+ maxs = mins + terms;
	22	+#else
	23	+ mins = alloca(terms * sizeof(*mins));
	24	+ maxs = alloca(terms * sizeof(*maxs));
	25	+#endif
	26	+
	27	+ k = this->constants[terms];
	28	+
	29	+ // sum the minima and the maxima of the terms
	30	+ min = max = 0;
	31	+
	32	+ for (i = 0; i < terms; ++i)
	33	+ {
	34	+ int vl, vh;
	35	+
	36	+ if (this->constants[i] > 0)
	37	+ {
	38	+ vl = mins[i] = _fd_var_min(VAR(this, i));
	39	+ vh = maxs[i] = _fd_var_max(VAR(this, i));
	40	+ }
	41	+ else
	42	+ {
	43	+ vl = maxs[i] = _fd_var_max(VAR(this, i));
	44	+ vh = mins[i] = _fd_var_min(VAR(this, i));
	45	+ }
	46	+
	47	+ min += this->constants[i] * vl;
	48	+ max += this->constants[i] * vh;
	49	+ }
	50	+
	51	+ if (min > k \|\| max < k)
	52	+ return FD_NOSOLUTION;
	53	+
	54	+ if (min == max)
	55	+ return FD_OK;
	56	+
	57	+ // XXX: poor man's propagation
	58	+ if (min == k)
	59	+ {
	60	+ for (i = 0; i < terms; ++i)
	61	+ {
	62	+ int c = this->constants[i];
	63	+
	64	+ if (c && mins[i] != maxs[i])
	65	+ {
	66	+ if (c > 0)
	67	+ _fd_var_del_gt(mins[i], VAR(this, i));
	68	+ else
	69	+ _fd_var_del_lt(maxs[i], VAR(this, i));
	70	+
	71	+ // check needed for when variables occur more than once
	72	+ // and with opposite signs
	73	+ if (fd_domain_empty(VAR(this, i)))
	74	+ return FD_NOSOLUTION;
	75	+
	76	+ _fd_revise_connected(this, VAR(this, i));
	77	+ }
	78	+ }
	79	+
	80	+ fd__constraint_set_entailed(this);
	81	+ }
	82	+ else if (max == k)
	83	+ {
	84	+ for (i = 0; i < terms; ++i)
	85	+ {
	86	+ int c = this->constants[i];
	87	+
	88	+ if (c && mins[i] != maxs[i])
	89	+ {
	90	+ if (c > 0)
	91	+ _fd_var_del_lt(maxs[i], VAR(this, i));
	92	+ else
	93	+ _fd_var_del_gt(mins[i], VAR(this, i));
	94	+
	95	+ if (fd_domain_empty(VAR(this, i)))
	96	+ return FD_NOSOLUTION;
	97	+
	98	+ _fd_revise_connected(this, VAR(this, i));
	99	+ }
	100	+ }
	101	+
	102	+ fd__constraint_set_entailed(this);
	103	+ }
	104	+ else
	105	+ {
	106	+ if (max > k)
	107	+ for (i = 0; i < terms; ++i)
	108	+ {
	109	+ int c = this->constants[i];
	110	+ int xmin = mins[i];
	111	+ int xmax = maxs[i];
	112	+
	113	+ if (c > 0)
	114	+ {
	115	+ // enforce sum{j!=i}(c[j] * min[j]) + c[i] * max[i] <= k
	116	+
	117	+ if ((xmax - xmin) * c > k - min)
	118	+ {
	119	+ // xmax = floor((k - min) / c) + xmin
	120	+ xmax = (k - min) / c + xmin;
	121	+
	122	+ _fd_var_del_gt(xmax, VAR(this, i));
	123	+
	124	+ if (fd_domain_empty(VAR(this, i)))
	125	+ return FD_NOSOLUTION;
	126	+
	127	+ _fd_revise_connected(this, VAR(this, i));
	128	+ }
	129	+ }
	130	+ else if (c < 0)
	131	+ {
	132	+ // enforce sum{j!=i}(c[j] * min[j]) + c[i] * min[i] <= k
	133	+
	134	+ if ((xmin - xmax) * c > k - min)
	135	+ {
	136	+ // xmin = ceil((k - min) / c) + xmax
	137	+ xmin = (k - min) / c + xmax;
	138	+
	139	+ _fd_var_del_lt(xmin, VAR(this, i));
	140	+
	141	+ if (fd_domain_empty(VAR(this, i)))
	142	+ return FD_NOSOLUTION;
	143	+
	144	+ _fd_revise_connected(this, VAR(this, i));
	145	+ }
	146	+ }
	147	+ }
	148	+
	149	+ if (min < k)
	150	+ for (i = 0; i < terms; ++i)
	151	+ {
	152	+ int c = this->constants[i];
	153	+ int xmin = mins[i];
	154	+ int xmax = maxs[i];
	155	+
	156	+ if (c > 0)
	157	+ {
	158	+ // enforce sum{j!=i}(c[j] * max[j]) + c[i] * min[i] >= k
	159	+
	160	+ if ((xmax - xmin) * c > max - k)
	161	+ {
	162	+ // xmin = ceil((k - max) / c) + xmax
	163	+ xmin = (k - max) / c + xmax;
	164	+
	165	+ _fd_var_del_lt(xmin, VAR(this, i));
	166	+
	167	+ if (fd_domain_empty(VAR(this, i))) // domains may have holes
	168	+ return FD_NOSOLUTION;
	169	+
	170	+ _fd_revise_connected(this, VAR(this, i));
	171	+ }
	172	+ }
	173	+ else if (c < 0)
	174	+ {
	175	+ // enforce sum{j!=i}(c[j] * max[j]) + c[i] * max[i] >= k
	176	+
	177	+ if ((xmax - xmin) * c < k - max)
	178	+ {
	179	+ // xmax = floor((k - max) / c) + xmin
	180	+ xmax = (k - max) / c + xmin;
	181	+
	182	+ _fd_var_del_gt(xmax, VAR(this, i));
	183	+
	184	+ if (fd_domain_empty(VAR(this, i))) // domains may have holes
	185	+ return FD_NOSOLUTION;
	186	+
	187	+ _fd_revise_connected(this, VAR(this, i));
	188	+ }
	189	+ }
	190	+ }
	191	+ }
	192	+
	193	+#ifdef CONSTRAINT_TEMPS
	194	+ // save values
	195	+ *base = min;
	196	+ *(base + 1) = max;
	197	+
	198	+ fd__constraint_remember(this);
	199	+#endif
	200	+
	201	+ return FD_OK;
	202	+}
	203	+
	204	+static int fd_poly_eq_k_propagate2(fd_constraint this, fd_int culprit)
	205	+{
	206	+#ifdef CONSTRAINT_TEMPS
	207	+ int k;
	208	+ int min, max;
	209	+ int terms = this->nvariables;
	210	+ int mins, maxs;
	211	+ int i;
	212	+ int *base;
	213	+ int x, c;
	214	+ int nmin, nmin_x, nmax, nmax_x;
	215	+
	216	+ if (!fd__constraint_data_valid(this))
	217	+ return fd_poly_eq_k_filter(this); // ignores culprit
	218	+
	219	+ // bounds filtering
	220	+
	221	+ base = constraint_memory[this->index];
	222	+
	223	+ mins = base + 2;
	224	+ maxs = mins + terms;
	225	+
	226	+ min = *base;
	227	+ max = *(base + 1);
	228	+
	229	+ k = this->constants[terms];
	230	+
	231	+ // the culprit appears in one of the terms, find out which one(s)
	232	+ for (x = 0; culprit != VAR(this, x); ++x)
	233	+ ;
	234	+
	235	+ nmin_x = _fd_var_min(VAR(this, x));
	236	+ nmax_x = _fd_var_max(VAR(this, x));
	237	+
	238	+ if (nmin_x == mins[x] && nmax_x == maxs[x])
	239	+ return FD_OK;
	240	+
	241	+ nmin = min;
	242	+ nmax = max;
	243	+
	244	+ do
	245	+ {
	246	+ c = this->constants[x];
	247	+
	248	+ if (c > 0)
	249	+ {
	250	+ nmin = nmin + (nmin_x - mins[x]) * c;
	251	+ nmax = nmax - (maxs[x] - nmax_x) * c;
	252	+ }
	253	+ else if (c < 0)
	254	+ {
	255	+ nmin = nmin - (maxs[x] - nmax_x) * c;
	256	+ nmax = nmax + (nmin_x - mins[x]) * c;
	257	+ }
	258	+
	259	+ if (nmin > k \|\| nmax < k)
	260	+ return FD_NOSOLUTION;
	261	+
	262	+ mins[x] = nmin_x;
	263	+ maxs[x] = nmax_x;
	264	+
	265	+ while (++x < terms && culprit != VAR(this, x))
	266	+ ;
	267	+ }
	268	+ while (x < terms);
	269	+
	270	+ // XXX: poor man's propagation
	271	+ if (nmin == k)
	272	+ {
	273	+ for (i = 0; i < terms; ++i)
	274	+ {
	275	+ int c = this->constants[i];
	276	+
	277	+ if (c && mins[i] != maxs[i])
	278	+ {
	279	+ if (c > 0)
	280	+ {
	281	+ _fd_var_del_gt(mins[i], VAR(this, i));
	282	+
	283	+ maxs[i] = mins[i];
	284	+ }
	285	+ else
	286	+ {
	287	+ _fd_var_del_lt(maxs[i], VAR(this, i));
	288	+
	289	+ mins[i] = maxs[i];
	290	+ }
	291	+
	292	+ if (fd_domain_empty(VAR(this, i)))
	293	+ return FD_NOSOLUTION;
	294	+
	295	+ _fd_revise_connected(this, VAR(this, i));
	296	+ }
	297	+ }
	298	+
	299	+ fd__constraint_set_entailed(this);
	300	+ }
	301	+ else if (nmax == k)
	302	+ {
	303	+ for (i = 0; i < terms; ++i)
	304	+ {
	305	+ int c = this->constants[i];
	306	+
	307	+ if (c && mins[i] != maxs[i])
	308	+ {
	309	+ if (c > 0)
	310	+ {
	311	+ _fd_var_del_lt(maxs[i], VAR(this, i));
	312	+
	313	+ mins[i] = maxs[i];
	314	+ }
	315	+ else
	316	+ {
	317	+ _fd_var_del_gt(mins[i], VAR(this, i));
	318	+
	319	+ maxs[i] = mins[i];
	320	+ }
	321	+
	322	+ if (fd_domain_empty(VAR(this, i)))
	323	+ return FD_NOSOLUTION;
	324	+
	325	+ _fd_revise_connected(this, VAR(this, i));
	326	+ }
	327	+ }
	328	+
	329	+ fd__constraint_set_entailed(this);
	330	+ }
	331	+
	332	+#if 0 // this is too expensive
	333	+ {
	334	+ if (nmax < max && nmax > k)
	335	+ for (i = 0; i < terms; ++i)
	336	+ {
	337	+ int c = this->constants[i];
	338	+ int xmin = mins[i];
	339	+ int xmax = maxs[i];
	340	+
	341	+ if (c > 0)
	342	+ {
	343	+ // enforce sum{j!=i}(c[j] * min[j]) + c[i] * max[i] <= k
	344	+
	345	+ if ((xmax - xmin) * c > k - nmin)
	346	+ {
	347	+ // xmax = floor((k - nmin) / c) + xmin
	348	+ xmax = (k - nmin) / c + xmin;
	349	+
	350	+ if (_fd_var_del_gt(xmax, VAR(this, i)))
	351	+ {
	352	+ if (fd_domain_empty(VAR(this, i)))
	353	+ return FD_NOSOLUTION;
	354	+
	355	+ _fd_revise_connected(this, VAR(this, i));
	356	+ }
	357	+
	358	+ maxs[i] = xmax;
	359	+ }
	360	+ }
	361	+ else if (c < 0)
	362	+ {
	363	+ // enforce sum{j!=i}(c[j] * min[j]) + c[i] * min[i] <= k
	364	+
	365	+ if ((xmin - xmax) * c > k - nmin)
	366	+ {
	367	+ // xmin = ceil((k - nmin) / c) + xmax
	368	+ xmin = (k - nmin) / c + xmax;
	369	+
	370	+ if (_fd_var_del_lt(xmin, VAR(this, i)))
	371	+ {
	372	+ if (fd_domain_empty(VAR(this, i)))
	373	+ return FD_NOSOLUTION;
	374	+
	375	+ _fd_revise_connected(this, VAR(this, i));
	376	+ }
	377	+
	378	+ mins[i] = xmin;
	379	+ }
	380	+ }
	381	+ }
	382	+
	383	+ if (nmin > min && nmin < k)
	384	+ for (i = 0; i < terms; ++i)
	385	+ {
	386	+ int c = this->constants[i];
	387	+ int xmin = mins[i];
	388	+ int xmax = maxs[i];
	389	+
	390	+ if (c > 0)
	391	+ {
	392	+ // enforce sum{j!=i}(c[j] * max[j]) + c[i] * min[i] >= k
	393	+
	394	+ if ((xmax - xmin) * c > nmax - k)
	395	+ {
	396	+ // xmin = ceil((k - nmax) / c) + xmax
	397	+ xmin = (k - nmax) / c + xmax;
	398	+
	399	+ if (_fd_var_del_lt(xmin, VAR(this, i)))
	400	+ {
	401	+
	402	+ if (fd_domain_empty(VAR(this, i))) // domains may have holes
	403	+ return FD_NOSOLUTION;
	404	+
	405	+ _fd_revise_connected(this, VAR(this, i));
	406	+ }
	407	+
	408	+ mins[i] = xmin;
	409	+ }
	410	+ }
	411	+ else if (c < 0)
	412	+ {
	413	+ // enforce sum{j!=i}(c[j] * max[j]) + c[i] * max[i] >= k
	414	+
	415	+ if ((xmax - xmin) * c < k - nmax)
	416	+ {
	417	+ // xmax = floor((k - nmax) / c) + xmin
	418	+ xmax = (k - nmax) / c + xmin;
	419	+
	420	+ if (_fd_var_del_gt(xmax, VAR(this, i)))
	421	+ {
	422	+ if (fd_domain_empty(VAR(this, i))) // domains may have holes
	423	+ return FD_NOSOLUTION;
	424	+
	425	+ _fd_revise_connected(this, VAR(this, i));
	426	+ }
	427	+
	428	+ maxs[i] = xmax;
	429	+ }
	430	+ }
	431	+ }
	432	+ }
	433	+#endif
	434	+
	435	+ *base = nmin;
	436	+ *(base + 1) = nmax;
	437	+
	438	+ return FD_OK;
	439	+#else /* CONSTRAINT_TEMPS */
	440	+ return fd_poly_eq_k_filter(this); // ignores culprit
	441	+#endif /* CONSTRAINT_TEMPS */
	442	+}
	443	+
	444	+fd_constraint fd_poly_eq_k(int cs[], fd_int xs[], int nterms, int k)
	445	+{
	446	+ fd_constraint c = _fd_constraint_new(nterms, nterms + 1);
	447	+ int i;
	448	+
	449	+ if (c)
	450	+ {
	451	+ for (i = 0; i < nterms; ++i)
	452	+ c->variables[i] = FD_INT2C_VAR(xs[i]);
	453	+ for (i = 0; i < nterms; ++i)
	454	+ c->constants[i] = cs[i];
	455	+ c->constants[nterms] = k;
	456	+#ifdef CONSTRAINT_CLASS
	457	+ c->kind = FD_CONSTR_POLY_EQ_K;
	458	+#else /* CONSTRAINT_CLASS */
	459	+ c->propagator2 = fd_poly_eq_k_propagate2;
	460	+#endif /* CONSTRAINT_CLASS */
	461	+
	462	+ for (i = 0; i < c->nvariables; ++i)
	463	+ _fd_var_add_constraint(VAR(c, i), c);
	464	+
	465	+ _fd_add_constraint(c);
	466	+ }
	467	+
	468	+ return c;
	469	+}
...	...
...	...	@@ -0,0 +1,250 @@
	1	+/* poly-ne-k(C, X, k) == sum(C . X) != k */
	2	+
	3	+static int fd_poly_ne_k_filter(fd_constraint this)
	4	+{
	5	+ int k;
	6	+ int min, max;
	7	+ int terms = this->nvariables;
	8	+ int mins, maxs;
	9	+ int unset = 0; // non-singleton variables
	10	+ int i;
	11	+#ifdef CONSTRAINT_TEMPS
	12	+ int *base;
	13	+
	14	+ assert(!fd__constraint_data_valid(this));
	15	+
	16	+ if (!constraint_memory[this->index])
	17	+ constraint_memory[this->index] = malloc((2 * terms + 3) * sizeof(int));
	18	+
	19	+ base = constraint_memory[this->index];
	20	+
	21	+ mins = base + 3;
	22	+ maxs = mins + terms;
	23	+#else
	24	+ mins = alloca(terms * sizeof(*mins));
	25	+ maxs = alloca(terms * sizeof(*maxs));
	26	+#endif
	27	+
	28	+ k = this->constants[terms];
	29	+
	30	+ // sum the minima and the maxima of the terms
	31	+ min = max = 0;
	32	+
	33	+ for (i = 0; i < terms; ++i)
	34	+ {
	35	+ int vl, vh;
	36	+ int c = this->constants[i];
	37	+
	38	+ if (c > 0)
	39	+ {
	40	+ vl = mins[i] = _fd_var_min(VAR(this, i));
	41	+ vh = maxs[i] = _fd_var_max(VAR(this, i));
	42	+ }
	43	+ else
	44	+ {
	45	+ vl = maxs[i] = _fd_var_max(VAR(this, i));
	46	+ vh = mins[i] = _fd_var_min(VAR(this, i));
	47	+ }
	48	+
	49	+ if (c && vl != vh)
	50	+ unset++;
	51	+
	52	+ min += c * vl;
	53	+ max += c * vh;
	54	+ }
	55	+
	56	+ if (min > k \|\| max < k)
	57	+ {
	58	+ fd__constraint_set_entailed(this);
	59	+
	60	+ return FD_OK;
	61	+ }
	62	+
	63	+ if (min == max)
	64	+ {
	65	+ if (min == k)
	66	+ return FD_NOSOLUTION;
	67	+
	68	+ fd__constraint_set_entailed(this);
	69	+
	70	+ return FD_OK;
	71	+ }
	72	+
	73	+ if (unset == 1)
	74	+ {
	75	+ int m, v;
	76	+
	77	+ // find out which is the non-singleton variable
	78	+ for (i = terms - 1; i >= 0; --i)
	79	+ if (this->constants[i] && mins[i] != maxs[i])
	80	+ break;
	81	+
	82	+ if (this->constants[i] > 0)
	83	+ m = min - this->constants[i] * mins[i];
	84	+ else
	85	+ m = min - this->constants[i] * maxs[i];
	86	+
	87	+ v = (k - m) / this->constants[i];
	88	+
	89	+ if (v * this->constants[i] + m == k)
	90	+ fd_update_domain_and_check(del_val, v, VAR(this, i));
	91	+
	92	+ fd__constraint_set_entailed(this);
	93	+ }
	94	+
	95	+#ifdef CONSTRAINT_TEMPS
	96	+ // save values
	97	+ *base = min;
	98	+ *(base + 1) = max;
	99	+ *(base + 2) = unset;
	100	+
	101	+ fd__constraint_remember(this);
	102	+#endif
	103	+
	104	+ return FD_OK;
	105	+}
	106	+
	107	+static int fd_poly_ne_k_propagate2(fd_constraint this, fd_int culprit)
	108	+{
	109	+#ifdef CONSTRAINT_TEMPS
	110	+ int k;
	111	+ int min, max;
	112	+ int terms = this->nvariables;
	113	+ int mins, maxs;
	114	+ int unset;
	115	+ int i;
	116	+ int *base;
	117	+ int x, c;
	118	+ int nmin, nmin_x, nmax, nmax_x;
	119	+
	120	+ if (!fd__constraint_data_valid(this))
	121	+ return fd_poly_ne_k_filter(this); // ignores culprit
	122	+
	123	+ // bounds filtering
	124	+
	125	+ base = constraint_memory[this->index];
	126	+
	127	+ mins = base + 3;
	128	+ maxs = mins + terms;
	129	+
	130	+ k = this->constants[terms];
	131	+
	132	+ min = *base;
	133	+ max = *(base + 1);
	134	+
	135	+ unset = *(base + 2);
	136	+
	137	+ // the culprit appears in one of the terms, find out which one(s)
	138	+ for (x = 0; culprit != VAR(this, x); ++x)
	139	+ ;
	140	+
	141	+ nmin_x = _fd_var_min(VAR(this, x));
	142	+ nmax_x = _fd_var_max(VAR(this, x));
	143	+
	144	+ if (nmin_x == mins[x] && nmax_x == maxs[x])
	145	+ return FD_OK;
	146	+
	147	+ nmin = min;
	148	+ nmax = max;
	149	+
	150	+ do
	151	+ {
	152	+ c = this->constants[x];
	153	+
	154	+ if (c && nmin_x == nmax_x)
	155	+ unset--;
	156	+
	157	+ if (c > 0)
	158	+ {
	159	+ nmin = nmin + (nmin_x - mins[x]) * c;
	160	+ nmax = nmax - (maxs[x] - nmax_x) * c;
	161	+ }
	162	+ else if (c < 0)
	163	+ {
	164	+ nmin = nmin - (maxs[x] - nmax_x) * c;
	165	+ nmax = nmax + (nmin_x - mins[x]) * c;
	166	+ }
	167	+
	168	+ mins[x] = nmin_x;
	169	+ maxs[x] = nmax_x;
	170	+
	171	+ while (++x < terms && culprit != VAR(this, x))
	172	+ ;
	173	+ }
	174	+ while (x < terms);
	175	+
	176	+ if (nmin > k \|\| nmax < k)
	177	+ {
	178	+ fd__constraint_set_entailed(this);
	179	+
	180	+ return FD_OK;
	181	+ }
	182	+
	183	+ if (nmin == nmax)
	184	+ {
	185	+ if (nmin == k)
	186	+ return FD_NOSOLUTION;
	187	+
	188	+ fd__constraint_set_entailed(this);
	189	+
	190	+ return FD_OK;
	191	+ }
	192	+
	193	+ if (unset == 1)
	194	+ {
	195	+ int m, v;
	196	+
	197	+ // find out which is the non-singleton variable
	198	+ for (i = terms - 1; i >= 0; --i)
	199	+ if (this->constants[i] && mins[i] != maxs[i])
	200	+ break;
	201	+
	202	+ if (this->constants[i] > 0)
	203	+ m = min - this->constants[i] * mins[i];
	204	+ else
	205	+ m = min - this->constants[i] * maxs[i];
	206	+
	207	+ v = (k - m) / this->constants[i];
	208	+
	209	+ if (v * this->constants[i] + m == k)
	210	+ fd_update_domain_and_check(del_val, v, VAR(this, i));
	211	+
	212	+ fd__constraint_set_entailed(this);
	213	+ }
	214	+
	215	+ *base = nmin;
	216	+ *(base + 1) = nmax;
	217	+ *(base + 2) = unset;
	218	+
	219	+ return FD_OK;
	220	+#else /* CONSTRAINT_TEMPS */
	221	+ return fd_poly_ne_k_filter(this); // ignores culprit
	222	+#endif /* CONSTRAINT_TEMPS */
	223	+}
	224	+
	225	+fd_constraint fd_poly_ne_k(int cs[], fd_int xs[], int nterms, int k)
	226	+{
	227	+ fd_constraint c = _fd_constraint_new(nterms, nterms + 1);
	228	+ int i;
	229	+
	230	+ if (c)
	231	+ {
	232	+ for (i = 0; i < nterms; ++i)
	233	+ c->variables[i] = FD_INT2C_VAR(xs[i]);
	234	+ for (i = 0; i < nterms; ++i)
	235	+ c->constants[i] = cs[i];
	236	+ c->constants[nterms] = k;
	237	+#ifdef CONSTRAINT_CLASS
	238	+ c->kind = FD_CONSTR_POLY_NE_K;
	239	+#else /* CONSTRAINT_CLASS */
	240	+ c->propagator2 = fd_poly_ne_k_propagate2;
	241	+#endif /* CONSTRAINT_CLASS */
	242	+
	243	+ for (i = 0; i < c->nvariables; ++i)
	244	+ _fd_var_add_constraint(VAR(c, i), c);
	245	+
	246	+ _fd_add_constraint(c);
	247	+ }
	248	+
	249	+ return c;
	250	+}
...	...
...	...	@@ -0,0 +1,288 @@
	1	+/* poly-ne(C, X, y) == sum(C . X) != y */
	2	+
	3	+static int fd_poly_ne_filter(fd_constraint this)
	4	+{
	5	+ int ub, lb;
	6	+ int min, max;
	7	+ int terms = this->nconstants;
	8	+ int mins, maxs;
	9	+ int unset = 0; // non-singleton variables
	10	+ int i;
	11	+#ifdef CONSTRAINT_TEMPS
	12	+ int *base;
	13	+
	14	+ assert(!fd__constraint_data_valid(this));
	15	+
	16	+ if (!constraint_memory[this->index])
	17	+ constraint_memory[this->index] = malloc((2 * terms + 5) * sizeof(int));
	18	+
	19	+ base = constraint_memory[this->index];
	20	+
	21	+ mins = base + 5;
	22	+ maxs = mins + terms;
	23	+#else
	24	+ mins = alloca(terms * sizeof(*mins));
	25	+ maxs = alloca(terms * sizeof(*maxs));
	26	+#endif
	27	+
	28	+ lb = _fd_var_min(VAR(this, this->nvariables - 1)); // lower bound
	29	+ ub = _fd_var_max(VAR(this, this->nvariables - 1)); // upper bound
	30	+
	31	+ // sum the minima and the maxima of the terms
	32	+ min = max = 0;
	33	+
	34	+ for (i = 0; i < terms; ++i)
	35	+ {
	36	+ int vl, vh;
	37	+ int c = this->constants[i];
	38	+
	39	+ if (c > 0)
	40	+ {
	41	+ vl = mins[i] = _fd_var_min(VAR(this, i));
	42	+ vh = maxs[i] = _fd_var_max(VAR(this, i));
	43	+ }
	44	+ else
	45	+ {
	46	+ vl = maxs[i] = _fd_var_max(VAR(this, i));
	47	+ vh = mins[i] = _fd_var_min(VAR(this, i));
	48	+ }
	49	+
	50	+ if (c && vl != vh)
	51	+ unset++;
	52	+
	53	+ min += c * vl;
	54	+ max += c * vh;
	55	+ }
	56	+
	57	+ if (min > ub \|\| max < lb)
	58	+ {
	59	+ fd__constraint_set_entailed(this);
	60	+
	61	+ return FD_OK;
	62	+ }
	63	+
	64	+ if (min == max)
	65	+ {
	66	+ fd_update_domain_and_check(del_val, min, VAR(this, this->nvariables - 1));
	67	+
	68	+ fd__constraint_set_entailed(this);
	69	+ }
	70	+ else if (unset == 1 && lb == ub)
	71	+ {
	72	+ int m, v;
	73	+
	74	+ // find out which is the non-singleton variable
	75	+ for (i = terms - 1; i >= 0; --i)
	76	+ if (this->constants[i] && mins[i] != maxs[i])
	77	+ break;
	78	+
	79	+ if (this->constants[i] > 0)
	80	+ m = min - this->constants[i] * mins[i];
	81	+ else
	82	+ m = min - this->constants[i] * maxs[i];
	83	+
	84	+ v = (lb - m) / this->constants[i];
	85	+
	86	+ if (v * this->constants[i] + m == lb)
	87	+ fd_update_domain_and_check(del_val, v, VAR(this, i));
	88	+
	89	+ fd__constraint_set_entailed(this);
	90	+ }
	91	+
	92	+#ifdef CONSTRAINT_TEMPS
	93	+ // save values
	94	+ *base = lb;
	95	+ *(base + 1) = ub;
	96	+ *(base + 2) = min;
	97	+ *(base + 3) = max;
	98	+ *(base + 4) = unset;
	99	+
	100	+ fd__constraint_remember(this);
	101	+#endif
	102	+
	103	+ return FD_OK;
	104	+}
	105	+
	106	+static int fd_poly_ne_propagate2(fd_constraint this, fd_int culprit)
	107	+{
	108	+#ifdef CONSTRAINT_TEMPS
	109	+ int ub, lb;
	110	+ int min, max;
	111	+ int terms = this->nconstants;
	112	+ int mins, maxs;
	113	+ int unset;
	114	+ int i;
	115	+ int *base;
	116	+ int x, c;
	117	+ int nmin, nmin_x, nmax, nmax_x;
	118	+
	119	+ if (!fd__constraint_data_valid(this))
	120	+ return fd_poly_ne_filter(this); // ignores culprit
	121	+
	122	+ // bounds filtering
	123	+
	124	+ base = constraint_memory[this->index];
	125	+
	126	+ mins = base + 5;
	127	+ maxs = mins + terms;
	128	+
	129	+ lb = *base;
	130	+ ub = *(base + 1);
	131	+
	132	+ min = *(base + 2);
	133	+ max = *(base + 3);
	134	+
	135	+ unset = *(base + 4);
	136	+
	137	+ if (culprit == VAR(this, this->nvariables - 1))
	138	+ {
	139	+ // the culprit is the sum variable
	140	+ int nlb, nub;
	141	+
	142	+ nlb = _fd_var_min(culprit);
	143	+ nub = _fd_var_max(culprit);
	144	+
	145	+ if (nlb == lb && nub == ub)
	146	+ return FD_OK;
	147	+
	148	+ if (min > nub \|\| max < nlb)
	149	+ fd__constraint_set_entailed(this);
	150	+
	151	+ if (unset == 1 && nlb == nub)
	152	+ {
	153	+ int m, v;
	154	+
	155	+ // find out which is the non-singleton variable
	156	+ for (i = this->nvariables - 2; i >= 0; --i)
	157	+ if (this->constants[i] && mins[i] != maxs[i])
	158	+ break;
	159	+
	160	+ if (this->constants[i] > 0)
	161	+ m = min - this->constants[i] * mins[i];
	162	+ else
	163	+ m = min - this->constants[i] * maxs[i];
	164	+
	165	+ v = (nlb - m) / this->constants[i];
	166	+
	167	+ if (v * this->constants[i] + m == nlb)
	168	+ fd_update_domain_and_check(del_val, v, VAR(this, i));
	169	+
	170	+ fd__constraint_set_entailed(this);
	171	+ }
	172	+
	173	+ *base = nlb;
	174	+ *(base + 1) = nub;
	175	+
	176	+ return FD_OK;
	177	+ }
	178	+
	179	+ // the culprit appears in one of the terms, find out which one(s)
	180	+ for (x = 0; culprit != VAR(this, x); ++x)
	181	+ ;
	182	+
	183	+ nmin_x = _fd_var_min(VAR(this, x));
	184	+ nmax_x = _fd_var_max(VAR(this, x));
	185	+
	186	+ if (nmin_x == mins[x] && nmax_x == maxs[x])
	187	+ return FD_OK;
	188	+
	189	+ nmin = min;
	190	+ nmax = max;
	191	+
	192	+ do
	193	+ {
	194	+ c = this->constants[x];
	195	+
	196	+ if (c && nmin_x == nmax_x)
	197	+ unset--;
	198	+
	199	+ if (c > 0)
	200	+ {
	201	+ nmin = nmin + (nmin_x - mins[x]) * c;
	202	+ nmax = nmax - (maxs[x] - nmax_x) * c;
	203	+ }
	204	+ else if (c < 0)
	205	+ {
	206	+ nmin = nmin - (maxs[x] - nmax_x) * c;
	207	+ nmax = nmax + (nmin_x - mins[x]) * c;
	208	+ }
	209	+
	210	+ mins[x] = nmin_x;
	211	+ maxs[x] = nmax_x;
	212	+
	213	+ while (++x < terms && culprit != VAR(this, x))
	214	+ ;
	215	+ }
	216	+ while (x < terms);
	217	+
	218	+ if (nmin > ub \|\| nmax < lb)
	219	+ {
	220	+ fd__constraint_set_entailed(this);
	221	+
	222	+ return FD_OK;
	223	+ }
	224	+
	225	+ if (nmin == nmax)
	226	+ {
	227	+ fd_update_domain_and_check(del_val, nmin, VAR(this, this->nvariables - 1));
	228	+
	229	+ fd__constraint_set_entailed(this);
	230	+ }
	231	+ else if (unset == 1 && lb == ub)
	232	+ {
	233	+ int m, v;
	234	+
	235	+ // find out which is the non-singleton variable
	236	+ for (i = terms - 1; i >= 0; --i)
	237	+ if (this->constants[i] && mins[i] != maxs[i])
	238	+ break;
	239	+
	240	+ if (this->constants[i] > 0)
	241	+ m = min - this->constants[i] * mins[i];
	242	+ else
	243	+ m = min - this->constants[i] * maxs[i];
	244	+
	245	+ v = (lb - m) / this->constants[i];
	246	+
	247	+ if (v * this->constants[i] + m == lb)
	248	+ fd_update_domain_and_check(del_val, v, VAR(this, i));
	249	+
	250	+ fd__constraint_set_entailed(this);
	251	+ }
	252	+
	253	+ *(base + 2) = nmin;
	254	+ *(base + 3) = nmax;
	255	+ *(base + 4) = unset;
	256	+
	257	+ return FD_OK;
	258	+#else /* CONSTRAINT_TEMPS */
	259	+ return fd_poly_ne_filter(this); // ignores culprit
	260	+#endif /* CONSTRAINT_TEMPS */
	261	+}
	262	+
	263	+fd_constraint fd_poly_ne(int cs[], fd_int xs[], int nterms, fd_int y)
	264	+{
	265	+ fd_constraint c = _fd_constraint_new(nterms + 1, nterms);
	266	+ int i;
	267	+
	268	+ if (c)
	269	+ {
	270	+ for (i = 0; i < nterms; ++i)
	271	+ c->variables[i] = FD_INT2C_VAR(xs[i]);
	272	+ c->variables[nterms] = FD_INT2C_VAR(y);
	273	+ for (i = 0; i < nterms; ++i)
	274	+ c->constants[i] = cs[i];
	275	+#ifdef CONSTRAINT_CLASS
	276	+ c->kind = FD_CONSTR_POLY_NE;
	277	+#else /* CONSTRAINT_CLASS */
	278	+ c->propagator2 = fd_poly_ne_propagate2;
	279	+#endif /* CONSTRAINT_CLASS */
	280	+
	281	+ for (i = 0; i < c->nvariables; ++i)
	282	+ _fd_var_add_constraint(VAR(c, i), c);
	283	+
	284	+ _fd_add_constraint(c);
	285	+ }
	286	+
	287	+ return c;
	288	+}
...	...