#include <stdio.h>
#include <string.h>
#include <fcntl.h>
#include <errno.h>

#include "fdc_int.h"

// (CSPLib 019)

#define NMAX 32

static int N = 5;		// actual order of the square

int main(int argc, char *argv[])
{
  fd_int square[NMAX][NMAX], all[NMAX * NMAX];
  fd_int numbers[NMAX];
  int sum;
  int gecode_symmetries = 0, old_symmetries = 0;
  int solutions = 0, one_solution = 1;
  int l, c, i, j;

  fd_init(&argc, &argv);

  for (i = 1; i < argc; ++i)
    if (!strcmp(argv[i], "--all"))
      one_solution = 0;
    else if (!strcmp(argv[i], "--gecode"))
      gecode_symmetries = 1;
    else if (!strcmp(argv[i], "--old"))
      old_symmetries = 1;
    else
      {
	N = atoi(argv[i]);

	break;
      }

  if (N > NMAX)
    _fd_fatal("magic-square: maximum order exceeded");

  j = 0;
  for (++i; i < argc; ++i)
    {
      int k0, k1;
      char *p;

      k0 = k1 = atoi(argv[i]) - 1;
      if (p = strchr(argv[i], '-'))
	k1 = atoi(p + 1) - 1;

      all[j++] = fd_new(k0, k1);
    }

  for (; j < N * N; ++j)
    all[j] = fd_new(0, N * N - 1);

  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      square[l][c] = all[l * N + c];

#if 01
  fd_all_different(all, N * N);
#else
  for (i = 0; i < N * N; ++i)
    fd_exactly(all, N * N, 1, i);	// fd_exactly_one(all, N * N, i);
#endif

  sum = (N * N - 1) * N * N / 2 / N;

#if 01
  // constrain lines
  for (l = 0; l < N; l++)
    fd_sum(square[l], N, sum);

  // constrain columns
  for (c = 0; c < N; c++)
    {
      for (l = 0; l < N; l++)
	numbers[l] = square[l][c];

      fd_sum(numbers, N, sum);
    }

  // constrain diagonals
  for (l = 0; l < N; l++)
    numbers[l] = square[l][l];

  fd_sum(numbers, N, sum);

  for (l = 0; l < N; l++)
    numbers[l] = square[l][N - 1 - l];

  fd_sum(numbers, N, sum);
#else
  fd_int sv = fd_new(sum, sum);

  // constrain lines
  for (l = 0; l < N; l++)
    fd_sum2(square[l], N, sv);

  // constrain columns
  for (c = 0; c < N; c++)
    {
      for (l = 0; l < N; l++)
	numbers[l] = square[l][c];

      fd_sum2(numbers, N, sv);
    }

  // constrain diagonals
  for (l = 0; l < N; l++)
    numbers[l] = square[l][l];

  fd_sum2(numbers, N, sv);

  for (l = 0; l < N; l++)
    numbers[l] = square[l][N - 1 - l];

  fd_sum2(numbers, N, sv);
#endif

  if (one_solution)
    {
      // break some symmetries

      if (gecode_symmetries)	  // from Gecode
	{
	  fd_gt(square[0][0], square[0][N - 1]);
	  fd_gt(square[0][0], square[N - 1][0]);
	}
      else if (old_symmetries)
	{
	  // disallow rotations and some reflections
	  fd_lt(square[0][0], square[0][N - 1]);
	  fd_lt(square[N - 1][0], square[0][N - 1]);
	  fd_lt(square[0][0], square[N - 1][N - 1]);
	}
      else
	{
	  // disallow reflections and rotations
#if 0
	  fd_lt(square[0][0], square[0][N - 1]);
	  fd_lt(square[0][N - 1], square[N - 1][0]);
#else
	  fd_lt(square[0][0], square[N - 1][0]);
	  fd_lt(square[N - 1][0], square[0][N - 1]);
#endif
	  fd_lt(square[0][0], square[N - 1][N - 1]);
	}
    }

  while (fd_solve() == FD_OK)
    {
      printf("solution %d:\n", ++solutions);

#if 0
      for (l = 0; l < N; ++l)
	for (c = 0; c < N; ++c)
	  if (!fd_var_single(square[l][c], NULL))
	    _fd_fatal("solution contains non-singleton variable");
#endif

      for (l = 0; l < N; ++l)
	{
	  for (c = 0; c < N; ++c)
	    printf(" %2d", fd_var_value(square[l][c]) + 1);

	  putchar('\n');
	}

#if !(defined(LOCAL_SEARCH) || defined(DISTRIBUTED_SOLVER))
      if (one_solution)
#endif
	break;
    }

  if (solutions)
    printf("%d solutions found\n", solutions);
  else
    printf("inconsistent CSP\n");

  fd_end();

  return !solutions;
}