/*
 *  File:    LatinSquaresLT.c
 *  Author:  S Harry White
 *  Created: 2018-12-27
 *  Updated: 2022-01-10
 *    Replace & with && and | with || where necessary. Tidy code.
 *  Updated: 2023-01-18
 *    Change to compile with g++ for macOS.
 */

/*
 *  Adapted from a program written by Johan Öfverstedt that uses
 *  "Constraint-Based Local Search" and "Tabu Search" techniques:
 *
 *        "Water Retention on Magic Squares Solver",
 *        http://sourceforge.net/projects/wrmssolver/.
 *
 *  This program attempts to make Latin squares.
 *
 *  The type of Latin square can be specified as:
 *
 *    Latin, diagonal Latin, associative, or pandiagonal
 *
 *  Pandiagonal succeeds only for small orders.
 *
 */

#include <ctype.h>
#include <curses.h>
#include <errno.h>
#include <fcntl.h>
//#include <limits.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/stat.h>
//#include <sys/types.h>
#include <time.h>
#include <unistd.h>

int N,
    M,         // N-1
    O,         // N+1
    Nd2,       // N/2
    NN,        // N*N
    NNm1,      // NN-1
    NNm2,      // NN-2
    NNd2,      // NN/2
    SN;        // line sum=N*M/2

int *Q, *col, **colNums, *diagLnum, *diagRnum,
    *diagsL, *diagsR, **diagsLnums, **diagsRnums, *numbers, *tabu, tabuLen;
const bool F=false, T=true; bool pls, *diagL, *diagR, *all, allocDiags, allocPanDiags;

enum magicType { LS, DLS, ADLS, PLS, WPLS };
const enum magicType firstType=LS, lastType=WPLS; enum magicType squareType;
const char *squareTypeName[]={ "LS", "DLS", "ADLS", "PLS", "WPLS" };

void seed_rand() { srand((unsigned int)time(NULL)); }
//   ---------

int random_(int x) { return rand()%x; }
//  -------

const int bufSize=1024;
bool bell=T; char fmt[bufSize];
bool reportError(const char *msg) {
//   -----------
  printf("%sError: %s.\n", bell?"\a\n":"", msg); return bell=F;
} // reportError
//================================================== store ====================================================

void freeMatrix(int*** matrix, const int size) {
//   ----------
  if (*matrix!=NULL) { for(int i=0; i<size; i++) free((*matrix)[i]); free(*matrix); *matrix=NULL; }
} // freeMatrix

void freeStore() {
//   ---------
  if (Q!=NULL) { free(Q); Q=NULL; } if (col!=NULL) { free(col); col=NULL; } freeMatrix(&colNums, N);
  if (diagLnum!=NULL) { free(diagLnum); diagLnum=NULL; } if (diagRnum!=NULL) { free(diagRnum); diagRnum=NULL; }
  if (diagsL!=NULL) { free(diagsL); diagsL=NULL; } if (diagsR!=NULL) { free(diagsR); diagsR=NULL; }
  freeMatrix(&diagsLnums, N); freeMatrix(&diagsRnums, N); if (numbers!=NULL) { free(numbers); numbers=NULL; }
  if (tabu!=NULL) { free(tabu); tabu=NULL; } if (diagL!=NULL) { free(diagL); diagL=NULL; }
  if (diagR!=NULL) { free(diagR); diagR=NULL; } if (all!=NULL) { free(all); all=NULL; }
}

bool allocateMatrix(int*** matrix, const int size) {
//   --------------
  bool ok; *matrix=(int**) malloc(size*sizeof(int*)); ok=(*matrix!=NULL);
  if (ok) {
    int numAllocated=0;
    for (int i=0; i<size; i++) {
      int *p=(int*) malloc(size*sizeof(int)); (*matrix)[i]=p; if (p==NULL) { numAllocated=i; ok=F; break; }
    }
    if (!ok) freeMatrix(matrix, numAllocated);
  }
  return ok;
} // allocateMatrix

bool allocateInts(int **x, const int size) { return (*x=(int*)malloc(size*sizeof(int)))!=NULL; }

bool allocateBools(bool **x, const int size) { return (*x=(bool*)malloc(size*sizeof(bool)))!=NULL; }

bool allocateStore() {
//   -------------
  bool ok=T;
  if ((ok=allocateMatrix(&colNums, N)))
    if ((ok=allocateInts(&Q, NN)))
        if ((ok=allocateInts(&col, NN)))
          if ((ok=allocateInts(&tabu, NN)))
            if ((ok=allocateInts(&diagLnum, N)))
              if ((ok=allocateInts(&diagRnum, N)))
                if ((ok=allocateBools(&all, N))) ok=allocateInts(&numbers, N);
  if (((ok)&(allocDiags))) if ((ok=allocateBools(&diagL, NN))) ok=allocateBools(&diagR, NN);
  if (((ok)&(allocPanDiags)))
    if ((ok=allocateMatrix(&diagsLnums, N)))
      if ((ok=allocateMatrix(&diagsRnums, N)))
        if ((ok=allocateInts(&diagsL, NN))) ok=allocateInts(&diagsR, NN);
  if (!ok) { freeStore(); reportError("Storage allocation failed"); }
  return ok;
} // allocateStore
//=========================================== input =============================================

void clearLine(int c) { while (c!='\n') c=getchar(); }
//   ---------

