#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stream.h>
#include <math.h>

#define Pi3(c, i)	((c/ThreeUp[i])%3)
#define Pi(c, i)	((c/NUp[i])%N)
#define PiNp(c, i)	((c/NpUp[i])%(N+1))
#define IotaPlus(c, i)	((N*(c/NUp[i])+N-1)*NUp[i]+c%NUp[i])

typedef char little;		// intended to hold values up to N*N
typedef short medium;		// intended to hold values up to N*2^N
typedef int large;		// intended to hold values up to N^2*N^N

struct clash_data {
  large b,c;
};

struct value {
  double val;
  int clashes;
};

little N=3;
medium TwoUpN=8;
medium *TwoUp;			// {1,2,4,8,...,2^N}
large ThreeUpN=27;
large *ThreeUp;			// {1,3,9,...,3^N}
large NmUpN=8;
large *NmUp;			// @N=3: {1,2,4,8}
large NUpN=27;
large *NUp;			// @N=3: {1,3,9,27}
large NpUpN=64;
large *NpUp;			// @N=3: {1,4,16,64}

little *SumOfDigitsModNp;	// Sum Of Digits of a base N+1 number of
				// N digits.

int **Binomial;			// The binomial coefficients

large *order;			// The order type of all N^N sums of
				// edge-lengths.

int verbose=1;			// verbose level
double T=1.;			// simulated annealing temperature
double CF=1.;			// Cooling Factor
int ladino_only=1;		// Restrict to ladino squares?
int compute_ladino_failure=0;	// Compute Ladino Failure?
int minimal_clash_dimension=0;  // Begin testing at clashes of this dimension.

large *default_order();
void read_cube(little **cube);
void print_cube(little **cube);
void copy_cube(little **src, little **dest);
void compute_coords(little **cube, large **corners, large **tops);
void print_corners(large **corners);
void print_tops(large **tops);
void make_checklist(
  large **corners, large **tops, int *cl_lens,
  large ***cl_low, large ***cl_high, clash_data **cl_cdata
);
value evaluate(
  int *cl_lens, large ***cl_low, large ***cl_high, clash_data **cl_cdata,
  int pc
);
void box_rotate(little **cube, large c, little a, little b);
void plane_swap(little **cube, little ind, little a, little b);
large ouroboros_move(little **cube, large start, little a, little b);

