/* sum_prod */
/* sum_prod(X, Y, k) == X . Y = k */

// XXX: assumes all values are non-negative

static int fd_sum_prod_filter(fd_constraint this)
{
#ifdef CONSTRAINT_TEMPS
  int k;
  int min, max;
  int terms = this->nvariables / 2;
  int *mins, *maxs;
  int *base;
  int i;

  assert(!fd__constraint_data_valid(this));

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

  base = constraint_memory[this->index];

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

  k = this->constants[0];

  // 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 > k)
	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 < k)
    return FD_NOSOLUTION;

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

	  if (ymin != 0 && maxs[i] > xmin && _fd_var_del_gt(xmin, VAR(this, i)))
	    _fd_revise_connected(this, VAR(this, i));

	  if (xmin != 0 && maxs[terms + i] > ymin &&
	      _fd_var_del_gt(ymin, VAR(this, terms + i)))
	    _fd_revise_connected(this, VAR(this, terms + i));
	}

      fd__constraint_set_entailed(this);
    }
  else if (max == k)
    {
      for (i = 0; i < terms; ++i)
	{
	  int xmax = maxs[i];
	  int ymax = maxs[terms + i];

	  if (xmax != 0 && ymax != 0)
	    {
	      if (mins[i] < xmax && _fd_var_del_lt(xmax, VAR(this, i)))
		_fd_revise_connected(this, VAR(this, i));

	      if (mins[terms + i] < ymax &&
		  _fd_var_del_lt(ymax, VAR(this, terms + i)))
		_fd_revise_connected(this, VAR(this, terms + i));
	    }
	}

      fd__constraint_set_entailed(this);
    }
  else	// min < k < max
    {
      for (i = 0; i < terms; ++i)
	{
	  int xmin = mins[i];
	  int xmax = maxs[i];
	  int ymin = mins[terms + i];
	  int ymax = maxs[terms + i];
	  int changed_x = 0, changed_y = 0;

	  if (min + (xmax - xmin) * ymin > k)
	    {
	      int nxmax = (k - min) / ymin + xmin;

	      changed_x = _fd_var_del_gt(nxmax, VAR(this, i));
	    }

	  if (min + xmin * (ymax - ymin) > k)
	    {
	      int nymax = (k - min) / xmin + ymin;

	      changed_y = _fd_var_del_gt(nymax, VAR(this, terms + i));
	    }

	  if (max - (xmax - xmin) * ymax < k)
	    {
	      int nxmin = (k - max) / ymax + xmax;

	      changed_x |= _fd_var_del_lt(nxmin, VAR(this, i));
	    }

	  if (max - xmax * (ymax - ymin) < k)
	    {
	      int nymin = (k - max) / xmax + ymax;

	      changed_y |= _fd_var_del_lt(nymin, VAR(this, terms + i));
	    }

	  if (changed_x)
	    assert(!fd_domain_empty(VAR(this, i))),
	    _fd_revise_connected(NULL, VAR(this, i));

	  if (changed_y)
	    assert(!fd_domain_empty(VAR(this, terms + i))),
	    _fd_revise_connected(NULL, VAR(this, terms + i));
	}
    }

  // save values
  *base = min;
  *(base + 1) = max;

  fd__constraint_remember(this);

  return FD_OK;
#else /* CONSTRAINT_TEMPS */
  int k;
  int min, max;
  int terms = this->nvariables / 2;
  int i;

  k = this->constants[0];

  // sum the minima of the terms
  min = 0;

  for (i = 0; i < terms; ++i)
    {
      min += _fd_var_min(VAR(this, i)) * _fd_var_min(VAR(this, terms + i));

      if (min > k)
	return FD_NOSOLUTION;
    }

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

  for (i = 0; i < terms; ++i)
    max += _fd_var_max(VAR(this, i)) * _fd_var_max(VAR(this, terms + i));

  if (max < k)
    return FD_NOSOLUTION;

  // XXX: poor man's propagation
  if (min == k)
    for (i = 0; i < terms; ++i)
      {
	int xmin = _fd_var_min(VAR(this, i));
	int ymin = _fd_var_min(VAR(this, terms + i));

	if (ymin != 0 && _fd_var_del_gt(xmin, VAR(this, i)))
	  _fd_revise_connected(this, VAR(this, i));

	if (xmin != 0 && _fd_var_del_gt(ymin, VAR(this, terms + i)))
	  _fd_revise_connected(this, VAR(this, terms + i));
      }
  else if (max == k)
    for (i = 0; i < terms; ++i)
      {
	int xmax = _fd_var_max(VAR(this, i));
	int ymax = _fd_var_max(VAR(this, terms + i));

	if (xmax != 0 && ymax != 0)
	  {
	    if (_fd_var_del_lt(xmax, VAR(this, i)))
	      _fd_revise_connected(this, VAR(this, i));

	    if (_fd_var_del_lt(ymax, VAR(this, terms + i)))
	      _fd_revise_connected(this, VAR(this, terms + i));
	  }
      }

  return FD_OK;
#endif /* CONSTRAINT_TEMPS */
}

