#include <string.h>
#include <alloca.h>
#include <stdlib.h>

#include <assert.h>

#include "fdc_int.h"
#include "variables.h"
#include "values.h"
#include "matching.h"

#define UNSEEN (-2)		// value not seen in any domain

#define WHITE  (-1)
#define GREY   0
#define BLACK  1

#define VARIABLE 0
#define VALUE    1

/*
  Find out which edges, not belonging to the matching, are in
  M-alternating cycles [HCP].

  The cycles correspond to the strongly connected components of the
  value graph, viewed as a directed graph with the matching edges
  going from left (the variables) to right (the values), and the other
  edges going in the opposite direction [HCP].
*/
static void _fd_mark_cycles(fd_constraint c, int8_t *vgraph, int r_verts,
			    int *l_edge, int *r_edge, int *l_scc, int *r_scc,
			    int min)
{
  int nvariables = c->nvariables;
  int8_t *l_colour, *r_colour;
  int *stack, top;
  int8_t *type;			// type of node in the stack: VARIABLE or VALUE
#define push(x, t) (top++, type[top] = (t), stack[top] = (x))
#define pop() stack[top--]
#define top() stack[top]
#define last_branch stack[top + 1]	// the last branch taken when
					// leaving the top() node
  int *l_finish, *r_finish, time;	// finishing times of the DFS algorithm
  int *variables;			// variables ordered by finishing time
  int scc;
  int i, j, v;

  /* Strongly connected components found through the [...] algorithm
     [Cormen]. */

  l_colour = alloca(nvariables * sizeof(*l_colour));
  r_colour = alloca(r_verts * sizeof(*r_colour));
  memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
  memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

  stack = alloca((nvariables + r_verts) * sizeof(*stack));
  type = alloca((nvariables + r_verts) * sizeof(*type));

  variables = alloca(nvariables * sizeof(*variables));

  // forward DFS

  top = -1;

  l_finish = alloca(nvariables * sizeof(*l_finish));
  r_finish = alloca(r_verts * sizeof(*r_finish));
  memset(l_finish, -1, nvariables * sizeof(*l_finish));
  memset(r_finish, -1, r_verts * sizeof(*r_finish));

  time = -1;

  for (i = 0, j = nvariables; i < nvariables; ++i)
    {
      if (l_colour[i] != WHITE)		// already visited this variable
	continue;

      push(i, VARIABLE);
      l_colour[i] = GREY;

      while (top >= 0)
	{
	  switch (type[top])
	    {
	    case VARIABLE:
	      v = l_edge[top()]; // matching edges go from left to right

	      if (r_colour[v] != WHITE)
		{
		  // matched value already visited

		  //l_colour[top()] = BLACK;
		  l_finish[top()] = ++time;
		  variables[--j] = pop();

		  break;
		}

	      push(v, VALUE);
	      last_branch = -1;
	      r_colour[v] = GREY;

	      break;
	    case VALUE:
	      for (v = last_branch + 1; v < nvariables; ++v)
		if (vgraph[v * r_verts + top()] == 1 &&
		    l_colour[v] == WHITE)
		  break;

	      if (v == nvariables)	// no unvisited successor found
		{
		  //r_colour[top()] = BLACK;
		  r_finish[top()] = ++time;
		  pop();

		  break;
		}

	      push(v, VARIABLE);	// sets last_branch to v
	      l_colour[v] = GREY;

	      break;
	    }
	}
    }

  // backward DFS

  top = -1;

  int *l_parent, *r_parent;		// depth-first forest

  l_parent = alloca(nvariables * sizeof(*l_parent));
  r_parent = alloca(r_verts * sizeof(*r_parent));
  memset(l_parent, -1, nvariables * sizeof(*l_parent));
  memset(r_parent, -1, r_verts * sizeof(*r_parent));

  memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
  memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

  for (j = 0; j < nvariables; ++j)
    {
      i = variables[j];

      if (l_colour[i] != WHITE)		// already visited this variable
	continue;

      push(i, VARIABLE);
      last_branch = -1;
      l_colour[i] = GREY;

      while (top >= 0)
	{
	  switch (type[top])
	    {
	    case VARIABLE:
	      for (v = last_branch + 1; v < r_verts; ++v)
		if (vgraph[top() * r_verts + v] == 1 &&
		    r_colour[v] == WHITE &&
		    r_finish[v] > -1 && r_finish[v] < l_finish[i])
		  break;

	      if (v == r_verts)		// no unvisited successor found
		{
		  //l_colour[top()] = BLACK;
		  pop();

		  break;
		}

	      r_parent[v] = top();

	      push(v, VALUE);		// sets last_branch to v
	      r_colour[v] = GREY;

	      break;
	    case VALUE:
	      v = r_edge[top()]; // matching edges go from right to left

	      if (l_colour[v] != WHITE)
		{
		  // matched variable already visited

		  //r_colour[top()] = BLACK;
		  pop();

		  break;
		}

	      l_parent[v] = top();

	      push(v, VARIABLE);
	      last_branch = -1;
	      l_colour[v] = GREY;

	      break;
	    }
	}
    }

  // identify strongly connected components nodes

  int8_t *is_root;
  int *l_children, *l_nc, *r_child;

  //is_root = alloca(nvariables * sizeof(*is_root));
  is_root = l_colour;					// XXX: reuse memory
  l_children = alloca(nvariables * r_verts * sizeof(*l_children));
  //l_nc = alloca(nvariables * sizeof(*l_nc));
  l_nc = stack;						// XXX: reuse memory
  r_child = alloca(r_verts * sizeof(*r_child));

  memset(is_root, 1, nvariables * sizeof(*is_root));
  memset(l_nc, 0, nvariables * sizeof(*l_nc));
  memset(r_child, -1, r_verts * sizeof(*r_child));

  /* identify strongly connected components nodes */

  // reverse edges of the depth-first search forest
  for (i = 0; i < nvariables; ++i)
    {
      if (l_parent[i] == -1)
	  continue;

      is_root[i] = 0;

      r_child[l_parent[i]] = i;
    }

  for (v = 0; v < r_verts; ++v)
    {
      if (r_parent[v] == -1)
	continue;

      i = r_parent[v];

      l_children[i * r_verts + l_nc[i]] = v;
      l_nc[i]++;
    }

  // mark values as used
  int *l_queue, l_head, l_tail, *r_queue, r_head, r_tail;

  //l_queue = alloca(nvariables * sizeof(*l_queue));
  l_queue = l_finish;					// XXX: reuse memory
  //r_queue = alloca(r_verts * sizeof(*r_queue));
  r_queue = r_finish;					// XXX: reuse memory

  memset(l_scc, -1, nvariables * sizeof(*l_scc));
  memset(r_scc, -1, r_verts * sizeof(*r_scc));
  scc = -1;

  for (i = 0; i < nvariables; ++i)
    {
      if (!is_root[i] || l_nc[i] == 0)
	continue;

      scc++;

      l_head = l_tail = 0;
      l_queue[l_tail] = i;
      r_head = 0; r_tail = -1;

      while (l_head <= l_tail)
	{
	  do
	    {
	      int n = l_queue[l_head++];

	      l_scc[n] = scc;	// record the scc the variable belongs to

	      for (v = 0; v < l_nc[n]; ++v)
		r_queue[++r_tail] = l_children[n * r_verts + v];
	    }
	  while (l_head <= l_tail);

	  while (r_head <= r_tail)
	    {
	      int n = r_queue[r_head++];

	      r_scc[n] = scc;	// record the scc the value belongs to

	      if (r_child[n] != -1)
		l_queue[++l_tail] = r_child[n];
	    }
	}

      // all the nodes from one strongly connected component are in
      // l_queue (the variables) and r_queue (the values)

      for (j = 0; j < l_head; ++j)
	{
	  int var = l_queue[j];

	  for (v = 0; v < r_head; ++v)
	    {
	      int val = r_queue[v];

	      if (vgraph[var * r_verts + val] == 1)
		{
//		  _fd_debug("marking %d in %d's domain\n", val + min, VAR(c, var)->index);

		  vgraph[var * r_verts + val] = 3;
	      }
	    }
	}

      assert(l_head == r_head);		// bug catching assert
    }
#undef push
#undef pop
#undef top
#undef last_branch
}