int main(int argc, char *argv[]) {

  int L=0;			// Number of attempts.
  int read_stdin=0;		// Read cube from stdin?
  int plane_swap_priority=10;
  int ouroboros_move_priority=10;
  int use_xor_configuration = 0;
  int use_Samorodnitsky = 0;
  double target_energy = 0.;
  int statistics_interval = 1;
  int print_clashes = 0;

  for (int carg=1; carg<argc; ++carg)
  if (0==strcmp(argv[carg], "-h")) {
    printf("Usage: ag [option]...
      -h       Print this message and exit successfully.
      -N num   Set the value of N to num.
      -2,-3,-4,-5,-6,-7  Set N to the corresponding value.
      -L num   Set the value of L to num.
      -T num   Set the value of T to num.
      -CF num  Set the Cooling Factor to num.
      -r       Read a cube from stdin.
      +ls      Set the ladino_only option to 1 (default).
      -ls      Set the ladino_only option to 0.
      +clf     Set compute_ladino_failure to 1.
      -clf     Set compute_ladino_failure to 0 (default)
      -ps num  Set the priority of a plane swap to num (default: 10).
      -om num  Set the priority of an ouroboros move to num (default: 10).
      -mcd num Begin testing at clashes of dimension num.
      -xor     Use the xor starting configuration (sets N=4).
      -as      Use Alex Samorodnitsky's solution (sets N=4).
      -te num  Set target_energy to num (default 0.).
      -si num  Set statistics_interval to num (default 1).
      -pc num  Set print_clashes to num (default 0).
      -s num   Set the seed of the random number generator(s).
      \n");
    return 0;
  } else if (0==strcmp(argv[carg], "-N")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -N missing\n");
      return 1;
    } else N=atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-2")) N=2;
  else if (0==strcmp(argv[carg], "-3")) N=3;
  else if (0==strcmp(argv[carg], "-4")) N=4;
  else if (0==strcmp(argv[carg], "-5")) N=5;
  else if (0==strcmp(argv[carg], "-6")) N=6;
  else if (0==strcmp(argv[carg], "-7")) N=7;
  else if (0==strcmp(argv[carg], "-L")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -L missing\n");
      return 1;
    } else L=atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-T")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -T missing\n");
      return 1;
    } else T=atof(argv[carg]);
  } else if (0==strcmp(argv[carg], "-CF")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -CF missing\n");
      return 1;
    } else CF=atof(argv[carg]);
  } else if (0==strcmp(argv[carg], "-r")) read_stdin=1;
  else if (0==strcmp(argv[carg], "+ls")) ladino_only = 1;
  else if (0==strcmp(argv[carg], "-ls")) ladino_only = 0;
  else if (0==strcmp(argv[carg], "+clf")) compute_ladino_failure = 1;
  else if (0==strcmp(argv[carg], "-clf")) compute_ladino_failure = 0;
  else if (0==strcmp(argv[carg], "-ps")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -ps missing\n");
      return 1;
    } else plane_swap_priority=atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-om")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -om missing\n");
      return 1;
    } else ouroboros_move_priority=atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-mcd")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -mcd missing\n");
      return 1;
    } else minimal_clash_dimension=atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-xor")) {
    use_xor_configuration = 1;
    N=4;
  } else if (0==strcmp(argv[carg], "-as")) {
    use_Samorodnitsky = 1;
    N=4;
  } else if (0==strcmp(argv[carg], "-te")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -te missing\n");
      return 1;
    } else target_energy = atof(argv[carg]);
  } else if (0==strcmp(argv[carg], "-si")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -si missing\n");
      return 1;
    } else statistics_interval = atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-pc")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -pc missing\n");
      return 1;
    } else print_clashes = atoi(argv[carg]);
  } else if (0==strcmp(argv[carg], "-s")) {
    if (argc<=++carg) {
      fprintf(stderr, "ag: argument to -s missing\n");
      return 1;
    } else {
      int s = atoi(argv[carg]);
      srand(s); srand48(s);
    }
  } else {
    fprintf(stderr, "ag: '%s' - unrecognized option.\n", argv[carg]);
    return 1;
  }
  if (verbose) {
    cout << "Running";
    for (int carg=0; carg<argc; ++carg) cout << " " << argv[carg];
    cout << "\n";
  }

  CF = exp(log(CF)/L);
  TwoUp = new medium[N+1]; TwoUpN = TwoUp[0] =1;
  ThreeUp = new large[N+1]; ThreeUpN = ThreeUp[0] =1;
  NmUp = new large[N+1]; NmUpN = NmUp[0] = 1;
  NUp = new large[N+1]; NUpN = NUp[0] = 1;
  NpUp = new large[N+1]; NpUpN = NpUp[0] = 1;
  for (little i=1; i<=N ; ++i) {
    TwoUp[i] = TwoUpN*=2;
    ThreeUp[i] = ThreeUpN*=3;
    NmUp[i] = NmUpN*=(N-1);
    NUp[i] = NUpN*=N;
    NpUp[i] = NpUpN*=(N+1);
  }
  if (verbose) printf("N=%d, TwoUpN=%d, NUpN=%d\n", N, TwoUpN, NUpN);
  SumOfDigitsModNp = new little[NpUpN];
  for (large c=0; c<=NpUpN; ++c) {
    little sum=0;
    for (little i=0; i<N; ++i) sum += PiNp(c,i);
    SumOfDigitsModNp[c] = sum;
  }
  Binomial = new (int *)[N+1];
  for (little n=0; n<=N; ++n) {
    Binomial[n] = new int[n+1];
    Binomial[n][0] = Binomial[n][n] = 1;
    for (little k=1; k<n ; ++k)
      Binomial[n][k] = Binomial[n-1][k-1] + Binomial[n-1][k];
  }
  order=default_order();

  little **cube;
  cube = new (little *)[NUpN];
  for (large c=0; c<NUpN; ++c) {
    cube[c] = new little[N];
    int s=0;
    for (little i=0; i<N; ++i) s += Pi(c,i);
    for (little i=0; i<N; ++i) cube[c][i] = (s+i)%N;
  }
  if (N==4) {
    plane_swap(cube, 3, 0, 1);
    plane_swap(cube, 3, 1, 3);
    if (use_xor_configuration) for (large c=0; c<NUpN; ++c) {
      int s=0;
      for (little i=0; i<N; ++i) s ^= Pi(c,i);
      for (little i=0; i<N; ++i) cube[c][i] = s^i;
    }
    if (use_Samorodnitsky) {
      int i0,t0,t1,c00,c10,c20,c30;
      for (large c=0; c<NUpN; ++c) for (little i=0; i<N; ++i) {
        t0 = 2&(Pi(c,0)^Pi(c,1));
        t1 = 2&(Pi(c,2)^Pi(c,3));
        if (t0) {
          c00 = 1&Pi(c,1);
          c10 = 1&Pi(c,0);
        } else {
          c00 = 1&Pi(c,0);
          c10 = 1&Pi(c,1);
        }
        if (t1) {
          c20 = 1&Pi(c,3);
          c30 = 1&Pi(c,2);
        } else {
          c20 = 1&Pi(c,2);
          c30 = 1&Pi(c,3);
        }
        i0 = 1&i;
        if (2&i) cube[c][i] = t1^(i0<<1); else cube[c][i] = t0^(i0<<1);
        if (cube[c][i]) cube[c][i] |= (i>>1)^c10^c30;
        else cube[c][i] |= (i>>1)^c00^c20;
      }
    }
  }
  if (N==5) {
    plane_swap(cube, 4, 0, 1);
    plane_swap(cube, 4, 1, 4);
  }
  if (read_stdin) read_cube(cube);
  print_cube(cube);

  large **corners;
  corners = new (large *)[NUpN];
  for (large c=0; c<NUpN; ++c) corners[c] = new large[N];
  large **tops;
  if (compute_ladino_failure) {
    tops = new (large *)[N];
    for (little i=0; i<N; ++i) tops[i] = new large[NUp[N-1]];
  }
  compute_coords(cube, corners, tops);
  // print_corners(corners);
  if (compute_ladino_failure) print_tops(tops);

  int *cl_lens;
  large ***cl_low, ***cl_high;
  clash_data **cl_cdata;
  cl_lens = new int[N];
  cl_low = new (large **)[N];
  cl_high = new (large **)[N];
  cl_cdata = new (clash_data *)[N];
  make_checklist(corners, tops, cl_lens, cl_low, cl_high, cl_cdata);
  value currval = evaluate(cl_lens, cl_low, cl_high, cl_cdata, 0);
  printf(
    "currval.val = %f, currval.clashes = %d\n", currval.val, currval.clashes
  );

  if (L>0) {
    cout << "Starting loop...\n";
    little **copy, **best;
    copy = new (little *)[NUpN];
    best = new (little *)[NUpN];
    for (large c=0; c<NUpN; ++c) {
      copy[c] = new little[N];
      best[c] = new little[N];
    }
    copy_cube(cube, best);
    value bestval = currval;
    int plane_swaps=0;
    int ouroboros_moves=0;
    int l=0;

    do {
      copy_cube(cube, copy);
      int do_plane_swap = 0;
      int do_ouroboros_move = 0;
      int do_box_rotate = 0;
      if (ladino_only) {
        if (
          rand() % (
            plane_swaps + plane_swap_priority + ouroboros_moves
            + ouroboros_move_priority
          ) < plane_swaps + plane_swap_priority
        )
          do_plane_swap = 1;
        else do_ouroboros_move = 1;
      } else do_box_rotate = 1;
      value newval;

      if (do_box_rotate) {
        large c = rand()%NUpN;
        little a = rand()%N;
        little b;
        do {b = rand()%N;} while (b==a);
        box_rotate(copy, c, a, b);
        compute_coords(copy, corners, tops);
        newval = evaluate(cl_lens, cl_low, cl_high, cl_cdata, 0);
        if (!(l%statistics_interval)) printf(
          "%d(%f,%f,%f): box rotate at %d, %d, %d: val=%f.\n",
          l, bestval.val, currval.val, T, c, a, b, newval.val
        );
      }

      if (do_plane_swap) {
        little ind = rand()%N;
        little a = rand()%N;
        little b;
        do {b = rand()%N;} while (b==a);
        plane_swap(copy, ind, a, b);
        compute_coords(copy, corners, tops);
        newval = evaluate(cl_lens, cl_low, cl_high, cl_cdata, 0);
        if (!(l%statistics_interval)) printf(
          "%d(%f,%f,%f): plane swap at %d, %d, %d: val=%f.\n",
          l, bestval.val, currval.val, T, ind, a, b, newval.val
        );
      } 

      if (do_ouroboros_move) {
        large start;
        little a,b;
        int done=0;
        int olen;
        do {
          start = rand()%NUpN;
          a = rand()%N;
          do {b = rand()%N;} while (b==a);
          olen = ouroboros_move(copy, start, a, b);
          if (rand()%olen < 5) done = 1;
          else {
            olen=ouroboros_move(copy, start, a, b);
          }
        } while (!done);
        compute_coords(copy, corners, tops);
        newval = evaluate(cl_lens, cl_low, cl_high, cl_cdata, 0);
        if (!(l%statistics_interval)) {
          printf(
            "%d(%f,%f,%f): ouroboros move at %d, %d, %d: length=%d, val=%f.\n",
            l, bestval.val, currval.val, T, start, a, b, olen, newval.val
          );
        }
      }

      if (
        newval.val<currval.val || drand48()<exp(-(newval.val-currval.val)/T)
      ) {
        plane_swaps += do_plane_swap;
        ouroboros_moves += do_ouroboros_move;
        if (newval.val != currval.val || !(l%statistics_interval)) {
          printf(
            "%d: moved to %f; plane_swaps=%d, ouroboros_moves=%d.\n",
            l, newval.val, plane_swaps, ouroboros_moves
          );
          if (print_clashes && newval.val == bestval.val)
            evaluate(cl_lens, cl_low, cl_high, cl_cdata, print_clashes);
        }
        currval = newval;
        copy_cube(copy, cube);
        if (currval.val < bestval.val) {
          copy_cube(cube, best);
          bestval = currval;
        }
      }
      T *= CF;
      ++l;
    } while (l < L && currval.val > target_energy);
    printf("Current cube is\n");
    print_cube(cube);
    printf("bestval.val=%f attained at\n", bestval.val);
    print_cube(best);
  }
}