bool getY() {
//   ----
  int c; do { c=getchar(); } while ((c==' ')|(c=='\t')|(c=='\n'));
  clearLine(c); return (c=='Y')|(c=='y');
} // getY

int getInt(const char *prompt) {
//  ------
  int n=0; printf("%s",prompt); scanf("%d", &n); clearLine(getchar());
  if (n<=0) reportError("Invalid input"); return n;
} // getInt
 
char *getLine(char *s, const char nl) {
//    -------
  int c, i=0; char *sin=s; do { c=getchar(); } while ((c==' ')|(c=='\t')|(c==nl));
  *s=c; if (c=='\n') return s; while (i++<bufSize) if ((*++s=getchar())=='\n') break;
  if (*s!='\n') { printf("\nType list too long.\n"); clearLine(*s); return NULL; }
  while ((*--s==' ')||(*s=='\t')); // strip any trailing whitespace
  *++s='\0'; s=sin; const int length=strlen(s);
  for (int i=0; i<length; ++i) s[i]=tolower(s[i]); return sin;
} // getLine

bool getType() {
//   -------
  printf("Type of square: Latin, diagonal Latin, associative, or pandiagonal, (as l, d, a, or p)? ");
  char buf[bufSize], *s=getLine(buf, '\n'); pls=(N%6==1)|(N%6==5); 
  if ((s!=NULL)&&(*s!='\0')) {
    while ((*s==' ')|(*s=='\t')) ++s;
    if (strlen(s)>=1) {
      switch (s[0]) {
        case 'a': squareType=ADLS; return T;
        case 'd': squareType=DLS; return T;
        case 'l': squareType=LS; return T;
        case 'p': squareType=pls?PLS:WPLS; return T;
        default: break;
      }
    }
  }
  printf("\a\nUnknown type: \"%s\".\n", s); return F;
} // getType
//============================================= output ==============================================

char OutputFileName[bufSize];
FILE *openOutput() {
//    ----------
  char buf[bufSize], defaultName[bufSize]; int sub=0; FILE *wfp=NULL;
  snprintf(buf, bufSize, "%s", squareTypeName[squareType]);
  snprintf(defaultName, bufSize, "./%dx%d%s", N, N, buf);
  snprintf(buf, bufSize, "%s.txt", defaultName);
  do {
    if (access(buf, F_OK)!=0) break;
    snprintf(buf, bufSize, "%s_%d.txt", defaultName, ++sub);
  } while (T);
  if ((wfp=fopen(buf, "w"))==NULL) {
    char msg[bufSize+50]; snprintf(msg, bufSize+50, "\a\nCan't open for write %s", buf); perror(msg);
  } else {
    printf("\nOutput file is %s\n", buf); strcpy(OutputFileName, buf);
  }
  return wfp;
} // openOutput

int fieldWidth() {
//  ----------
  if (N<=10) return 1; int i=N, width=1; while ((i=i/10)!=0) ++width; return width;
} // fieldWidth

typedef bool(*t_fprintFW)(FILE *fp, const int i);

bool fprintFW1(FILE *fp, const int i) { return fprintf(fp, "%1d", i)>0; }
bool fprintFW2(FILE *fp, const int i) { return fprintf(fp, "%2d", i)>0; }
bool fprintFW3(FILE *fp, const int i) { return fprintf(fp, "%3d", i)>0; }
bool fprintFW4(FILE *fp, const int i) { return fprintf(fp, "%4d", i)>0; }
bool fprintFW5(FILE *fp, const int i) { return fprintf(fp, "%5d", i)>0; }
bool fprintFW6(FILE *fp, const int i) { return fprintf(fp, "%6d", i)>0; }
bool fprintFW7(FILE *fp, const int i) { return fprintf(fp, "%7d", i)>0; }
bool fprintFW8(FILE *fp, const int i) { return fprintf(fp, "%8d", i)>0; }
bool fprintFW9(FILE *fp, const int i) { return fprintf(fp, "%9d", i)>0; }
bool fprintFWa(FILE *fp, const int i) { return fprintf(fp, "%10d", i)>0; }

static t_fprintFW fprintFW[]={ NULL,
  fprintFW1, fprintFW2, fprintFW3, fprintFW4, fprintFW5,
  fprintFW6, fprintFW7, fprintFW8, fprintFW9, fprintFWa
};
const int maxFieldWidth=10; int FW0, FW;

bool filePrint(int *x, FILE *wfp) {
//   ---------
  for (int i=0; i<N; i++) {
    if (!fprintFW[FW0](wfp, x[i*N])) return F;
    for (int j=1; j<N; j++) if (!fprintFW[FW](wfp, x[i*N+j])) return F;
    if (fputc('\n', wfp)==EOF) return F;
  }
  return fputc('\n', wfp)!=EOF;
} // filePrint

bool consolePrint(int *x) { return filePrint(x, stdout); }
//================================================= make =================================================

void initializeDiags() {
//   ---------------
  for (int i=0; i<N; ++i) {
    int k=i, l=i, c=N;
    while (c--) { diagsL[k]=i; k+=k==(NNm1-(i*N)) ? 1 : O; diagsR[l]=i; l+=l==(i*N) ? N+M : M; }
  }
} // initializeDiags

