/* knapsack2(X, Y, Z) == min(Z) <= sum(X . Y) <= max(Z) */

// XXX: assumes all values are non-negative

static int fd_knapsack2_filter(fd_constraint this)
{
  int ub, lb;
  int min, max;
  int terms = (this->nvariables - 1) / 2;
  int *mins, *maxs;	// XXX: trouble possible if a variable appears repeated
  int i;
#ifdef CONSTRAINT_TEMPS
  int *base;

  assert(!fd__constraint_data_valid(this));

  if (!constraint_memory[this->index])
    constraint_memory[this->index] = malloc((2 * 2 * terms + 4) * sizeof(int));

  base = constraint_memory[this->index];

  mins = base + 4;
  maxs = mins + 2 * terms;
#else
  mins = alloca(2 * terms * sizeof(*mins));
  maxs = alloca(2 * 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 of the terms
  min = 0;

  for (i = 0; i < terms; ++i)
    {
      mins[i] = _fd_var_min(VAR(this, i));
      mins[terms + i] = _fd_var_min(VAR(this, terms + i));

      min += mins[i] * mins[terms + i];

      if (min > ub)
	return FD_NOSOLUTION;
    }

  // sum the maxima of the terms
  max = 0;

  for (i = 0; i < terms; ++i)
    {
      maxs[i] = _fd_var_max(VAR(this, i));
      maxs[terms + i] = _fd_var_max(VAR(this, terms + i));

      max += maxs[i] * maxs[terms + i];
    }

  if (max < lb)
    return FD_NOSOLUTION;

  // XXX: poor man's propagation
  if (min == ub)
    for (i = 0; i < terms; ++i)
      {
	int min_x = mins[i];
	int min_y = mins[terms + i];

	if (min_x != 0 || min_y != 0)
	  {
	    if (min_y != 0 && _fd_var_del_gt(min_x, VAR(this, i)))
	      _fd_revise_connected(this, VAR(this, i));

	    if (min_x != 0 && _fd_var_del_gt(min_y, VAR(this, terms + i)))
	      _fd_revise_connected(this, VAR(this, terms + i));
	  }
      }
  else if (max == lb)
    for (i = 0; i < terms; ++i)
      {
	int max_x = maxs[i];
	int max_y = maxs[terms + i];

	if (max_x != 0 && max_y != 0)
	  {
	    if (_fd_var_del_lt(max_x, VAR(this, i)))
	      _fd_revise_connected(this, VAR(this, i));

	    if (_fd_var_del_lt(max_y, VAR(this, terms + i)))
	      _fd_revise_connected(this, VAR(this, terms + i));
	  }
      }
  else if (max > ub)
    for (i = 0; i < terms; ++i)
      {
	int min_x = mins[i];
	int max_x = maxs[i];
	int min_y = mins[terms + i];
	int max_y = maxs[terms + i];

	if (min_x * (max_y - min_y) > ub - min)
	  {
	    max_y = (ub - min + min_x * min_y) / min_x;

	    _fd_var_del_gt(max_y, VAR(this, terms + i));

	    _fd_revise_connected(this, VAR(this, terms + i));
	  }

	if ((max_x - min_x) * min_y > ub - min)
	  {
	    max_x = (ub - min + min_x * min_y) / min_y;

	    _fd_var_del_gt(max_x, VAR(this, i));

	    _fd_revise_connected(this, VAR(this, i));
	  }
      }

  if ((min > lb && _fd_var_del_lt(min, VAR(this, this->nvariables - 1))) |
      (max < ub && _fd_var_del_gt(max, VAR(this, this->nvariables - 1))))
    _fd_revise_connected(this, VAR(this, this->nvariables - 1));

#ifdef CONSTRAINT_TEMPS
  // save values
  *base = lb;
  *(base + 1) = ub;
  *(base + 2) = min;
  *(base + 3) = max;

  fd__constraint_remember(this);
#endif

  return FD_OK;
}

static int fd_knapsack2_propagate2(fd_constraint this, fd_int culprit)
{
#ifdef CONSTRAINT_TEMPS
  int ub, lb;
  int min, max;
  int terms = (this->nvariables - 1) / 2;
  int *mins, *maxs;	// XXX: trouble possible if a variable appears repeated
  int i;
  int *base;
  int x, y;
  int nmin, nmin_x, nmax, nmax_x;

  if (!fd__constraint_data_valid(this))
    return fd_knapsack2_filter(this);	// ignores culprit

  // bounds filtering

  base = constraint_memory[this->index];

  mins = base + 4;
  maxs = mins + 2 * terms;

  lb = *base;
  ub = *(base + 1);

  min = *(base + 2);
  max = *(base + 3);

  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 (nlb != lb)
	{
	  if (max < nlb)
	    return FD_NOSOLUTION;

	  if (max == nlb)
	    for (i = 0; i < terms; ++i)
	      {
		int max_x = maxs[i];
		int max_y = maxs[terms + i];

		if (max_x != 0 && max_y != 0)
		  {
		    if (_fd_var_del_lt(max_x, VAR(this, i)))
		      {
			if (fd_domain_empty(VAR(this, i)))
			  return FD_NOSOLUTION;

			_fd_revise_connected(this, VAR(this, i));
		      }

		    if (_fd_var_del_lt(max_y, VAR(this, terms + i)))
		      {
			if (fd_domain_empty(VAR(this, terms + i)))
			  return FD_NOSOLUTION;

			_fd_revise_connected(this, VAR(this, terms + i));
		      }
		  }
	      }

	  *base = nlb;
	}

      if (nub != ub)
	{
	  if (min > nub)
	    return FD_NOSOLUTION;

	  if (min == nub)
	    for (i = 0; i < terms; ++i)
	      {
		int min_x = mins[i];
		int min_y = mins[terms + i];

		if (min_x != 0 || min_y != 0)
		  {
		    if (min_y != 0 && _fd_var_del_gt(min_x, VAR(this, i)))
		      {
			if (fd_domain_empty(VAR(this, i)))
			  return FD_NOSOLUTION;

			_fd_revise_connected(this, VAR(this, i));
		      }

		    if (min_x != 0 && _fd_var_del_gt(min_y, VAR(this, terms + i)))
		      {
			if (fd_domain_empty(VAR(this, terms + i)))
			  return FD_NOSOLUTION;

			_fd_revise_connected(this, VAR(this, terms + i));
		      }
		  }
	      }

	  *(base + 1) = nub;
	}

      return FD_OK;
    }

  // the culprit appears in one of the terms, find out which one
  for (x = 0; culprit != VAR(this, x); ++x)
    ;

  y = (x < terms) ? x + terms : x - terms; // the other variable in the term

  nmin_x = _fd_var_min(VAR(this, x));
  nmin = min + (nmin_x - mins[x]) * mins[y];

  nmax_x = _fd_var_max(VAR(this, x));
  nmax = max - (maxs[x] - nmax_x) * maxs[y];

  if (nmin > ub || nmax < lb)
    return FD_NOSOLUTION;

  if ((nmin > lb && _fd_var_del_lt(nmin, VAR(this, this->nvariables - 1))) |
      (nmax < ub && _fd_var_del_gt(nmax, VAR(this, this->nvariables - 1))))
    {
      if (fd_domain_empty(VAR(this, this->nvariables - 1)))
	return FD_NOSOLUTION;

      _fd_revise_connected(this, VAR(this, this->nvariables - 1));
    }

  mins[x] = nmin_x;
  maxs[x] = nmax_x;

  *(base + 2) = nmin;
  *(base + 3) = nmax;

  return FD_OK;
#else /* CONSTRAINT_TEMPS */
  return fd_knapsack2_filter(this);	// ignores culprit
#endif /* CONSTRAINT_TEMPS */
}

fd_constraint fd_knapsack2(fd_int xs[], fd_int ys[], int nvariables, fd_int limits)
{
  fd_constraint c = fd__constraint_new(2 * nvariables + 1, 0);
  int i;

  if (c)
    {
      for (i = 0; i < nvariables; ++i)
	c->variables[i] = FD_INT2C_VAR(xs[i]);
      for (; i < 2 * nvariables; ++i)
	c->variables[i] = FD_INT2C_VAR(ys[i - nvariables]);
      c->variables[2 * nvariables] = FD_INT2C_VAR(limits);

      c->kind = FD_CONSTR_KNAPSACK2;

      for (i = 0; i < c->nvariables; ++i)
	_fd_var_add_constraint(VAR(c, i), c);

      _fd_add_constraint(c);
    }

  return c;
}