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