bool initialize() {
//   ----------
  M=N-1; O=N+1; Nd2=N/2; NN=N*N; NNm1=NN-1; NNm2=NN-2; NNd2=NN/2; SN=N*M/2;
  if ((squareType==LS)|(squareType==DLS)) tabuLen=2;
  else if (squareType==ADLS) tabuLen=N/2; else tabuLen=(N<9)?4:N/2;
  FW0=fieldWidth(); FW=FW0+1; if (FW>maxFieldWidth) return reportError("Output format overflow.");
  enum magicType t=squareType; allocDiags=(t!=LS); allocPanDiags=(t==PLS)|(t==WPLS);
  if (!allocateStore()) return F; for (int i=0; i<N; ++i) numbers[i]=i;
  for (int i=0; i<NN; ++i) {
    const int idN=i/N, imN=i%N; col[i]=imN;
    if (squareType!=LS) { diagL[i]=(idN==imN); diagR[i]=((idN+imN)==M); }
  }
  if (allocPanDiags) initializeDiags(); return T;
} // initialize

int lastViolation;
int updateViolationLS(const int i1, const int i2) {
//  -----------------
  const int c1=col[i1], c2=col[i2], x1=Q[i1], x2=Q[i2]; // swapped
  int c1x1=colNums[c1][x1], c1x2=colNums[c1][x2], c2x1=colNums[c2][x1], c2x2=colNums[c2][x2];
  const int oldv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x1; --c1x2; --c2x1; ++c2x2;
  const int newv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++colNums[c1][x1]; --colNums[c1][x2]; --colNums[c2][x1]; ++colNums[c2][x2];
  return lastViolation+=(newv-oldv);
} // updateViolationLS

int updateViolationDLS(const int i1, const int i2) {
//  ------------------
  const int c1=col[i1], c2=col[i2], x1=Q[i1], x2=Q[i2]; // swapped

  int c1x1=colNums[c1][x1], c1x2=colNums[c1][x2], c2x1=colNums[c2][x1], c2x2=colNums[c2][x2],
      oldv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x1; --c1x2; --c2x1; ++c2x2;
  int newv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++colNums[c1][x1]; --colNums[c1][x2]; --colNums[c2][x1]; ++colNums[c2][x2];
  if (diagL[i1]|diagL[i2]) {
    int dLx1=diagLnum[x1], dLx2=diagLnum[x2]; oldv+=((dLx1==0)?1:dLx1-1)+((dLx2==0)?1:dLx2-1);
    if (diagL[i1]) { ++dLx1; --dLx2; ++diagLnum[x1]; --diagLnum[x2]; }
    else { --dLx1; ++dLx2; --diagLnum[x1]; ++diagLnum[x2]; }
    newv+=((dLx1==0)?1:dLx1-1)+((dLx2==0)?1:dLx2-1); 
  }
  if (diagR[i1]|diagR[i2]) {
    int dRx1=diagRnum[x1], dRx2=diagRnum[x2]; oldv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
    if (diagR[i1]) { ++dRx1; --dRx2; ++diagRnum[x1]; --diagRnum[x2]; }
    else { --dRx1; ++dRx2; --diagRnum[x1]; ++diagRnum[x2]; }
    newv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1); 
  }
  return lastViolation+=(newv-oldv);
} // updateViolationDLS

int updateViolationSym(int i1, int i2) {
//  ------------------
  int c1=col[i1], c2=col[i2], x1=Q[i1], x2=Q[i2]; // swapped
  int c1x1=colNums[c1][x1], c1x2=colNums[c1][x2], c2x1=colNums[c2][x1], c2x2=colNums[c2][x2],
      oldv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x1; --c1x2; --c2x1; ++c2x2;
  int newv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++colNums[c1][x1]; --colNums[c1][x2]; --colNums[c2][x1]; ++colNums[c2][x2];
  if (diagR[i1]|diagR[i2]) {
    int dRx1=diagRnum[x1], dRx2=diagRnum[x2]; oldv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
    if (diagR[i1]) { ++dRx1; --dRx2; ++diagRnum[x1]; --diagRnum[x2]; }
    else { --dRx1; ++dRx2; --diagRnum[x1]; ++diagRnum[x2]; }
    newv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1); 
  }
  i1=NNm1-i1; i2=NNm1-i2; c1=col[i1]; c2=col[i2]; x1=Q[i1]; x2=Q[i2];
  c1x1=colNums[c1][x1]; c1x2=colNums[c1][x2]; c2x1=colNums[c2][x1]; c2x2=colNums[c2][x2],
  oldv+=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x1; --c1x2; --c2x1; ++c2x2;
  newv+=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++colNums[c1][x1]; --colNums[c1][x2]; --colNums[c2][x1]; ++colNums[c2][x2];
  if (diagR[i1]|diagR[i2]) {
    int dRx1=diagRnum[x1], dRx2=diagRnum[x2]; oldv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
    if (diagR[i1]) { ++dRx1; --dRx2; ++diagRnum[x1]; --diagRnum[x2]; }
    else { --dRx1; ++dRx2; --diagRnum[x1]; ++diagRnum[x2]; }
    newv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1); 
  }
  return lastViolation+=(newv-oldv);
} // updateViolationSym