large *default_order() {
  large *inv;
  inv = new large[NpUpN];
  large i,j;
  for (i=0; i<NpUpN; ++i) inv[i] = (SumOfDigitsModNp[i]*NpUpN)+i; 
  large t;
  for (j=NpUpN-1; j>0; --j) for (i=0; i<j; ++i) if (inv[i]>inv[i+1]) {
    t=inv[i]; inv[i]=inv[i+1]; inv[i+1]=t;
  }
  large *order;
  order = new large[NpUpN];
  for (i=0; i<NpUpN; ++i) order[inv[i] % NpUpN] = i;
  return order;
}

void read_cube(little **cube) {
  char x;
  for (large c=0; c<NUpN; ++c) for (little i=0; i<N; ++i) {
    do {cin >> x;} while (!(x>='0' && x<'0'+N));
    cube[c][i] = x - '0';
  }
}

void print_cube(little **cube) {
  for (large c=0; c<NUpN; ++c) {
    little i;
    if (c>0) {
      i=0;
      if (Pi(c,i) != 0) cout << " ";
      else if (Pi(c, ++i) != 0) cout << "  ";
      else do {
        cout << "\n";
      } while ((++i)<N && Pi(c,i)==0 );
    }
    for (i=0; i<N; ++i) cout << (int) cube[c][i];
  }
  cout << "\n";
}