/*
  Recompute the strongly connected components wrt the updated
  matching.
*/
static void _fd_remark_cycles(fd_constraint c, int culprit, int8_t *vgraph,
			      int r_verts, int *l_edge, int *r_edge,
			      int *l_scc, int *r_scc, int min)
{
  int nvariables = c->nvariables;
  int8_t *l_colour, *r_colour;
  int *stack, top;
  int8_t *type;			// type of node in the stack: VARIABLE or VALUE
#define push(x, t) (top++, type[top] = (t), stack[top] = (x))
#define pop() stack[top--]
#define top() stack[top]
#define last_branch stack[top + 1]	// the last branch taken when
					// leaving the top() node
  int *l_finish, *r_finish, time;	// finishing times of the DFS algorithm
  int *variables;			// variables ordered by finishing time
  int i, j, v;

  int scc;				// the scc culprit belonged to
  int max_scc;

  /* Strongly connected components found through the [...] algorithm
     [Cormen]. */

  max_scc = scc = l_scc[culprit];

  for (i = 0; i < nvariables; ++i)
    if (l_scc[i] > max_scc)
      max_scc = l_scc[i];

  for (i = 0; i < r_verts; ++i)
    if (r_scc[i] == scc)
      r_scc[i] = -1;

  l_colour = alloca(nvariables * sizeof(*l_colour));
  r_colour = alloca(r_verts * sizeof(*r_colour));
  memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
  memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

  stack = alloca((nvariables + r_verts) * sizeof(*stack));
  type = alloca((nvariables + r_verts) * sizeof(*type));

  variables = alloca(nvariables * sizeof(*variables));

  // forward DFS

  top = -1;

  l_finish = alloca(nvariables * sizeof(*l_finish));
  r_finish = alloca(r_verts * sizeof(*r_finish));
  memset(l_finish, -1, nvariables * sizeof(*l_finish));
  memset(r_finish, -1, r_verts * sizeof(*r_finish));

  time = -1;

  for (i = 0, j = nvariables; i < nvariables; ++i)
    {
      if (l_scc[i] != scc)		// only need to update the culprit's scc
	continue;

      l_scc[i] = -1;			// a little cleaning up

      if (l_colour[i] != WHITE)		// already visited this variable
	continue;

      push(i, VARIABLE);
      l_colour[i] = GREY;

      while (top >= 0)
	{
	  switch (type[top])
	    {
	    case VARIABLE:
	      v = l_edge[top()]; // matching edges go from left to right

	      if (r_colour[v] != WHITE)
		{
		  // matched value already visited

		  //l_colour[top()] = BLACK;
		  l_finish[top()] = ++time;
		  variables[--j] = pop();

		  break;
		}

	      r_scc[v] = -1;	// value belongs to a new SCC (XXX: check!)

	      push(v, VALUE);
	      last_branch = -1;
	      r_colour[v] = GREY;

	      break;
	    case VALUE:
	      for (v = last_branch + 1; v < nvariables; ++v)
		if (vgraph[v * r_verts + top()] == 1 &&
		    l_colour[v] == WHITE)
		  break;

	      if (v == nvariables)	// no unvisited successor found
		{
		  //r_colour[top()] = BLACK;
		  r_finish[top()] = ++time;
		  pop();

		  break;
		}

	      push(v, VARIABLE);	// sets last_branch to v
	      l_colour[v] = GREY;

	      break;
	    }
	}
    }

  // backward DFS

  top = -1;

  int *l_parent, *r_parent;		// depth-first forest

  l_parent = alloca(nvariables * sizeof(*l_parent));
  r_parent = alloca(r_verts * sizeof(*r_parent));
  memset(l_parent, -1, nvariables * sizeof(*l_parent));
  memset(r_parent, -1, r_verts * sizeof(*r_parent));

  memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
  memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

  for (; j < nvariables; ++j)	// j's value corresponds to the last variable
				// finished above
    {
      i = variables[j];

      if (l_colour[i] != WHITE)		// already visited this variable
	continue;

      push(i, VARIABLE);
      last_branch = -1;
      l_colour[i] = GREY;

      while (top >= 0)
	{
	  switch (type[top])
	    {
	    case VARIABLE:
	      for (v = last_branch + 1; v < r_verts; ++v)
		if (vgraph[top() * r_verts + v] == 1 &&
		    r_colour[v] == WHITE &&
		    r_finish[v] > -1 && r_finish[v] < l_finish[i])
		  break;

	      if (v == r_verts)		// no unvisited successor found
		{
		  //l_colour[top()] = BLACK;
		  pop();

		  break;
		}

	      r_parent[v] = top();

	      push(v, VALUE);		// sets last_branch to v
	      r_colour[v] = GREY;

	      break;
	    case VALUE:
	      v = r_edge[top()]; // matching edges go from right to left

	      if (l_colour[v] != WHITE)
		{
		  // matched variable already visited

		  //r_colour[top()] = BLACK;
		  pop();

		  break;
		}

	      l_parent[v] = top();

	      push(v, VARIABLE);
	      last_branch = -1;
	      l_colour[v] = GREY;

	      break;
	    }
	}
    }

  // identify strongly connected components nodes

  int8_t *is_root;
  int *l_children, *l_nc, *r_child;

  //is_root = alloca(nvariables * sizeof(*is_root));
  is_root = l_colour;					// XXX: reuse memory
  l_children = alloca(nvariables * r_verts * sizeof(*l_children));
  //l_nc = alloca(nvariables * sizeof(*l_nc));
  l_nc = stack;						// XXX: reuse memory
  r_child = alloca(r_verts * sizeof(*r_child));

  memset(is_root, 1, nvariables * sizeof(*is_root));
  memset(l_nc, 0, nvariables * sizeof(*l_nc));
  memset(r_child, -1, r_verts * sizeof(*r_child));

  /* identify strongly connected components nodes */

  // reverse edges of the depth-first search forest
  for (i = 0; i < nvariables; ++i)
    {
      if (l_parent[i] == -1)
	  continue;

      is_root[i] = 0;

      r_child[l_parent[i]] = i;
    }

  for (v = 0; v < r_verts; ++v)
    {
      if (r_parent[v] == -1)
	continue;

      i = r_parent[v];

      l_children[i * r_verts + l_nc[i]] = v;
      l_nc[i]++;
    }

  // mark values as used
  int *l_queue, l_head, l_tail, *r_queue, r_head, r_tail;

  //l_queue = alloca(nvariables * sizeof(*l_queue));
  l_queue = l_finish;					// XXX: reuse memory
  //r_queue = alloca(r_verts * sizeof(*r_queue));
  r_queue = r_finish;					// XXX: reuse memory

  int mscc = max_scc;

  for (i = 0; i < nvariables; ++i)
    {
      if (!is_root[i] || l_nc[i] == 0)
	continue;

      max_scc++;

      l_head = l_tail = 0;
      l_queue[l_tail] = i;
      r_head = 0; r_tail = -1;

      while (l_head <= l_tail)
	{
	  do
	    {
	      int n = l_queue[l_head++];

	      l_scc[n] = max_scc;    // record the scc the variable belongs to

	      for (v = 0; v < l_nc[n]; ++v)
		r_queue[++r_tail] = l_children[n * r_verts + v];
	    }
	  while (l_head <= l_tail);

	  while (r_head <= r_tail)
	    {
	      int n = r_queue[r_head++];

	      r_scc[n] = max_scc;	// record the scc the value belongs to

	      if (r_child[n] != -1)
		l_queue[++l_tail] = r_child[n];
	    }
	}

      // all the nodes from one strongly connected component are in
      // l_queue (the variables) and r_queue (the values)

      for (j = 0; j < l_head; ++j)
	{
	  int var = l_queue[j];

	  for (v = 0; v < r_head; ++v)
	    {
	      int val = r_queue[v];

	      if (vgraph[var * r_verts + val] == 1)
		{
//		  _fd_debug("marking %d in %d's domain\n", val + min, VAR(c, var)->index);

		  vgraph[var * r_verts + val] = 3;
	      }
	    }
	}

      assert(l_head == r_head);		// bug catching assert
    }

  // re-mark edges from unaffected strongly connected components
  for (v = 0; v < r_verts; ++v)
    if (r_scc[v] > -1 && r_scc[v] <= mscc)
      for (i = 0; i < nvariables; ++i)
	if (l_scc[i] == r_scc[v] && vgraph[i * r_verts + v] == 1)
	  {
//	    _fd_debug("marking %d in %d's domain\n", v + min, VAR(c, i)->index);

	    vgraph[i * r_verts + v] = 3;
	  }
#undef push
#undef pop
#undef top
#undef last_branch
}