int updateViolationPLS(const int i1, const int i2) {
//  ------------------
  const int c1=col[i1], c2=col[i2], l1=diagsL[i1], r1=diagsR[i1],
            l2=diagsL[i2], r2=diagsR[i2], x1=Q[i1], x2=Q[i2]; // swapped

  int a1x1=colNums[c1][x1], a1x2=colNums[c1][x2], a2x1=colNums[c2][x1], a2x2=colNums[c2][x2],
      oldv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x1; --a1x2; --a2x1; ++a2x2;
  int newv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++colNums[c1][x1]; --colNums[c1][x2]; --colNums[c2][x1]; ++colNums[c2][x2];
  a1x1=diagsLnums[l1][x1]; a1x2=diagsLnums[l1][x2]; a2x1=diagsLnums[l2][x1]; a2x2=diagsLnums[l2][x2];
  oldv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x1; --a1x2; --a2x1; ++a2x2;
  newv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++diagsLnums[l1][x1]; --diagsLnums[l1][x2]; --diagsLnums[l2][x1]; ++diagsLnums[l2][x2];
  a1x1=diagsRnums[r1][x1]; a1x2=diagsRnums[r1][x2]; a2x1=diagsRnums[r2][x1]; a2x2=diagsRnums[r2][x2];
  oldv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x1; --a1x2; --a2x1; ++a2x2;
  newv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++diagsRnums[r1][x1]; --diagsRnums[r1][x2]; --diagsRnums[r2][x1]; ++diagsRnums[r2][x2];
  return lastViolation+=(newv-oldv);
} // updateViolationPLS

int updateViolationWPLS(const int i1, const int i2) {
//  -------------------
  const int c1=col[i1], c2=col[i2], l1=diagsL[i1], r1=diagsR[i1],
            l2=diagsL[i2], r2=diagsR[i2], x1=Q[i1], x2=Q[i2], d=x1-x2; // swapped
  int a1x1=colNums[c1][x1], a1x2=colNums[c1][x2], a2x1=colNums[c2][x1], a2x2=colNums[c2][x2],
      oldv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x1; --a1x2; --a2x1; ++a2x2;
  int newv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++colNums[c1][x1]; --colNums[c1][x2]; --colNums[c2][x1]; ++colNums[c2][x2];
  oldv+=abs(diagsLnums[l1][0]-SN)+abs(diagsLnums[l2][0]-SN); diagsLnums[l1][0]+=d; diagsLnums[l2][0]-=d;
  newv+=abs(diagsLnums[l1][0]-SN)+abs(diagsLnums[l2][0]-SN);
  oldv+=abs(diagsRnums[r1][0]-SN)+abs(diagsRnums[r2][0]-SN); diagsRnums[r1][0]+=d; diagsRnums[r2][0]-=d;
  newv+=abs(diagsRnums[r1][0]-SN)+abs(diagsRnums[r2][0]-SN); return lastViolation+=(newv-oldv);
} // updateViolationWPLS

int violationLS() {
//  -----------
  int violation=0;
  for (int c=0; c<N; ++c) for (int v=0; v<N; ++v) {
    const int num=colNums[c][v]; violation+=(num==0)?1:num-1;
  }
  return lastViolation=violation;
} // violationLS

int violationDLS() {
//  ------------
  int violation=0;
  for (int c=0; c<N; ++c) {
    int num; for (int v=0; v<N; ++v) { num=colNums[c][v]; violation+=(num==0)?1:num-1; }
    num=diagLnum[c]; violation+=(num==0)?1:num-1; num=diagRnum[c]; violation+=(num==0)?1:num-1;
  }
  return lastViolation=violation;
} // violationDLS

int violationSym() {
//  ------------
  int violation=0;
  for (int c=0; c<N; ++c) {
    int num; for (int v=0; v<N; ++v) { num=colNums[c][v]; violation+=(num==0)?1:num-1; }
    num=diagRnum[c]; violation+=(num==0)?1:num-1;
  }
  return lastViolation=violation;
} // violationSym

int violationPLS() {
//  ------------
  int violation=0; // Column, pandiagonal constraints
  for (int i=0; i<N; ++i) {
    for (int j=0; j<N; ++j) {
      int num=colNums[i][j]; violation+=(num==0)?1:num-1; num=diagsLnums[i][j];
      violation+=(num==0)?1:num-1; num=diagsRnums[i][j]; violation+=(num==0)?1:num-1;
    }
  }
  return lastViolation=violation;
} // violationPLS

int violationWPLS() {
//  -------------
  int violation=0; // Column, pandiagonal constraints
  for (int i=0; i<N; ++i) {
    int num=0, Lsum=0, Rsum=0, *pL=Q+i, *pR=Q+i;
    for (int j=0; j<N; ++j) {
      num=colNums[i][j]; violation+=(num==0)?1:num-1;
      Lsum+=*pL; Rsum+=*pR; pL+=(i+j)==M?1:O; pR+=i==j?N+M:M;
    }
    diagsLnums[i][0]=Lsum; diagsRnums[i][0]=Rsum; violation+=abs(Lsum-SN)+abs(Rsum-SN);
  }
  return lastViolation=violation;
} // violationWPLS

void randomFillLS(int *X) {
//   ------------
  for (int i=0; i<N; ++i) X[i]=i; // nfr
  for (int r=1; r<N; ++r) { int si=r*N, ei=si+N;
    for (int i=0; i<N; ++i) {
      const int ri=i+random_(N-i), v=numbers[i]; numbers[i]=numbers[ri]; numbers[ri]=v;
    }
    for (int i=si, j=0; i<ei; ++i, ++j) X[i]=numbers[j];
  }
  for (int c=0; c<N; ++c) for (int i=0; i<N; ++i) colNums[c][i]=0;
  for (int i=0; i<NN; ++i) ++colNums[col[i]][X[i]];
} // randomFillLS