void copy_cube(little **src, little **dest) {
  for (large c=0; c<NUpN; ++c) for (little i=0; i<N; ++i) dest[c][i]=src[c][i];
}

void compute_coords(little **cube, large **corners, large **tops) {
  large c,b;
  little i;
  for (c=0; c<NUpN; ++c) for (i=0; i<N; ++i) if (Pi(c,i)==0) corners[c][i]=0;
  else {
    b=c-NUp[i];
    corners[c][i] = corners[b][i]+NpUp[cube[b][i]];
  }
  if (compute_ladino_failure) for (i=0; i<N; ++i) for (c=0; c<NUp[N-1]; ++c) {
    b = IotaPlus(c, i);
    tops[i][c] = order[corners[b][i]+NpUp[cube[b][i]]];
  }
  for (c=0; c<NUpN; ++c) for (i=0; i<N; ++i)
    corners[c][i] = order[corners[c][i]];
}

void print_corners(large **corners) {
  for (large c=0; c<NUpN; ++c) {
    little i;
    if (c>0) {
      i=0;
      if (Pi(c,i) != 0) cout << " ";
      else if (Pi(c, ++i) != 0) cout << "  ";
      else do {
        cout << "\n";
      } while ((++i)<N && Pi(c,i)==0 );
    }
    for (i=0; i<N; ++i) {
      if (i>0) cout << "|";
      cout << (int) corners[c][i];
    }
  }
  cout << "\n";
}