static void _fd_mark_paths(fd_constraint c, int8_t *vgraph, int r_verts,
			   int *l_edge, int *r_edge, int min)
{
#define vgraph(l, c) vgraph[(l) * r_verts + (c)]
  int nvariables = c->nvariables;
  int *queue, head, tail;	// for the breadth-first search
  int8_t *l_colour, *r_colour;
  int8_t *type;			// type of node in queue: VARIABLE or VALUE
  int i;

  queue = alloca((nvariables + r_verts) * sizeof(*queue));
  type = alloca((nvariables + r_verts) * sizeof(*type));

  l_colour = alloca(nvariables * sizeof(*l_colour));
  r_colour = alloca(r_verts * sizeof(*r_colour));

  memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
  memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

  // find all M-alternating paths starting at an unmatched value
  for (i = 0; i < r_verts; ++i)
    {
      if (r_edge[i] != -1)	// skip values matched or not in any domain
	continue;

      r_colour[i] = GREY;

      queue[head = tail = 0] = i;
      type[head] = VALUE;

      while (head <= tail)
	{
	  switch (type[head])
	    {
	      case VARIABLE:		// node is a variable
		{
		  int v;

		  v = l_edge[queue[head]];

		  if (r_colour[v] == WHITE)
		    {
		      r_colour[v] = GREY;
		      queue[++tail] = v;
		      type[tail] = VALUE;
		    }

		  //l_colour[queue[head]] = BLACK;

		  head++;
		}

		break;
	      case VALUE:		// node is a value (domain element)
		{
		  int v;

		  for (v = 0; v < nvariables; ++v)
		    {
		      if (!vgraph(v, queue[head]) ||
			  vgraph(v, queue[head]) == 2)
			continue;

#if DEBUG_MATCH > 1
		      _fd_debug("marking %d in %d's domain\n",
				queue[head] + min, VAR(c, v)->index);
#endif
		      // this edge is in an M-alternating path
		      vgraph(v, queue[head]) += 4;

		      if (l_colour[v] == WHITE)
			{
			  l_colour[v] = GREY;
			  queue[++tail] = v;
			  type[tail] = VARIABLE;
			}
		    }

		  //r_colour[queue[head]] = BLACK;

		  head++;
		}

		break;
	    }
	}
    }
#undef vgraph
}