void randomFillDLS(int *X) {
//   -------------
  randomFillLS(X); for (int i=0; i<N; ++i) { diagLnum[i]=0; diagRnum[i]=0; }
  for (int i=0; i<NN; ++i) { if (diagL[i]) ++diagLnum[X[i]]; if (diagR[i]) ++diagRnum[X[i]]; }
} // randomFillDLS

void fillSym(int *X) {
//   -------
  for (int r=0; r<Nd2; ++r) { const int si=r*N, ei=si+N;
    for (int i=0; i<N; ++i) {
      const int ri=i+random_(N-i), v=numbers[i]; numbers[i]=numbers[ri]; numbers[ri]=v;
    }
    for (int i=si, j=0; i<ei; ++i, ++j) X[i]=numbers[j];
    for (int i=si; i<ei; ++i) if (X[i]==r) { X[i]=X[si+r]; X[si+r]=r; break; }
    for (int i=si; i<ei; ++i) X[NNm1-i]=M-X[i];
  }
  if (N&1) {
    for (int i=0; i<Nd2; ++i) numbers[i]=i;
    for (int i=0; i<Nd2; ++i) {
      const int ri=i+random_(Nd2-i), v=numbers[i]; numbers[i]=numbers[ri]; numbers[ri]=random_(2)?v:M-v;
    }
    for (int i=0; i<Nd2; ++i) numbers[M-i]=M-numbers[i]; numbers[Nd2]=Nd2;
    const int si=Nd2*N, ei=si+N; for (int i=si, j=0; i<ei; ++i, ++j) X[i]=numbers[j];
  }
} // fillSym

void randomFillSym(int *X) {
//   -------------
  fillSym(X);
  for (int c=0; c<N; ++c) for (int i=0; i<N; ++i) colNums[c][i]=0;
  for (int i=0; i<N; ++i) diagRnum[i]=0; 
  for (int i=0; i<NN; ++i) { ++colNums[col[i]][X[i]]; if (diagR[i]) ++diagRnum[X[i]]; }
} // randomFillSym

void randomFillPLS(int *X) {
//   -------------
  randomFillLS(X);
  for (int d=0; d<N; ++d) for (int i=0; i<N; ++i) { diagsLnums[d][i]=0; diagsRnums[d][i]=0; }
  for (int i=0; i<NN; ++i) { ++diagsLnums[diagsL[i]][X[i]]; ++diagsRnums[diagsR[i]][X[i]]; }
} // randomFillPLS

void randomFillWPLS(int *X) {
//   --------------
  randomFillLS(X); for (int d=0; d<N; ++d) { diagsLnums[d][0]=0; diagsRnums[d][0]=0; }
  for (int i=0; i<NN; ++i) { diagsLnums[diagsL[i]][0]+=X[i]; diagsRnums[diagsR[i]][0]+=X[i]; }
} // randomFillWPLS

int restartLS(int *X) { randomFillLS(X); return violationLS(); }
//  ---------

int restartDLS(int *X) { randomFillDLS(X); return violationDLS(); }
//  ----------

int restartSym(int *X) { randomFillSym(X); return violationSym(); }
//  -----------

int restartPLS(int *X) { randomFillPLS(X); return violationPLS(); }
//  ----------

int restartWPLS(int *X) { randomFillWPLS(X); return violationWPLS(); }
//  -----------

void doSwapLS(int *X, const int i1, const int i2) { const int t=X[i1]; X[i1]=X[i2]; X[i2]=t; }
//   --------

void doSwapSym(int *X, int i1, int i2) {
//   ---------
  int tmp=X[i1]; X[i1]=X[i2]; X[i2]=tmp; i1=NNm1-i1; i2=NNm1-i2; tmp=X[i1]; X[i1]=X[i2]; X[i2]=tmp;
} // doSwapSym

int swapDeltaLS(const int *X, const int i1, const int i2) {
//  -----------
  const int c1=col[i1], c2=col[i2], x1=X[i1], x2=X[i2];
  int c1x1=colNums[c1][x1], c1x2=colNums[c1][x2], c2x1=colNums[c2][x1], c2x2=colNums[c2][x2];
  const int oldv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x2; --c1x1; --c2x2; ++c2x1;
  const int newv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  return newv-oldv;
} // swapDeltaLS

int swapDeltaDLS(const int *X, const int i1, const int i2) {
//  ------------
  const int c1=col[i1], c2=col[i2], x1=X[i1], x2=X[i2];
  int c1x1=colNums[c1][x1], c1x2=colNums[c1][x2],c2x1=colNums[c2][x1], c2x2=colNums[c2][x2],
      oldv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x2; --c1x1; --c2x2; ++c2x1;
  int newv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  if (diagL[i1]|diagL[i2]) {
    int dLx1=diagLnum[x1], dLx2=diagLnum[x2]; oldv+=((dLx1==0)?1:dLx1-1)+((dLx2==0)?1:dLx2-1);
    if (diagL[i1]) { --dLx1; ++dLx2; } else { ++dLx1; --dLx2; }
    newv+=((dLx1==0)?1:dLx1-1)+((dLx2==0)?1:dLx2-1);
  }
  if (diagR[i1]|diagR[i2]) {
    int dRx1=diagRnum[x1], dRx2=diagRnum[x2]; oldv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
    if (diagR[i1]) { --dRx1; ++dRx2; } else { ++dRx1; --dRx2; }
    newv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
  }
  return newv-oldv;
} // swapDeltaDLS