void print_tops(large **tops) {
  large c;
  little i,j;
  for (i=0; i<N; ++i) {
    printf("top %d:\n", i);
    for (c=0; c<NUp[N-1]; ++c) {
      if (c>0) {
        j=0;
        if (Pi(c,j) != 0) cout << " ";
        else if (Pi(c, ++j) != 0) cout << "  ";
        else do {
          cout << "\n";
        } while ((++j)<N-1 && Pi(c,j)==0 );
      }
      cout << (int) tops[i][c];
    }
    cout << "\n";
  }
}

void make_checklist(
  large **corners, large **tops, int *cl_lens,
  large ***cl_low, large ***cl_high, clash_data **cl_cdata
) {
  little d;
  large ***low = new (large **)[N];
  large ***high = new (large **)[N];
  clash_data **cdata = new (clash_data *)[N];
  for (d=0; d<=N-2; ++d) {
    int len = cl_lens[d] = Binomial[N][d]*NUp[d]*NmUp[N-d]*TwoUp[N-d-1];
    low[d] = cl_low[d] = new (large *)[(N-d)*len];
    high[d] = cl_high[d] = new (large *)[(N-d)*len];
    cdata[d] = cl_cdata[d] = new clash_data[len];
  }
  for (large pert=0; pert<ThreeUpN; ++pert) {
    d=0;
    little fnz = 1;
    for (little i=0; i<N; ++i)
      if (Pi3(pert,i)==1) ++d; else if (fnz==1) fnz=Pi3(pert,i);
    if (d<=N-2 && fnz==0) for (large b=0; b<NUpN; ++b) {
      int ok=1;
      large c=b;
      for (little i=0; ok && i<N; ++i) if (Pi3(pert,i)==0)
        {if (ok = Pi(b,i)>0) c -= NUp[i];}
      else
        {if (Pi3(pert,i)==2) if (ok = Pi(b,i)<N-1) c += NUp[i];}
      if (ok) {
        cdata[d]->b = b;
        cdata[d]->c = c;
        cdata[d]++;
        for (little i=0; i<N; ++i) {
          if (Pi3(pert,i)==0) {
            *(low[d]++) = corners[b]+i;
            *(high[d]++) = corners[c+NUp[i]]+i;
          } else if (Pi3(pert,i)==2) {
            *(low[d]++) = corners[c]+i;
            *(high[d]++) = corners[b+NUp[i]]+i;
          }
        }
      }
    }
  }
  if (compute_ladino_failure) {
    cl_lens[N-1] = NUpN;
    cl_low[N-1] = new (large *)[NUpN];
    cl_high[N-1] = new (large *)[NUpN];
    for (little i=0; i<N; ++i) for (large c=0; c<NUp[N-1]; ++c) {
      cl_low[N-1][NUp[N-1]*i+c] = order+(NpUpN-1)/N;
      cl_high[N-1][NUp[N-1]*i+c] = tops[i]+c;
    }
  }

  if (verbose) {
    int total=0;
    for (d=0; d<=N-2; ++d) {
      total += cl_lens[d];
      printf(
        "cl_lens[%d]=%d, real: %d.\n",
        d, cl_lens[d], (high[d]-cl_high[d])/(N-d)
      );
    }
    if (compute_ladino_failure) {
      total += cl_lens[N-1];
      printf("cl_lens[%d]=%d.\n", N-1, cl_lens[N-1]);
    }
    printf("total = %d\n", total);
  }
}