static int fd_sum_prod_propagate2(fd_constraint this, fd_int culprit)
{
#ifdef CONSTRAINT_TEMPS
  int k;
  int min, max;
  int terms = this->nvariables / 2;
  int *mins, *maxs;
  int *base;
  int nmin, nmax, nmin_v, nmax_v;
  int i;

  if (!fd__constraint_data_valid(this))
    return fd_sum_prod_filter(this);

  // fetch values
  base = constraint_memory[this->index];

  min = *base;
  max = *(base + 1);
  mins = base + 2;
  maxs = base + 2 + 2 * terms;

  k = this->constants[0];

  // get the new bounds of culprit's domain
  nmin_v = _fd_var_min(culprit);
  nmax_v = _fd_var_max(culprit);

  // find out which term(s) culprit belongs to
  for (i = 0; culprit != VAR(this, i); ++i)
    ;

  if (nmin_v == mins[i] && nmax_v == maxs[i])
    return FD_OK;	// nothing has (meaningfully) changed

  nmin = min;
  nmax = max;

  do
    {
      int j = (i + terms) % this->nvariables;

      nmin += (nmin_v - mins[i]) * mins[j];
      nmax -= (maxs[i] - nmax_v) * maxs[j];

      mins[i] = nmin_v;
      maxs[i] = nmax_v;

      // find out which other terms culprit belongs to
      for (++i; i < this->nvariables && culprit != VAR(this, i); ++i)
	;
    }
  while (i < this->nvariables);

  if (nmin > k || nmax < k)
    return FD_NOSOLUTION;

  if (nmin == min && nmax == max)
    return FD_OK;	// nothing has (meaningfully) changed

  // XXX: poor man's propagation
  if (nmin == k)
    {
      for (i = 0; i < terms; ++i)
	{
	  int xmin = mins[i];
	  int ymin = mins[terms + i];

	  if (ymin != 0 && maxs[i] > xmin && _fd_var_del_gt(xmin, VAR(this, i)))
	    {
	      if (fd_domain_empty(VAR(this, i)))
		return FD_NOSOLUTION;

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

	  if (xmin != 0 && maxs[terms + i] > ymin &&
	      _fd_var_del_gt(ymin, VAR(this, terms + i)))
	    {
	      if (fd_domain_empty(VAR(this, terms + i)))
		return FD_NOSOLUTION;

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

      fd__constraint_set_entailed(this);
    }
  else if (nmax == k)
    {
      for (i = 0; i < terms; ++i)
	{
	  int xmax = maxs[i];
	  int ymax = maxs[terms + i];

	  if (xmax != 0 && ymax != 0)
	    {
	      if (mins[i] < xmax && _fd_var_del_lt(xmax, VAR(this, i)))
		{
		  if (fd_domain_empty(VAR(this, i)))
		    return FD_NOSOLUTION;

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

	      if (mins[terms + i] < ymax &&
		  _fd_var_del_lt(ymax, VAR(this, terms + i)))
		{
		  if (fd_domain_empty(VAR(this, terms + i)))
		    return FD_NOSOLUTION;

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

      fd__constraint_set_entailed(this);
    }
  else	// nmin < k < nmax
    {
      for (i = 0; i < terms; ++i)
	{
	  int xmin = mins[i];
	  int xmax = maxs[i];
	  int ymin = mins[terms + i];
	  int ymax = maxs[terms + i];
	  int changed_x = 0, changed_y = 0;

	  if (nmin + (xmax - xmin) * ymin > k)
	    {
	      int nxmax = (k - nmin) / ymin + xmin;

	      changed_x = _fd_var_del_gt(nxmax, VAR(this, i));
	    }

	  if (nmin + xmin * (ymax - ymin) > k)
	    {
	      int nymax = (k - nmin) / xmin + ymin;

	      changed_y = _fd_var_del_gt(nymax, VAR(this, terms + i));
	    }

	  if (nmax - (xmax - xmin) * ymax < k)
	    {
	      int nxmin = (k - nmax) / ymax + xmax;

	      changed_x |= _fd_var_del_lt(nxmin, VAR(this, i));
	    }

	  if (nmax - xmax * (ymax - ymin) < k)
	    {
	      int nymin = (k - nmax) / xmax + ymax;

	      changed_y |= _fd_var_del_lt(nymin, VAR(this, terms + i));
	    }

	  if (changed_x)
	    {
	      if (fd_domain_empty(VAR(this, i)))
		return FD_NOSOLUTION;

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

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

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

  // save new values
  *base = nmin;
  *(base + 1) = nmax;

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

fd_constraint fd_sum_prod(fd_int *xs, fd_int *ys, int nvariables, int k)
{
  fd_constraint c = _fd_constraint_new(2 * nvariables, 1);
  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->constants[0] = k;
#ifdef CONSTRAINT_CLASS
      c->kind = FD_CONSTR_SUM_PROD;
#else /* CONSTRAINT_CLASS */
      c->propagator2 = fd_sum_prod_propagate2;
#endif /* CONSTRAINT_CLASS */

      for (i = 0; i < 2 * nvariables; ++i)
	_fd_var_add_constraint(VAR(c, i), c);

      _fd_add_constraint(c);
    }

  return c;
}