int swapDeltaSym(const int *X, const int i1, const int i2) {
//  ------------
  const int c1=col[i1], c2=col[i2], x1=X[i1], x2=X[i2];
  int c1x1=colNums[c1][x1], c1x2=colNums[c1][x2], c2x1=colNums[c2][x1], c2x2=colNums[c2][x2], 
      oldv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  ++c1x2; --c1x1; --c2x2; ++c2x1;
  int newv=((c1x1==0)?1:c1x1-1)+((c1x2==0)?1:c1x2-1)+((c2x1==0)?1:c2x1-1)+((c2x2==0)?1:c2x2-1);
  if (diagR[i1]|diagR[i2]) {
    int dRx1=diagRnum[x1], dRx2=diagRnum[x2]; oldv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
    if (diagR[i1]) { --dRx1; ++dRx2; } else { ++dRx1; --dRx2; }
    newv+=((dRx1==0)?1:dRx1-1)+((dRx2==0)?1:dRx2-1);
  }
  return newv-oldv;
} // swapDeltaSym

int swapDeltaPLS(const int *X, const int i1, const int i2) {
//  ------------
  const int c1=col[i1], c2=col[i2], l1=diagsL[i1], r1=diagsR[i1], l2=diagsL[i2], r2=diagsR[i2],
            x1=X[i1], x2=X[i2];
  int a1x1=colNums[c1][x1], a1x2=colNums[c1][x2], a2x1=colNums[c2][x1], a2x2=colNums[c2][x2], 
      oldv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x2; --a1x1; --a2x2; ++a2x1;
  int newv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  a1x1=diagsLnums[l1][x1]; a1x2=diagsLnums[l1][x2]; a2x1=diagsLnums[l2][x1]; a2x2=diagsLnums[l2][x2]; 
  oldv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x2; --a1x1; --a2x2; ++a2x1;
  newv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  a1x1=diagsRnums[r1][x1]; a1x2=diagsRnums[r1][x2]; a2x1=diagsRnums[r2][x1]; a2x2=diagsRnums[r2][x2]; 
  oldv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x2; --a1x1; --a2x2; ++a2x1;
  newv+=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  return newv-oldv;
} // swapDeltaPLS

int swapDeltaWPLS(const int *X, const int i1, const int i2) {
//  -------------
  const int c1=col[i1], c2=col[i2], l1=diagsL[i1], r1=diagsR[i1], l2=diagsL[i2], r2=diagsR[i2],
            x1=X[i1], x2=X[i2];
  int a1x1=colNums[c1][x1], a1x2=colNums[c1][x2], a2x1=colNums[c2][x1], a2x2=colNums[c2][x2];
  const int oldv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  ++a1x2; --a1x1; --a2x2; ++a2x1;
  const int newv=((a1x1==0)?1:a1x1-1)+((a1x2==0)?1:a1x2-1)+((a2x1==0)?1:a2x1-1)+((a2x2==0)?1:a2x2-1);
  int x, delta=newv-oldv; const int y=x2-x1;
  x=diagsLnums[l1][0]-SN; delta+=(abs(x+y)-abs(x)); x=diagsLnums[l2][0]-SN; delta+=(abs(x-y)-abs(x));
  x=diagsRnums[r1][0]-SN; delta+=(abs(x+y)-abs(x)); x=diagsRnums[r2][0]-SN; delta+=(abs(x-y)-abs(x));
  return delta;
} // swapDeltaWPLS