value evaluate(
  int *cl_lens, large ***cl_low, large ***cl_high, clash_data **cl_cdata,
  int pc
) {
  value currval;
  currval.val=0.; currval.clashes=0;
  int maxd=N-2;
  if (compute_ladino_failure) maxd=N-1;
  for (little d=minimal_clash_dimension; d<=maxd; ++d) {
    large **low = cl_low[d];
    large **high = cl_high[d];
    for (int k=0; k<cl_lens[d]; ++k) {
      int clash=1;
      for (little i=0; i<N-d; ++i) if(**(low++) >= **(high++)) clash=0;
      if (clash) {
        currval.val += 1.;
        ++currval.clashes;
        if (currval.clashes<=pc) {
          printf("%d:",d);
          for (little i=0; i<N; ++i) printf("%d", Pi(cl_cdata[d][k].b, i));
          printf("-");
          for (little i=0; i<N; ++i) printf("%d", Pi(cl_cdata[d][k].c, i));
          printf(" ");
        }
      }
    }
  }
  if (pc) {
    if (currval.clashes > pc) printf("...");
    printf("\n");
  }
  return currval;
}

void box_rotate(little **cube, large c, little a, little b) {
  little t = cube[c][a];
  cube[c][a] = cube[c][b];
  cube[c][b] = t;
}

void plane_swap(little **cube, little ind, little a, little b) {
  large diff = NUp[ind]*(b-a);
  for (large c=0; c<NUpN; ++c) if (Pi(c, ind) == a) {
    little *t = cube[c];
    cube[c] = cube[c+diff];
    cube[c+diff] = t;
  }
}

large ouroboros_move(little **cube, large start, little a, little b) {
  large flips=0;
  large c=start;
  large c1;
  int found=1;
  little j=0;
  while (cube[c][j] != b) ++j;
  do {
    ++flips;
    little i=0;
    while (cube[c][i] != a) ++i;
    cube[c][j] = a;
    c1 = c - Pi(c,i)*NUp[i];
    j=0;
    while (j<N && cube[c1+j*NUp[i]][i] != b) ++j;
    cube[c][i] = b;
    if (j==N) found=0;
    c = c1+j*NUp[i];
    j=i;
  } while (found);
  return flips;
}