/* poly-eq(C, X, y) == sum(C . X) = y */

// caveat: may not work correctly when there are more than one
// occurrences of some variable

static int fd_poly_eq_filter(fd_constraint this)
{
  int ub, lb;
  int min, max;
  int terms = this->nconstants;
  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 + 4) * sizeof(int));

  base = constraint_memory[this->index];

  mins = base + 4;
  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;

      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 > ub || max < lb)
    return FD_NOSOLUTION;

  if (min > lb)

    {

      fd_update_domain_and_check(del_lt, min, VAR(this, this->nvariables - 1));

      lb = min;
    }

  if (max < ub)
    {
      fd_update_domain_and_check(del_gt, max, VAR(this, this->nvariables - 1));

      ub = max;

    }

  if (min == max)

    {
      fd__constraint_set_entailed(this);
    }
  else
    {
      // bounds domain filtering
      // XXX: number of iterations depends on the size of the domains

      int omin, omax;

      do {
	omin = min;
	omax = max;

	if (max > ub)

	  {

	    for (i = 0; i < terms; ++i)

	      {

		int c = this->constants[i];
		int imin = mins[i];
		int imax = maxs[i];

		if (c == 0 || imin == imax)  continue;

		if (c > 0)
		  {
		    // enforce min(sum{j!=i}(c[j] * x[j])) + c[i] * max[i] <= max[y]

		    if (min + (imax - imin) * c > ub)
		      {
			// imax = floor((ub - min) / c) + imin
			imax = (ub - min) / c + imin;

			fd_update_domain_and_check(del_gt, imax, VAR(this, i));

			max -= (maxs[i] - imax) * c;

			maxs[i] = imax;
		      }
		  }
		else
		  {
		    // enforce min(sum{j!=i}(c[j] * x[j])) + c[i] * min[i] <= max[y]

		    if (min + (imin - imax) * c > ub)
		      {
			// imin = ceil((ub - min) / c) + imax
			imin = (ub - min) / c + imax;

			fd_update_domain_and_check(del_lt, imin, VAR(this, i));

			max -= (mins[i] - imin) * c;

			mins[i] = imin;
		      }
		  }

	      }

	    if (max < lb)
	      return FD_NOSOLUTION;

	  }

	if (min < lb)
	  {
	    for (i = 0; i < terms; ++i)

	      {

		int c = this->constants[i];
		int imin = mins[i];
		int imax = maxs[i];

		if (c == 0 || imin == imax)  continue;

		if (c > 0)

		  {

		    // enforce max(sum{j!=i}(c[j] * x[j])) + c[i] * min[i] >= min[y]

		    if (max + (imin - imax) * c < lb)
		      {
			// imin = ceil((lb - max) / c) + imax
			imin = (lb - max) / c + imax;

			fd_update_domain_and_check(del_lt, imin, VAR(this, i));

			min += (imin - mins[i]) * c;

			mins[i] = imin;
		      }

		  }

		else

		  {

		    // enforce max(sum{j!=i}(c[j] * x[j])) + c[i] * max[i] >= min[y]

		    if (max + (imax - imin) * c < lb)
		      {
			// imax = floor((lb - max) / c) + imin
			imax = (lb - max) / c + imin;

			fd_update_domain_and_check(del_gt, imax, VAR(this, i));

			min += (imax - maxs[i]) * c;

			maxs[i] = imax;
		      }

		  }
	      }

	    if (min > ub)
	      return FD_NOSOLUTION;

	  }

	if (min > lb && _fd_var_del_lt(min, 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));

	    lb = min;
	  }

	if (max < ub && _fd_var_del_gt(max, 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));

	    ub = max;
	  }

      } while (min != max && (min != omin || max != omax));

      if (min == max)
	fd__constraint_set_entailed(this);

    }