int tabuSumsLS(int *X, const int iterations, int *violation) {
//  ----------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=NNm1; i>=N; --i) tabu[i]=0;
    while (it<iterations) {
      int i1=-1, i2=-1, bestDelta=0;
      for (int r=1; r<N; ++r) { int si=r*N, ej=si+M;
        for (int i=si; i<ej; ++i) if (tabu[i]<=it) {
          for (int j=i+1; j<=ej; ++j) {
            const int delta=swapDeltaLS(X, i, j);
            if ((i1<0)|(bestDelta>delta)) { i1=i; i2=j; bestDelta=delta; }
            else if ((delta==bestDelta)&&random_(2)) {  i1=i; i2=j; }
          } // for (int j ...
        } // for (int i ...
      } // for (int r ...
      ++it;
      if (i1>=0) {
        doSwapLS(X, i1, i2); tabu[i1]=it+tabuLen; if ((viol=updateViolationLS(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsLS

int tabuSumsDLS(int *X, const int iterations, int *violation) {
//  -----------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=NNm1; i>=N; --i) tabu[i]=0;
    while (it<iterations) {
      int i1=-1, i2=-1, bestDelta=0;
      for (int r=1; r<N; ++r) { int si=r*N, ej=si+M;
        for (int i=si; i<ej; ++i) if (tabu[i]<=it) {
          for (int j=i+1; j<=ej; ++j) {
            const int delta=swapDeltaDLS(X, i, j);
            if ((i1<0)|(bestDelta>delta)) { i1=i; i2=j; bestDelta=delta; }
            else if ((delta==bestDelta)&&random_(2)) { i1=i; i2=j; }
          } // for (int j ...
        } // for (int i ...
      } // for (int r ...
      ++it;
      if (i1>=0) {
        doSwapLS(X, i1, i2); tabu[i1]=it+tabuLen; if ((viol=updateViolationDLS(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsDLS

int tabuSumsSym(int *X, const int iterations, int *violation) {
//  -----------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=NNd2; i>=0; --i) tabu[i]=0;
    while (it<iterations) {
      int i1=-1, i2=-1, bestDelta=0; const int er=(N&1)?Nd2+1:Nd2;
      for (int r=0; r<er; ++r) { int si=r*N, ej=((N&1)&(r==(er-1)))?si+Nd2-1:si+M;
        for (int i=si; i<ej; ++i) if ((!diagL[i])&(tabu[i]<=it)) {
          for (int j=i+1; j<=ej; ++j) if (!diagL[j]) {
            const int delta=swapDeltaSym(X, i, j);
            if ((i1<0)|(bestDelta>delta)) { i1=i; i2=j; bestDelta=delta; }
            else if ((delta==bestDelta)&&random_(2)) { i1=i; i2=j; }
          } // for (int j ...
        } // for (int i ...
      } // for (int r ...
      ++it;
      if (i1>=0) {
        doSwapSym(X, i1, i2); tabu[i1]=it+tabuLen; if ((viol=updateViolationSym(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsSym

int tabuSumsPLS(int *X, const int iterations, int *violation) {
//  -----------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=NNm1; i>=N; --i) tabu[i]=0;
    while (it<iterations) {
      int i1=-1, i2=-1, bestDelta=0;
      for (int r=1; r<N; ++r) { int si=r*N, ej=si+M;
        for (int i=si; i<ej; ++i) if (tabu[i]<=it) {
          for (int j=i+1; j<=ej; ++j) {
            const int delta=swapDeltaPLS(X, i, j);
            if ((i1<0)|(bestDelta>delta)) { i1=i; i2=j; bestDelta=delta; }
            else if ((delta==bestDelta)&&random_(2)) { i1=i; i2=j; }
          } // for (int j ...
        } // for (int i ...
      } // for (int r ...
      ++it;
      if (i1>=0) {
        doSwapLS(X, i1, i2); tabu[i1]=it+tabuLen; if ((viol=updateViolationPLS(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsPLS

int tabuSumsWPLS(int *X, const int iterations, int *violation) {
//  ------------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=NNm1; i>=N; --i) tabu[i]=0;
    while (it<iterations) {
      int i1=-1, i2=-1, bestDelta=0;
      for (int r=1; r<N; ++r) { int si=r*N, ej=si+M;
        for (int i=si; i<ej; ++i) if (tabu[i]<=it) {
          for (int j=i+1; j<=ej; ++j) {
            const int delta=swapDeltaWPLS(X, i, j);
            if ((i1<0)|(bestDelta>delta)) { i1=i; i2=j; bestDelta=delta; }
            else if ((delta==bestDelta)&&random_(2)) { i1=i; i2=j; }
          } // for (int j ...
        } // for (int i ...
      } // for (int r ...
      ++it;
      if (i1>=0) {
        doSwapLS(X, i1, i2); tabu[i1]=it+tabuLen; if ((viol=updateViolationWPLS(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsWPLS
//============================================== main ================================================

typedef int(*t_iP)(int *); typedef int(*t_iPcIiP)(int *, const int, int *);
struct t_type { const char *name; t_iP restart; t_iPcIiP tabuSums; };
const int numTypes=lastType-firstType+1;
struct t_type squareTypes[numTypes]=
{
//  { squareTypeName[LS],   restartLS,   tabuSumsLS },
//  { squareTypeName[DLS],  restartDLS,  tabuSumsDLS },
//  { squareTypeName[ADLS], restartSym,  tabuSumsSym },
//  { squareTypeName[PLS],  restartPLS,  tabuSumsPLS },
//  { squareTypeName[WPLS], restartWPLS, tabuSumsWPLS }
  { "LS",   restartLS,   tabuSumsLS },
  { "DLS",  restartDLS,  tabuSumsDLS },
  { "ADLS", restartSym,  tabuSumsSym },
  { "PLS",  restartPLS,  tabuSumsPLS },
  { "WPLS", restartWPLS, tabuSumsWPLS }
};

const char *LSname="Latin", *DLSname="diagonal Latin", *assocName="associative",
           *plsName="pandiagonal", *wplsName="weakly pandiagonal", *ultraName="ultramagic",
           *errName="\"error\"";

bool isAssoc(int *X) { for (int i=0; i<NNd2; ++i) { if ((X[i]+X[NNm1-i]!=M)) return F; } return T; }
//   -------

bool isPLS(int *X) {
//   -----
  for (int i=0; i<N; ++i) { // Pandiagonals are transversals?
    int *pL=X+i, *pR=X+i; for (int k=0; k<N; ++k) all[k]=F;
    for (int j=0; j<N; ++j) { if (all[*pL]) return F; all[*pL]=T; pL+=(i+j)==M?1:O; }
    for (int k=0; k<N; ++k) all[k]=F;
    for (int j=0; j<N; ++j) { if (all[*pR]) return F; all[*pR]=T; pR+=i==j?N+M:M; }
  }
  return T;
} // isPLS

bool isWPLS(int *X) {
//   -----
  for (int i=0; i<N; ++i) { // Pandiagonal sums?
    int Lsum=0, Rsum=0, *pL=X+i, *pR=X+i;
    for (int j=0; j<N; ++j) { Lsum+=*pL; Rsum+=*pR; pL+=(i+j)==M?1:O; pR+=i==j?N+M:M; }
    if ((Lsum!=SN)|(Rsum!=SN)) return F;
  }
  return T;
} // isWPLS

bool isPan(int *X) { return pls ? isPLS(X) : isWPLS(X); }
//   -----

const char *checkDLS(int *X) {
//          --------
  if (isPan(X)) return isAssoc(X)?ultraName:pls?plsName:wplsName;
  else if (isAssoc(X)) return assocName; else return DLSname;
} // checkDLS

const char *checkLS(int *X) {
//          -------
  bool ok=T; int v; if ((!pls)&&isWPLS(X)) return wplsName;
  for (int i=0; i<N; ++i) all[i]=F;
  for (int i=0; i<N; ++i) { v=X[O*i]; if (all[v]) { ok=F; break; } all[v]=T; }
  if (ok) { for (int i=0; i<N; ++i) all[i]=F;
    for (int i=0; i<N; ++i) { v=X[M*(i+1)]; if (all[v]) { ok=F; break; } all[v]=T; }
  }
  return ok ? checkDLS(X) : LSname;
} // checkLS

const char *checkType(int *X) {
//          ---------
  if (N==1) return DLSname;
  switch (squareType) {
    case LS:   return checkLS(X);
    case DLS:  return checkDLS(X);
    case ADLS: return isPan(X)?ultraName:assocName;
    case PLS:  return isAssoc(X)?ultraName:plsName;
    case WPLS: return isAssoc(X)?ultraName:wplsName;
    default: break;
  }
  return errName;
} // checkType

bool checkN() {
//   ------
  if (N<=0) return reportError("Invalid order"); enum magicType t=squareType;
  if (((N==2)|(N==3))&(t!=LS)) {
    printf("\a\nThere are no %dx%d diagonal Latin squares.", N, N); return F;
  }
  if ((t==ADLS)&((N&3)==2))
    return reportError("There are no singly even associative Latin squares");
  if ((t==PLS)|(t==WPLS)) {
    if ((N&3)==2)
      return reportError("There are no singly even pandiagonal Latin squares");
    if (N==3) return reportError("There are no order 3 pandiagonal Latin squares");
  }
  return T;
} // checkN

//bool checkSquare(const int *X) {
////   -----------
//  for (int i=0; i<=NN; ++i) dupl[i]=F;
//  for (int i=0; i<NN; ++i) {
//    const int v=X[i];
//    if (dupl[v]) return reportError("Program error, duplicate value"); dupl[v]= T;
//  }
//  return T;
//} // checkSquare

void outputLocalTime() {
//   --------------
  time_t startTime=time(NULL); struct tm *local=localtime(&startTime);
  if (local!=NULL) {
    char dateTime[100];
    size_t slen=strftime(dateTime, 100, "%A %Y-%m-%d %X %Z\n\0", local); printf("\n%s", dateTime);
  }
} // outputLocalTime

int main() {
//  ----
  const int oneB=1000000000; bool ok=F; outputLocalTime(); seed_rand(); N=getInt("\nOrder? ");
  if (N>0) {
    if (getType()&&checkN()) { // checkN uses squareType
      if (initialize()) {
        struct t_type *sqType=&squareTypes[squareType];
        do {
          int runs=0, iterations=0, good=0, goodIt=0, BgoodIt=0, maxGoodIt=0;
          double elapsedTime=0., goodET=0., failET=0.;
          runs=getInt("Runs? "); if (runs<=0) break;
          iterations=getInt("Iterations? "); if (iterations<=0) break; FILE *wfp=openOutput();
          if (wfp!=NULL) {
            bool noSquare=T;
            for (int i=0; i<runs; ++i) {
              int violation=sqType->restart(Q), clock1=clock(),
                it=sqType->tabuSums(Q, iterations, &violation);
              double t=((clock()-clock1)/((double)CLOCKS_PER_SEC)); elapsedTime+=t;
              if (violation==0) {
                //if (!checkSquare(Q)) { consolePrint(Q); break; }
                ++good; goodIt+=it; if (goodIt>=oneB) { ++BgoodIt; goodIt-=oneB; }
                goodET+=t; if (it>maxGoodIt) maxGoodIt=it; noSquare=F;
                if ((i==0)|(i==(runs-1))) {
                  printf("\nRun: %d, iterations: %d, time: %.1fs.\n", i+1, it, t);
                  printf("Square is %s.\n\n", checkType(Q));
                  consolePrint(Q); if ((i==0)&(runs>1)) printf("   ...\n");
                }
                if (!filePrint(Q, wfp)) { reportError("Output failed"); break; } fflush(wfp);
              } else failET+=t;
            } // for (int i ...
            if (runs>0) {
              double avgT=elapsedTime/runs;
              if (good==0) { printf("\nRuns: %d, Squares: 0, Average Time: %.1fs", runs, avgT); }
              else {
                printf("\nRuns: %d, Squares: %d, Success Iterations: Avg %d Max %d\n", runs, good,
                  (int)(0.5+((double)(BgoodIt)*1E+9+(double)goodIt)/(double)good),maxGoodIt);
                printf("Average Time: %.1fs", avgT);
                if ((good<runs)&(avgT>=.1))
                  printf(", Avg Success Time: %.1fs, Avg Fail Time: %.1fs", goodET/good, failET/(runs-good));
              }
              printf("\n\n");
            }
            fclose(wfp); if (noSquare) remove(OutputFileName); else ok=T;
          } // if (wfp ...
          printf("Another? Y or N: "); if (!getY()) break;
        } while (T);
        freeStore();
      } // if (initialize ...
    } //if (getType ...
  }
  printf("\nHit return to close the console.");
  while (T) if (getchar()=='\n') break; return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main