#ifndef CONSTRAINT_TEMPS
int _fd_find_matching(fd_constraint c)
#else
int _fd_find_matching(fd_constraint c, int **memory)
#endif
{
  int nvariables = c->nvariables;
  int size;			// constraint memory size
  int8_t *vgraph;		// value graph
  int *l_edge, *r_edge;		// matches (matching edges)
  int *queue, head, tail;	// for the breadth-first search
  int8_t *l_colour, *r_colour;
  int8_t *type;			// type of node in queue: VARIABLE or VALUE
  int *l_parent, *r_parent;	// M-augmenting path
  int min, max, r_verts, nvalues = 0;
  int *pmin;
  int *l_scc, *r_scc;		// strongly connected components
  bool done;
  int i;

  // compute an upper bound on the size of the domains
  min = _fd_var_min(VAR(c, 0));
  max = _fd_var_max(VAR(c, 0));

  for (i = 1; i < nvariables; ++i)
    {
      int v;

      if ((v = _fd_var_min(VAR(c, i))) < min)
	min = v;

      if ((v = _fd_var_max(VAR(c, i))) > max)
	max = v;
    }

  r_verts = max - min + 1;

  // if there are less values in the domains than variables, there is
  // no maximum cardinality matching
  if (r_verts < nvariables)
    return 0;

#ifdef CONSTRAINT_TEMPS
  // memory requirement
  size = 4 * sizeof(int) +			// size, min, r_verts, nvalues
    nvariables * sizeof(*l_edge) + r_verts * sizeof(*r_edge) +
    nvariables * sizeof(*l_scc) + r_verts * sizeof(*r_scc) +
    nvariables * r_verts * sizeof(*vgraph);

  if (*memory == NULL)
    {
#ifdef DEBUG_MATCH
      _fd_debug("all-different memory size is %d\n", size);
#endif

      *memory = malloc(size);		// XXX: NULL

      // save memory block size
      **memory = size;
    }
  else if (**memory < size)
    {
      // the allocated memory is not enough
#ifdef DEBUG_MATCH
      _fd_debug("all-different new memory size is %d\n", size);
#endif

      *memory = realloc(*memory, size);	// XXX: NULL

      // save the new size
      **memory = size;
    }

  assert(sizeof(**memory) == sizeof(size));
  assert(sizeof(**memory) == sizeof(min));
  assert(sizeof(**memory) == sizeof(r_verts));
  assert(sizeof(**memory) == sizeof(nvalues));
  assert(sizeof(**memory) == sizeof(*l_edge));
  assert(sizeof(*l_edge) == sizeof(*r_edge));
  assert(sizeof(*r_edge) == sizeof(*l_scc));
  assert(sizeof(*l_scc) == sizeof(*r_scc));

  pmin = *memory + 1;
  l_edge = pmin + 3;
  r_edge = l_edge + nvariables;
  l_scc = r_edge + r_verts;
  r_scc = l_scc + nvariables;
  vgraph = (int8_t *) (r_scc + r_verts);
#else
  l_edge = alloca(nvariables * sizeof(*l_edge));
  r_edge = alloca(r_verts * sizeof(*r_edge));
  l_scc = alloca(nvariables * sizeof(*l_scc));
  r_scc = alloca(r_verts * sizeof(*r_scc));
  vgraph = alloca(nvariables * r_verts * sizeof(*vgraph));
#endif

  l_parent = alloca(nvariables * sizeof(*l_parent));
  l_colour = alloca(nvariables * sizeof(*l_colour));

  r_parent = alloca(r_verts * sizeof(*r_parent));
  r_colour = alloca(r_verts * sizeof(*r_colour));

  for (i = 0; i < r_verts; ++i)
    r_edge[i] = UNSEEN;

  // build the value graph
  memset(vgraph, 0, nvariables * r_verts * sizeof(*vgraph));

  for (i = 0; i < nvariables; ++i)
    {
      fd_iterator iterator = _fd_val_iterator(DOMAIN(VAR(c, i)));

      // XXX: could count values and just return if all domains have
      //      nvariables values (or ...)
      while (_fd_val_has_next(iterator))
	{
	  int v = _fd_val_next_element(iterator) - min;

	  vgraph[i * r_verts + v] = 1;

	  if (r_edge[v] == UNSEEN)
	    {
	      r_edge[v]++;
	      nvalues++;
	    }
	}

      _fd_val_iterator_dispose(iterator);
    }

  // now that the exact number of different values in the domains is
  // known, check that there are enough
  if (nvalues < nvariables)
    return 0;

  queue = alloca((nvariables + r_verts) * sizeof(*queue));
  type = alloca((nvariables + r_verts) * sizeof(*type));

  // try to find a maximum cardinality matching
  for (i = 0; i < nvariables; ++i)
    {
      // find an M-augmenting path from variable i to an unmatched value
      memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
      memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

      l_parent[i] = -1;
      l_colour[i] = GREY;

      queue[head = tail = 0] = i;
      type[head] = VARIABLE;

      done = false;

      while (!done && head <= tail)
	{
	  switch (type[head])
	    {
	      case VARIABLE:		// node is a variable
		{
		  int v;

		  for (v = 0; v < r_verts; ++v)
		    {
		      if (!vgraph[queue[head] * r_verts + v])
			continue;

		      // this is a little faster than enqueueing all
		      // the values in the variable domain right away
		      if (r_edge[v] == -1)
			{
			  // found an M-augmenting path
			  r_parent[v] = queue[head];

			  while (v != -1)
			    {
			      r_edge[v] = r_parent[v];
			      l_edge[r_parent[v]] = v;

			      v = l_parent[r_parent[v]];
			    }

			  done = true;

			  break;
			}

		      if (r_colour[v] == WHITE)
			{
			  r_colour[v] = GREY;
			  queue[++tail] = v;
			  type[tail] = VALUE;
			  r_parent[v] = queue[head];
			}
		    }

		  //l_colour[queue[head]] = BLACK;

		  head++;
		}

		break;
	      case VALUE:		// node is a value (domain element)
		// the check for whether it is not yet matched has
		// already been made above; the next edge belongs to
		// the current matching
		if (l_colour[r_edge[queue[head]]] == WHITE)
		  {
		    l_colour[r_edge[queue[head]]] = GREY;
		    l_parent[r_edge[queue[head]]] = queue[head];
		    queue[++tail] = r_edge[queue[head]];
		    type[tail] = VARIABLE;
		  }

		//r_colour[queue[head]] = BLACK;

		head++;

		break;
	    }
	}

      if (!done)
	return 0;
    }

  // found a maximum cardinality matching; mark its edges in the value
  // graph
  for (i = 0; i < nvariables; ++i)
    vgraph[i * r_verts + l_edge[i]] = 2;

  _fd_mark_cycles(c, vgraph, r_verts, l_edge, r_edge, l_scc, r_scc, min);

  if (nvalues != nvariables)
    _fd_mark_paths(c, vgraph, r_verts, l_edge, r_edge, min);

  // remove unusable values from the variables domains
  // XXX: there's a slight performance gain when this is moved to
  //	  the end of _fd_mark_cycles()
  for (i = 0; i < nvariables; ++i)
    {
      int changed = 0;
      int v;

      for (v = 0; v < r_verts; ++v)
	if (vgraph[i * r_verts + v] == 1)
	  {
#ifdef DEBUG_MATCH
	    _fd_debug("removing %d from %d's domain\n", v + min,
		      VAR(c, i)->index);
#endif
	    changed |= _fd_var_del_val(v + min, VAR(c, i));

	    vgraph[i * r_verts + v] = 0;
	  }
	else if (vgraph[i * r_verts + v] > 2)
	  vgraph[i * r_verts + v] = 1;		// XXX: make sure it doesn't
						//	come back to bite us
						//	later

      if (changed)
	_fd_revise_connected(c, VAR(c, i));
    }

#ifdef CONSTRAINT_TEMPS
  *pmin = min;
  *(pmin + 1) = r_verts;
  *(pmin + 2) = nvalues;
#endif

  return 1;
}