#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_poly_eq_propagate2(fd_constraint this, fd_int culprit)
{
#ifdef CONSTRAINT_TEMPS
  int ub, lb;
  int min, max;
  int terms = this->nconstants;
  int *mins, *maxs;
  int i;
  int *base;
  int x, c;
  int nmin, nmin_x, nmax, nmax_x;

  int nlb, nub;

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

  // bounds filtering

  base = constraint_memory[this->index];

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

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

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

  if (culprit == VAR(this, this->nvariables - 1))
    {
      // the culprit is the sum variable

      nlb = _fd_var_min(culprit);
      nub = _fd_var_max(culprit);

      if (nlb == lb && nub == ub)
	return FD_OK;

      if (min > nub || max < nlb)
	return FD_NOSOLUTION;
    }
  else
    {
      // 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;

      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 > ub || nmax < lb)
	    return FD_NOSOLUTION;

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

	  while (++x < terms && culprit != VAR(this, x))
	    ;
	}
      while (x < terms);

      if (nmin == min && nmax == max)
	return FD_OK;

      if (nmin > lb)
	{
	  fd_update_domain_and_check(del_lt, nmin, VAR(this, this->nvariables - 1));

	  nlb = nmin;
	}

      if (nmax < ub)
	{
	  fd_update_domain_and_check(del_gt, nmax, VAR(this, this->nvariables - 1));

	  nub = nmax;
	}
    }

  // (nmin, nlb, nub, nmax) != (min, lb, ub, max)

  if (nmin == nmax)
    {
      fd__constraint_set_entailed(this);
    }
  else
    {
      // bounds domain propagation
      // XXX: number of iterations depends on the size of the domains

      int onmin, onmax;

      do {
	onmin = nmin;
	onmax = nmax;

	if (nmin < nlb)
	  {

	    for (i = 0; i < terms; ++i)
	      {

		int c = this->constants[i];
		int imin = mins[i];
		int imax = maxs[i];

		if (c == 0 || imin == imax)  continue;

		if (c > 0)

		  {

		    // enforce max(sum{j!=i}(c[j] * x[j])) + c[i] * min[i] >= min[y]

		    if (nmax + (imin - imax) * c < nlb)

		      {

			// imin = ceil((nlb - nmax) / c) + imax
			imin = (nlb - nmax) / c + imax;

			fd_update_domain_and_check(del_lt, imin, VAR(this, i));

			nmin += (imin - mins[i]) * c;

			mins[i] = imin;

		      }
		  }

		else

		  {

		    // enforce max(sum{j!=i}(c[j] * x[j])) + c[i] * max[i] >= min[y]

		    if (nmax + (imax - imin) * c < nlb)

		      {

			// imax = floor((nlb - nmax) / c) + imin
			imax = (nlb - nmax) / c + imin;

			fd_update_domain_and_check(del_gt, imax, VAR(this, i));

			nmin += (imax - maxs[i]) * c;

			maxs[i] = imax;

		      }
		  }
	      }

	    if (nmin > nub)
	      return FD_NOSOLUTION;
	  }

	if (nmax > nub)
	  {
	    for (i = 0; i < terms; ++i)

	      {

		int c = this->constants[i];
		int imin = mins[i];
		int imax = maxs[i];

		if (c == 0 || imin == imax)  continue;

		if (c > 0)

		  {

		    // enforce min(sum{j!=i}(c[j] * x[j])) + c[i] * max[i] <= max[y]

		    if (nmin + (imax - imin) * c > nub)

		      {

			// imax = floor((nub - nmin) / c) + imin
			imax = (nub - nmin) / c + imin;

			fd_update_domain_and_check(del_gt, imax, VAR(this, i));

			nmax -= (maxs[i] - imax) * c;

			maxs[i] = imax;

		      }
		  }

		else

		  {

		    // enforce min(sum{j!=i}(c[j] * x[j])) + c[i] * min[i] <= max[y]

		    if (nmin + (imin - imax) * c > nub)

		      {

			// imin = ceil((nub - nmin) / c) + imax
			imin = (nub - nmin) / c + imax;

			fd_update_domain_and_check(del_lt, imin, VAR(this, i));

			nmax -= (mins[i] - imin) * c;

			mins[i] = imin;

		      }
		  }
	      }

	    if (nmax < nlb)
	      return FD_NOSOLUTION;
	  }

	if (nmin > nlb && _fd_var_del_lt(nmin, 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));

	    nlb = nmin;
	  }

	if (nmax < nub && _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));

	    nub = nmax;
	  }
      } while (nmin != nmax && (nmin != onmin || nmax != onmax));

      if (nmin == nmax)
	fd__constraint_set_entailed(this);

    }

  *(base + 0) = nlb;
  *(base + 1) = nub;

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

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

fd_constraint fd_poly_eq(int cs[], fd_int xs[], int nterms, fd_int y)
{

  fd_constraint fd_poly_eq_k(int[], fd_int[], int, int);
  fd_constraint c;
  int i, k;

  // specialise when y is a singleton
  if (fd_var_single(y, &k))
    return fd_poly_eq_k(cs, xs, nterms, k);

  c = fd__constraint_new(nterms + 1, nterms);

  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];

      c->kind = FD_CONSTR_POLY_EQ;

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

      _fd_add_constraint(c);
    }

  return c;
}