/* 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;