int _fd_update_matching(fd_constraint c, fd_int culprit, int *memory)
{
  int nvariables = c->nvariables;
  int8_t *vgraph;
  int *l_edge, *r_edge;
  int min, r_verts, nvalues;
  int *l_scc, *r_scc;
  int ci;			// culprit index in the constraint
  bool affects_scc;
  int i;

  // retrieve contents of the constraint memory
  min = *(memory + 1);		// *memory contains the size of the block
  r_verts = *(memory + 2);
  nvalues = *(memory + 3);
  l_edge = memory + 4;
  r_edge = l_edge + nvariables;
  l_scc = r_edge + r_verts;
  r_scc = l_scc + nvariables;
  vgraph = (int8_t *) (r_scc + r_verts);

  // find out the index of the changed variable
  for (ci = 0; VAR(c, ci) != culprit; ++ci)
    ;

  // update the value graph
  for (i = 0; i < r_verts; ++i)
    vgraph[ci * r_verts + i] = -vgraph[ci * r_verts + i];

  {
    fd_iterator iterator = _fd_val_iterator(DOMAIN(culprit));

    while (_fd_val_has_next(iterator))
      {
	int v = _fd_val_next_element(iterator) - min;

	vgraph[ci * r_verts + v] = 1;
      }

    _fd_val_iterator_dispose(iterator);
  }

  for (i = 0; i < r_verts && vgraph[ci * r_verts + i] >= 0; ++i)
    ;

  // check that culprit's domain has really changed wrt the last
  // update of the value graph
  if (i == r_verts)
    {
      vgraph[ci * r_verts + l_edge[ci]] = 2;

      return 1;		// it hasn't
    }

  vgraph[ci * r_verts + i] = 0;
  affects_scc = r_scc[i] != -1;
  for (++i; i < r_verts; ++i)
    if (vgraph[ci * r_verts + i] < 0)
      {
	vgraph[ci * r_verts + i] = 0;

	if (r_scc[i] != -1)
	  affects_scc = true;
      }

  // if the previously matched value is still in the variable domain,
  // there's no need to rebuild the matching
  if (vgraph[ci * r_verts + l_edge[ci]])
    {
      vgraph[ci * r_verts + l_edge[ci]] = 2;

      if (affects_scc)
  	_fd_remark_cycles(c, ci, vgraph, r_verts, l_edge, r_edge, l_scc, r_scc, min);
    }
  else
    {
      // must try to update the maximum cardinality matching

      int *queue, head, tail;	// for the breadth-first search
      int8_t *l_colour, *r_colour;
      int8_t *type;		// type of node in queue: VARIABLE or VALUE
      int *l_parent, *r_parent;	// M-augmenting path
      bool done;

      for (i = 0; i < nvariables; ++i)
	if (i != ci)
	  assert(vgraph[i * r_verts + l_edge[i]] == 2),
	  vgraph[i * r_verts + l_edge[i]] = 1;

      r_edge[l_edge[ci]] = -1;	// the previously matched value is now free

      l_parent = alloca(nvariables * sizeof(*l_parent));
      l_colour = alloca(nvariables * sizeof(*l_colour));

      r_parent = alloca(r_verts * sizeof(*r_parent));
      r_colour = alloca(r_verts * sizeof(*r_colour));

      queue = alloca((nvariables + r_verts) * sizeof(*queue));
      type = alloca((nvariables + r_verts) * sizeof(*type));

      // find an M-augmenting path from culprit to an unmatched value
      memset(l_colour, WHITE, nvariables * sizeof(*l_colour));
      memset(r_colour, WHITE, r_verts * sizeof(*r_colour));

      l_parent[ci] = -1;
      l_colour[ci] = GREY;

      queue[head = tail = 0] = ci;
      type[head] = VARIABLE;

      done = false;

      while (!done && head <= tail)
	{
	  switch (type[head])
	    {
	    case VARIABLE:		// node is a variable
	      {
		int v;

		for (v = 0; v < r_verts; ++v)
		  {
		    if (!vgraph[queue[head] * r_verts + v])
		      continue;

		    // this is a little faster than enqueueing all
		    // the values in the variable domain right away
		    if (r_edge[v] == -1)
		      {
			// found an M-augmenting path
			r_parent[v] = queue[head];

			while (v != -1)
			  {
			    r_edge[v] = r_parent[v];
			    l_edge[r_parent[v]] = v;

			    v = l_parent[r_parent[v]];
			  }

			done = true;

			break;
		      }

		    if (r_colour[v] == WHITE)
		      {
			r_colour[v] = GREY;
			queue[++tail] = v;
			type[tail] = VALUE;
			r_parent[v] = queue[head];
		      }
		  }

		//l_colour[queue[head]] = BLACK;

		head++;
	      }

	      break;
	    case VALUE:			// node is a value (domain element)
	      // the check for whether it is not yet matched has
	      // already been made above; the next edge belongs to
	      // the current matching
	      if (l_colour[r_edge[queue[head]]] == WHITE)
		{
		  l_colour[r_edge[queue[head]]] = GREY;
		  l_parent[r_edge[queue[head]]] = queue[head];
		  queue[++tail] = r_edge[queue[head]];
		  type[tail] = VARIABLE;
		}

	      //r_colour[queue[head]] = BLACK;

	      head++;

	      break;
	    }
	}

      if (!done)
	return 0;

      // found a maximum cardinality matching; mark its edges in the value
      // graph
      for (i = 0; i < nvariables; ++i)
	vgraph[i * r_verts + l_edge[i]] = 2;

      _fd_remark_cycles(c, ci, vgraph, r_verts, l_edge, r_edge, l_scc, r_scc, min);
    }

  if (nvalues != nvariables)
    _fd_mark_paths(c, vgraph, r_verts, l_edge, r_edge, min);

  // remove unusable values from the variables domains
  // XXX: there's a slight performance gain when this is moved to
  //	  the end of _fd_mark_cycles()
  for (i = 0; i < nvariables; ++i)
    {
      int changed = 0;
      int v;

      for (v = 0; v < r_verts; ++v)
	if (vgraph[i * r_verts + v] == 1 && /* {XXX */ r_scc[v] == -1 /* XXX} */)
	  {
#ifdef DEBUG_MATCH
	    _fd_debug("removing %d from %d's domain\n", v + min,
		      VAR(c, i)->index);
#endif
	    changed |= _fd_var_del_val(v + min, VAR(c, i));

	    vgraph[i * r_verts + v] = 0;
	  }
	else if (vgraph[i * r_verts + v] > 2)
	  vgraph[i * r_verts + v] = 1;		// XXX: make sure it doesn't
						//	come back to bite us
						//	later

      if (changed)
	{
	  // XXX: needed because of the imperfect view of the state of
	  //	  the variables the function has
	  if (fd_domain_empty(VAR(c, i)))
	    return 0;

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

  return 1;
}