/*
 *  File:    CompleteSquareSteps.cpp
 *  Author:  S Harry White
 *  Created: 2013-12-12
 *  Updated: 2021-12-20
 *    Tidy code.
 *  Updated: 2023-01-10
 *    Change to compile with g++ for macOS.
 */

/*
 *  Adapted from a program written by Johan Ofverstedt that uses
 *  "Constraint-Based Local Search" and "Tabu Search" techniques:
 *
 *        "Water Retention on Magic Squares Solver",
 *        http://sourceforge.net/projects/wrmssolver/.
 *
 *  Given an input rectangle with some cell values between 1 and RxC
 *  specified, (and other cell values 0), this program attempts to
 *  complete a magic rectangle by placing the remaining values in a
 *  specified maximum number of iterations. If successful, the
 *  starting rectangle and the rectangle after each iteration are output.
 *
 *  The type of rectangle can be specified as:
 *    semimagic, magic, associative
 *
 *  Note: If the rectangle is not square, magic equates to semimagic.
 */


#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>

#define min(a,b) (((a)<(b))?(a):(b))
#define max(a,b) (((a)>(b))?(a):(b))

int R, C,
    N,         // R==C ? R : sqrt(RC)
    Rm1, Cm1,  // R-1, C-1
    Rp1,       // R+1
    RC,        // R*C
    RCm1,      // RC-1
    RCd2,      // RC/2
    S2,        // RC+1
    S1,        // S2/2 for odd R, C
    SR,        // R*(RC+1)/2
    SC;        // C*(RC+1)/2

const int maxSteps=1000, maxSuccess=10, swapDataSize=9; const bool F=false, T=true;
int *xRectangle, *stepSqs, *stepSwap, *specific, *freeNum, *freeIx, *row, *col, *rowSum, *colSum,
    diagLsum, diagRsum, *tabu, tabuLen, Nfree, RCfree, RCfree_1, RCfree_2, numSteps, nextStep;

bool *diagL, *diagR, *numUsed, allZero, noZero, odd, singlyEven, fullSearch, retry; //, *dupl;

enum magicType { normalSemiMagic, normalMagic, Associative };
#define firstType normalSemiMagic
#define lastType Associative
enum magicType originalRectType, rectangleType;

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

bool reportError1(const char *msg, const int a) {
//   ------------
  snprintf(fmt, bufSize, "%sError: %s.\n", bell ? "\a\n" : "", msg); 
  printf(fmt, a); if (bell) bell=F; return F;
} // reportError1

bool reportError2(const char *msg, const int a, const int b) {
//   ------------
  snprintf(fmt, bufSize, "%sError: %s.\n", bell ? "\a\n" : "", msg); 
  printf(fmt, a, b); if (bell) bell=F; return F;
} // reportError2

bool reportError3(const char *msg, const int a, const int b, const int c) {
//   ------------
  snprintf(fmt, bufSize, "%sError: %s.\n", bell ? "\a\n" : "", msg); 
  printf(fmt, a, b, c); if (bell) bell=F; return F;
} // reportError3

bool reportError4(const char *msg, const int a, const int b, const int c, const int d) {
//   ------------
  snprintf(fmt, bufSize, "%sError: %s.\n", bell ? "\a\n" : "", msg); 
  printf(fmt, a, b, c, d); if (bell) bell=F; return F;
} // reportError4

void seed_rand() { srand((unsigned int)time(NULL)); }
//   ---------
int random_(const int x) { return rand()%x; }
//  -------
//============================================ store ==============================================

void freeStore() {
//   ---------
  if (xRectangle !=NULL) { free(xRectangle); xRectangle=NULL; }
  if (stepSqs    !=NULL) { free(stepSqs);    stepSqs   =NULL; }
  if (stepSwap   !=NULL) { free(stepSwap);   stepSwap  =NULL; }
  if (specific   !=NULL) { free(specific);   specific  =NULL; }
  if (freeNum    !=NULL) { free(freeNum);    freeNum   =NULL; }
  if (freeIx     !=NULL) { free(freeIx);     freeIx    =NULL; }
  if (row        !=NULL) { free(row);        row       =NULL; }
  if (col        !=NULL) { free(col);        col       =NULL; }
  if (diagL      !=NULL) { free(diagL);      diagL     =NULL; }
  if (diagR      !=NULL) { free(diagR);      diagR     =NULL; }
  if (rowSum     !=NULL) { free(rowSum);     rowSum    =NULL; }
  if (colSum     !=NULL) { free(colSum);     colSum    =NULL; }
  if (numUsed    !=NULL) { free(numUsed);    numUsed   =NULL; }
  if (tabu       !=NULL) { free(tabu);       tabu      =NULL; }
  //if (dupl       !=NULL) { free(dupl);       dupl      =NULL; }
} // freeStore

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; }
//   -------------

const char *allocFail="Storage allocation failed";
bool allocateStore() {
//   -------------
  bool ok=T;
  if ((ok=allocateInts(&xRectangle, RC)))
    if ((ok=allocateInts(&stepSqs, (maxSteps+1)*RC)))
      if ((ok=allocateInts(&stepSwap, (maxSteps+1)*swapDataSize)))
        if ((ok=allocateInts(&specific, RC)))
          if ((ok=allocateInts(&freeNum, RC)))
            if ((ok=allocateInts(&freeIx, RC)))
              if ((ok=allocateInts(&row, RC)))
                if ((ok=allocateInts(&col, RC)))
                  if ((ok=allocateInts(&rowSum, R)))
                    if ((ok=allocateInts(&colSum, C)))
                      if ((ok=allocateInts(&tabu, RC)))
                        ok=allocateBools(&numUsed, RC+1);
                          //ok=allocateBools(&dupl, RC+1);
  if (((ok)&(rectangleType==normalMagic))) {
    if ((ok=allocateBools(&diagL, RC))) ok=allocateBools(&diagR, RC);
  }
  if (!ok) { freeStore(); reportError(allocFail); } 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, const int min) {
//  -------
  printf("%s",prompt); int n=0; scanf("%d", &n); clearLine(getchar());
  if (n<min) reportError("Invalid input"); return n;
} // getInt

bool getInts(int *p, int *q, const int min) {
//   -------
  bool ok=F; int c; *p=*q=0; do { c=getchar(); } while ((c==' ')|(c=='\t'));
  if ( ('1'<=c)&(c<='9') ) {
    int i=c-'0'; while ( ('0'<=(c=getchar()))&&(c<='9') ) i=i*10+c-'0'; *p=i;
    if ((c==' ')|(c=='\t')) { do { c=getchar(); } while ((c==' ')|(c=='\t'));
      if ( ('1'<=c)&(c<='9') ) {
        int i=c-'0'; while ( ('0'<=(c=getchar()))&&(c<='9') ) i=i*10+c-'0'; *q=i;
      }
    }
    if (*q==0) *q=*p; ok=*p>=min;
  } 
  clearLine(c); if (!ok) return reportError("Invalid input"); return T;
} // getInts

bool getFileName(char *buf, const int size) {
//   -----------
  int c, i=0; char *s=buf;

  do { c=getchar(); } while ((c==' ')|(c=='\t')|(c=='\n') ); *s=c;
  while (i++<size) if ( (*++s=getchar())=='\n') break;
  if (*s!='\n') { printf("\nFile name too long.\n"); clearLine(*s); return F; }
  *s='\0'; return T;
} // getFileName

void adjustName(char *buf) {
//   ----------
  char *s=buf; bool txt=F;
  if ((*s!='.')&(*s!='/')) { // Assume in current folder.
    char tmp[bufSize]; snprintf(tmp, bufSize, "./%s", buf); strcpy(buf, tmp);
  }
  if (access(buf, R_OK)!=0) {
    while (*s++!='\0')
      ; --s;  if ((strlen(buf)>6)&&(*s--=='t')&&(*s--=='x')&&(*s--=='t')&&(*s=='.')) txt=T;
    if (!txt) strcat(buf, ".txt"); // no .txt entered, add it
  }
} // adjustName

FILE *openInput() {
//    ---------
  char buf[bufSize], *rFname=NULL; FILE *rfp=NULL;
  do {
    printf("\nEnter the name of the input file: ");
    if (getFileName(buf, bufSize-6)) { rFname=buf; break; }
    printf("\a\nCan't read the file name. Try again? y (yes) or n (no) "); if (!getY()) break;
  } while (T);
  if (rFname!=NULL) {
    adjustName(buf);
    if ((rfp=fopen(rFname, "r"))==NULL) {
      char msg[bufSize+50]; snprintf(msg, bufSize+50, "\a\nCan't open for read %s", buf); perror(msg);
    }
  }
  return rfp;
} // openInput

bool readInput(int *x) {
//   ---------
  FILE *rfp=openInput(); bool all0=T, no0=T; bool ok=(rfp!=NULL);
  if (ok) {
    for (int i=0; i<RC; ++i) {
      if ( fscanf(rfp, "%d", &x[i])==1) { if (x[i]==0) no0=F; else all0=F; }
      else { ok=reportError("Problem reading file"); break; }
    }
    fclose(rfp);
  }
  allZero=all0; noZero=no0; return ok;
} // readInput

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; // nothing
  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
//=============================================== output ============================================

bool openDir(char *dir, const int actualSteps, const char *rectName) {
//   -------
  char baseName[bufSize]; int fd=-1, sub=0;
  snprintf(baseName, bufSize, "./%dx%d%s%dSteps", R, C, rectName, actualSteps);
  strcpy(dir, baseName);
  do {
    if (mkdir(dir, 0775)==0) break;
    if (errno!=EEXIST) {
      char msg[bufSize+50]; snprintf(msg, bufSize+50, "Can't make folder %s", dir); perror(msg); return F;
    }
    snprintf(dir, bufSize, "%s_%d", baseName, ++sub);
  } while (T); 
  printf("Output files are in folder %s\n", dir); return T;
} // openDir

char runSummaryFile[bufSize];
FILE *openOutput(char *dir, const int step, const char *ext, const char *rectName) {
//    ----------
  FILE *wfp=NULL; bool saveName=F; int sub=0; char baseName[bufSize], buf[bufSize]; 
  if (step<0)
    if (dir==NULL) {
      saveName=T; snprintf(baseName, bufSize, "./%dx%d%sRuns", R,C,rectName);
    } else snprintf(baseName, bufSize, "%s/Swaps", dir);
  else snprintf(baseName, bufSize, "%s/%d", dir, step);
  snprintf(buf, bufSize, "%s.%s", baseName, ext);
  do {
    if (access(buf, F_OK)!=0) break;
    snprintf(buf, bufSize, "%s_%d.%s", baseName, ++sub, ext);
  } 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 { if (saveName) strcpy(runSummaryFile, buf); }
  return wfp;
} // openOutput

int fieldWidth() {
//  ----------
  if ((R==1)&(C==1)) return 1; int i=RC, 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<R; i++) { if (!fprintFW[FW0](wfp, x[i*C])) return F;
    for (int j=1; j<C; j++) if (!fprintFW[FW](wfp, x[i*C+j])) return F;
    if (fputc('\n', wfp)==EOF) return F;
  }
  return fputc('\n', wfp)!=EOF;
} // filePrint

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

void printSumsSemi(int *x) {
//   -------------
  int *pr=x;
	 // Row, column
  if (R==C) {
	   for (int i=0; i<R; ++i) {
	 	   int rSum=0, cSum=0; int *pc=x+i;
	 	   for (int j=0; j<C; ++j) { rSum+=*pr++; cSum+=*pc; pc+=C; }
	 	   printf("%5d:  row sum %8d   col sum %8d\n", i+1, rSum, cSum);
    }
    putchar('\n');
  } else {
    for (int i=0; i<R; ++i) {
	 	  int rSum=0; for (int j=0; j<C; ++j) rSum+=*pr++;
	 	  printf("%5d:  row sum %8d\n", i+1, rSum);
    }
    putchar('\n');
    for (int i=0; i<C; ++i) {
	 	  int cSum=0; const int *pc=x+i; for (int j=0; j<R; ++j) { cSum+=*pc; pc+=C; }
	 	  printf("%5d:  col sum %8d\n", i+1, cSum);
    }
    putchar('\n');
  }
} // printSumsSemi

void printSumsMagic(int *x) { // R==C
//   --------------
  int *pr=x; int Lsum=0, Rsum=0;
	// Row, column, diagonal
  for (int i=0; i<R; ++i) {
	  int rSum=0, cSum=0; int *pc=x+i, *pd=pr;
	  for (int j=0; j<R; ++j) { rSum+=*pr++; cSum+=*pc; pc+=R; }
	 	 printf("%5d:  row sum %8d   col sum %8d\n", i+1, rSum, cSum);
	 	Lsum+=pd[i]; Rsum+=pd[Rm1-i];
  }
  printf("      Ldiag sum %8d Rdiag sum %8d\n", Lsum, Rsum);
  putchar('\n');
} // printSumsMagic

enum t_ext { txt, html }; enum t_ext fileExt; const char *fileExtName[2]={ "txt", "html" };
bool getExt() {
//   ------
  printf("Output file type (txt, or html)? "); char buf[bufSize], *s=getLine(buf, '\n');
  if ((s!=NULL)&&(*s!='\0')) {
    while ((*s==' ')|(*s=='\t')) ++s;
    if (strlen(s)>=1) {
      if (s[0]=='h') { fileExt=html; return T; }
      if (s[0]=='t') { fileExt=txt;  return T; }
    }
  }
  printf("\a\nUnknown file type: \"%s\".\n", s); return F;
} // getExt

bool printStepSwaps(const int step, const int ia, const int a, const int ib, const int b,
//   --------------
                    const int ic, const int c, const int id, const int d, const int v, FILE *wfp) {
  if (step==0) { if (fprintf(wfp, "start ")<0) return F; }
  else {
    if (fprintf(wfp, "step %2d: swap s[%d][%d] %d, s[%d][%d] %d",
          step, row[ia]+1, col[ia]+1, a, row[ib]+1, col[ib]+1, b)<0) return F;
    if (ic<0) { if (fprintf(wfp, ": ")<0) return F; }
    else {
      if (fprintf(wfp, " and s[%d][%d] %d, s[%d][%d] %d: ",
            row[ic]+1, col[ic]+1, c, row[id]+1, col[id]+1, d)<0) return F;
    }
  }
  if (fprintf(wfp, "sums off %d\n", v)<0) return F;
  if ((step>0) &((step%5)==0)) if (fputc('\n', wfp)==EOF) return F; return T;
} // printStepSwaps

bool outSwaps(const int steps, FILE *wfp) {
//   --------
  bool ok=T;
  for (int i=0; i<=steps; ++i) {
    int *p=&stepSwap[i*swapDataSize];
    const int v=*p++, ia=*p++, a=*p++, ib=*p++, b=*p++, ic=*p++, c=*p++, id=*p++, d=*p;
    if (!printStepSwaps(i,ia,a,ib,b,ic,c,id,d,v,wfp)) { ok=F; break; }
  }
  if (ok) ok=fputc('\n', wfp)!=EOF; return ok;
} // outSwaps

bool writeHead(const int step, FILE *wfp) {
//   ---------
  return fprintf(wfp,
    "<!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML 1.0 Strict//EN\""
    " \"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd\">\n"
    "<html lang=\"en\" xml:lang=\"en\" "
    "xmlns=\"http://www.w3.org/1999/xhtml\">\n\n"
    "<head>\n"
    "<title>Step %d</title>\n"
    "<style type=\"text/css\" id=\"internalStyle\">\n"
    "  body { font-family : Verdana, Arial, Helvetica,"
    " sans-serif; color : #000080 }\n"
    "  td { font-size : small; font-weight : bold; text-align : center; }\n"
    "    .cellgg { background-color : #ccffcc; color : #008000 }\n"
    "    .cellrr { background-color : #ffcccc; color : #800000 }\n"
    "    .cellwb { background-color : #eeeeee; color : #000080 }\n"
    "</style>\n</head>\n\n<body>\n\n", step)>0;
} // writeHead

bool writeCell(FILE *wfp, const char *_class, const int value) {
//   ---------
  return fprintf(wfp, "<td class=\"%s\">%d</td>\n", _class, value)>0; }

int tableWidth() {
//  ----------
  const int cellWidth=C<10 ? 24 : C<32 ? 30 : C<100 ? 40 : 50; // ?? todo check
  return cellWidth*C;
} // tableWidth

bool writeRectangle(int *stepSq, int *swap, const int step, FILE *wfp) {
//   --------------
  const int width=tableWidth(); int *p=stepSq, *q=swap, a, b, c, d;
  q+=2; a=*q; q+=2; b=*q; q+=2; c=*q; q+=2; d=*q;

  if (fprintf(wfp,
        "<table width=\"%d\" summary=\"step %d\">"
        "<colgroup span=\"%d\" width=\"%d\"></colgroup>\n",
        width, step, C, width/C-2)<0) return F;  // ?? todo check

  for (int i=0; i<R; ++i) {  // ?? todo check
    if (fprintf(wfp, "<tr>\n")<0) return F;
   	for (int col=0; col<C; col++) {  // ?? todo check
      const int v=*p++;
      const char *s=((v==a)|(v==b)) ? "cellgg" : ((v==c)|(v==d)) ? "cellrr" : "cellwb";
      if (!writeCell(wfp, s, v)) return F;
   	}
    if (fprintf(wfp, "</tr>\n")<0) return F;
  }
  return fprintf(wfp, "</table>\n\n")>0;
} // writeRectangle

bool writeFoot(FILE *wfp) { return fprintf(wfp, "</body>\n</html>")>0; }
//   ---------

bool htmlPrint(int *stepSq, int *swap, const int step, FILE *wfp) {
//   ---------
  bool ok=writeHead(step, wfp);
  if (ok) ok=writeRectangle(stepSq, swap, step, wfp); if (ok) ok=writeFoot(wfp); return ok;
} // htmlPrint

bool outputSteps(const int actualSteps, const char *rectName) {
//   -----------
  char dir[bufSize]; bool ok=T; outSwaps(actualSteps, stdout);
  if (ok&&(ok=openDir(dir, actualSteps, rectName))) {
    FILE *wfp=openOutput(dir, -1, fileExtName[txt], rectName);
    if ((ok=outSwaps(actualSteps, wfp))) { fclose(wfp);
      for (int i=0; i<=actualSteps; ++i) {
        wfp=openOutput(dir, i, fileExtName[fileExt], rectName);
        if (fileExt==txt) ok=filePrint(stepSqs+(i*RC), wfp);
        else ok=htmlPrint(stepSqs+(i*RC), stepSwap+(i*swapDataSize), i, wfp); fclose(wfp);
      }
    }
  }
  return ok;
} // outputSteps

bool outputRun(const int run, const int it, const int viol, const double t,
//   ---------
               FILE *wfp, const char *rectName, const char *outType) {
  char buf[bufSize], tmp[bufSize]; if (run==1) putchar('\n');
  if (viol==0) snprintf(buf, bufSize, "Run %d, steps %d", run, it);
  else snprintf(buf, bufSize, "Run %d failed: sums off %d", run, viol);
  if (t>=0.05) { snprintf(tmp, bufSize, ", time: %.1fs", t); strcat(buf, tmp); }
  if (viol==0) {
    snprintf(tmp, bufSize, ": %s is %s", (R==C) ? "square" : "rectangle", outType);
    strcat(buf, tmp);
  }
  strcat(buf, "\n");
  printf("%s",buf); if (fprintf(wfp, "%s", buf)<0) return F;
  if (viol==0) {
    if (!(outputSteps(it, rectName)&&consolePrint(xRectangle)))
      return reportError("Output failed");
  }
  return T;
} // outputRun

bool outputSummary(const int runs, const int good, const int goodIt, const int maxGoodIt,
//   -------------
                   const double ET, const double goodET, const double failET, FILE *wfp) {
  char buf[bufSize]; const double avgT=ET/runs; const char *rect=(R==C) ? "Square" : "Rectangle";
  if (good==0) {
    snprintf(buf, bufSize, "\nRuns: %d, %ss: 0, Average Time: %.1fs", runs, rect, avgT);
  } else {
    char tmp[bufSize];
    snprintf(buf, bufSize, "\nRuns: %d, %ss: %d, Success Steps: Avg %d Max %d\n",
      runs, rect, good, goodIt/good, maxGoodIt);
    snprintf(tmp, bufSize, "Average Time: %.1fs", avgT); strcat(buf, tmp);
    if ((good<runs)&(avgT>=.05)) {
      snprintf(tmp, bufSize, ", Avg Success Time: %.1fs, Avg Fail Time: %.1fs",
        goodET/good, failET/(runs-good)); strcat(buf, tmp);
    }
  }
  strcat(buf, "\n"); printf("%s",buf); if (fprintf(wfp, "%s", buf)<0) return F; return T;
} // outputSummary
//=============================================== make =========================================

bool checkSpecific() {
//   -------------
  for (int i=0; i<=RC; ++i) numUsed[i]=F; bool ok=T;
  if (!allZero) {
    putchar('\n'); consolePrint(specific);
    for (int i=0; i<RC; ++i) {
      const int x=specific[i];
      if (x!=0) {
        if ((x<0)|(x>RC)) ok=reportError1("Number %d is out of range", x);
        else {
          if (numUsed[x]) ok=reportError1("%d specified more than once", x); numUsed[x]=T;
        }
      }
    }
  }
  return ok;
} // checkSpecific

bool getFree(void printSums(int *)) {
//   -------
  int next=0; if (noZero) allZero=T;
  if (allZero) {
    if (!noZero) printf("\nAll 0 input, numbers 1..%d free\n\n", RC);
    if (odd&(rectangleType==Associative)) {
      for (int i=0; i<RCd2; ++i) freeIx[next++]=i; for (int i=RCd2+1; i<RC; ++i) freeIx[next++]=i; }
    else { for (int i=0; i<RC; ++i) freeIx[next++]=i; }
    RCfree=next; next=0; for (int i=1; i<=RC; ++i) if (!numUsed[i]) freeNum[next++]=i;
  } else {
    printSums(specific); for (int i=0; i<RC; ++i) if (specific[i]==0) freeIx[next++]=i;
    RCfree=next; next=0; for (int i=1; i<=RC; ++i) if (!numUsed[i]) freeNum[next++]=i;
    if (RCfree!=0) {
      printf("\n\nfree:\n");
      for (int i=0; i<RCfree; ++i) { printf(" %4d", freeNum[i]); if ((i+1)%C==0) putchar('\n'); }
      putchar('\n'); putchar('\n');
    }
  }
  Nfree=(int)sqrt((double)RCfree+1.); RCfree_1=RCfree-1;
  RCfree_2=RCfree-2; tabuLen=((N==4)&allZero) ? 5 : (Nfree+Nfree)/3; return T;
} // getFree

bool initializeCommon() {
//   ----------------
  Rm1=R-1; Rp1=R+1; Cm1=C-1; RC=R*C; RCd2=RC/2; RCm1=RC-1;
  S2=RC+1; S1=S2/2; SR=C*S2/2; SC=R*S2/2; odd=R&1;
  singlyEven=((R&3)==2)&((C&3)==2); N=(R==C) ? R : (int)sqrt((double)RC);

  FW0=fieldWidth(); FW=FW0+1; if (FW>maxFieldWidth) return reportError("Output format overflow");

  if (!allocateStore()) return F;
  for (int i=0; i<RC; ++i) {
    const int idC=i/C, imC=i%C; row[i]=idC; col[i]=imC;
    if (rectangleType==normalMagic) { diagL[i]=(idC==imC); diagR[i]=((idC+imC)==Cm1); }
  }
  bool ok=readInput(specific)&&checkSpecific();
  originalRectType=rectangleType; return ok;
} // initializeCommon

bool checkSpecifiedRowsCols() {
//   ----------------------
  if (noZero) return T; bool ok=T; int *pr=specific;
	for (int i=0; i<R; ++i) {
	 	int rSum=0; bool rFull=T;
	  for (int j=0; j<C; ++j) { int rv=*pr++; if (rv==0) rFull=F; else rSum+=rv;}
	 	if ((rSum>SR)|(rFull&(rSum<SR))) {
      char gl=rSum>SR ? '>' : '<'; ok=reportError4("Row %d sum %d %c %d", i+1, rSum, gl, SR);
    }
  }
  for (int i=0; i<C; ++i) {
	  int cSum=0, *pc=specific+i; bool cFull=T;
	 	 for (int j=0; j<R; ++j) { const int cv=*pc; pc+=C; if (cv==0) cFull=F; else cSum+=cv; }
    if ((cSum>SC)|(cFull&(cSum<SC))) {
      char gl=cSum>SC ? '>' : '<'; ok=reportError4("Column %d sum %d %c %d", i+1, cSum, gl, SC);
    }
  }
  return ok;
} // checkSpecifiedRowsCols

bool checkSpecifiedRowsColsDiags() { // R==C
//   ---------------------------
  if (noZero) return T; bool ok=checkSpecifiedRowsCols();
  int *pd=specific, Lsum=0, RSum=0; bool Lfull=T, RFull=T;
	 for (int i=0; i<R; ++i) {
    const int d1v=pd[i], d2v=pd[Rm1-i];
    if (d1v==0) Lfull=F; else Lsum+=d1v; if (d2v==0) RFull=F; else RSum+=d2v; pd+=R;
  }
  if ((Lsum>SR)|(Lfull&(Lsum<SR))) {
    char gl=Lsum>SR ? '>' : '<'; ok=reportError3("Left diagonal sum %d %c %d", Lsum, gl, SR);
  }
  if ((RSum>SR)|(RFull&(RSum<SR))) {
    char gl=RSum>SR ? '>' : '<'; ok=reportError3("Right diagonal sum %d %c %d", RSum, gl, SR);
  }
  return ok;
} // checkSpecifiedRowsColsDiags

bool initializeSemi() {
//   --------------
  return initializeCommon()&&getFree(printSumsSemi)&&checkSpecifiedRowsCols();
} // initializeSemi

bool initializeMagic() {
//   ---------------
  return initializeCommon()&&getFree(printSumsMagic)&&checkSpecifiedRowsColsDiags();
} // initializeMagic

bool initializeAssoc() {
//   ---------------
  bool ok=initializeCommon();
  if (ok) {
    if (allZero) { if (odd) { specific[RCd2]=S1; numUsed[S1]=T; } }
    else {
      bool S1error=T, delta=F;
      for (int i=0; i<RC; ++i) {
        const int x=specific[i];
        if (x==0) { if (odd&(i==RCd2)) { specific[i]=S1; numUsed[S1]=T; delta=T; } }
        else {
          if (odd) {
            if (i==RCd2) {
              if (x!=S1) { ok=reportError2("Center value %d must be %d", x, S1); specific[i]=S1; }
            } else {
              if (S1error&(x==S1)) { ok=reportError1("%d must be center value", S1);
                S1error=F; // stop reporting again for complement cell
              }
            }
          }
          const int j=RCm1-i, y=specific[j];
          if (x!=(S2-y)) {
            if (y==0) { specific[j]=S2-x; numUsed[S2-x]=T; delta=T; }
            else { if (i<RCd2) ok=reportError2("Center symmetric %d and %d are not complementary", x, y); }
          }
        }
      }
      if (ok&(delta|noZero)) if (delta) ok=checkSpecific(); if (ok) ok=checkSpecifiedRowsCols();
    }
  }
  return ok ? getFree(R==C ? printSumsMagic : printSumsSemi) : F;
} // initializeAssoc

int lastViolation;
int updateViolationSemi(const int i1, const int i2) {
//  -------------------
  const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
	          x1=xRectangle[i1], x2=xRectangle[i2], d=x1-x2; // swapped

  int oldv=0, newv=0;
  if (r1!=r2) {
    oldv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR); rowSum[r1]+=d; rowSum[r2]-=d;
    newv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR);
  }
  if (c1!=c2) {
    oldv+=abs(colSum[c1]-SC)+abs(colSum[c2]-SC); colSum[c1]+=d; colSum[c2]-=d;
    newv+=abs(colSum[c1]-SC)+abs(colSum[c2]-SC);
  }
 	return lastViolation+=(newv-oldv);
} // updateViolationSemi

int updateViolationMagic(const int i1, const int i2) { // R==C
//  --------------------
  const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
	          x1=xRectangle[i1], x2=xRectangle[i2], d=x1-x2; // swapped

  int oldv=0, newv=0;
  if (r1!=r2) {
    oldv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR); rowSum[r1]+=d; rowSum[r2]-=d;
    newv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR);
  }
  if (c1!=c2) {
    oldv+=abs(colSum[c1]-SR)+abs(colSum[c2]-SR); colSum[c1]+=d; colSum[c2]-=d;
    newv+=abs(colSum[c1]-SR)+abs(colSum[c2]-SR);
  }
  if (r1==c1) {
    if (r2!=c2) { // First is in diag 1 but second is not.
      oldv+=abs(diagLsum-SR); diagLsum+=d; newv+=abs(diagLsum-SR);
    }
  } else if (r2==c2) { // First is not in diag 1 but second is.
    oldv+=abs(diagLsum-SR); diagLsum-=d; newv+=abs(diagLsum-SR);
  }

  if (r1==(Rm1-c1)) {
    if (r2!=(Rm1-c2)) { // First is in diag 2 but second is not.
      oldv+=abs(diagRsum-SR); diagRsum+=d; newv+=abs(diagRsum-SR);
    }
  } else if (r2==(Rm1-c2)) { // First is not in diag 2 but second is.
    oldv+=abs(diagRsum-SR); diagRsum-=d; newv+=abs(diagRsum-SR);
	 }
 	return lastViolation+=(newv-oldv);
} // updateViolationMagic

int updateViolationAssoc(int i1, int i2) {
//  --------------------
  int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2];
  const int x1=xRectangle[i1], x2=xRectangle[i2], d=x1-x2; // swapped

  int oldv=0, newv=0;
  if (r1!=r2) {
    oldv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR); rowSum[r1]+=d; rowSum[r2]-=d;
    newv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR);
  }
  if (c1!=c2) {
    oldv+=abs(colSum[c1]-SC)+abs(colSum[c2]-SC); colSum[c1]+=d; colSum[c2]-=d;
    newv+=abs(colSum[c1]-SC)+abs(colSum[c2]-SC);
  }

  i1=RCm1-i1;
  if (i1!=i2) { i2=RCm1-i2;
    r1=row[i1]; c1=col[i1]; r2=row[i2]; c2=col[i2];
    if (r1!=r2) {
      oldv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR); rowSum[r1]-=d; rowSum[r2]+=d;
      newv+=abs(rowSum[r1]-SR)+abs(rowSum[r2]-SR);
    }
    if (c1!=c2) {
      oldv+=abs(colSum[c1]-SC)+abs(colSum[c2]-SC); colSum[c1]-=d; colSum[c2]+=d;
      newv+=abs(colSum[c1]-SC)+abs(colSum[c2]-SC);
    }
  }
 	return lastViolation+=(newv-oldv);
} // updateViolationAssoc

int violationSemi() {
//  -------------
	int violation=0, *pr=xRectangle;
 // Row, column constraints
	for (int i=0; i<R; ++i) {
	  int rSum=0; for (int j=0; j<C; ++j) rSum+=*pr++; rowSum[i]=rSum; violation+=abs(rSum-SR);
  }
  for (int i=0; i<C; ++i) {
	 	 int cSum=0, *pc=xRectangle+i;
	 	 for (int j=0; j<R; ++j) { cSum+=*pc; pc+=C; } colSum[i]=cSum; violation+=abs(cSum-SC);
  }
 	return lastViolation=violation;
} // violationSemi

int violationMagic() { // R==C
//  --------------
  int violation=0, *pr=xRectangle, Lsum=0, Rsum=0;
  // Row, column, diagonal constraints
  for (int i=0; i<R; ++i) {
	  int rSum=0, cSum=0, *pc=xRectangle+i, *pd=pr;
	 	for (int j=0; j<R; ++j) { rSum+=*pr++; cSum+=*pc; pc+=R; }rowSum[i]=rSum; colSum[i]=cSum;
    violation+=abs(rSum-SR)+abs(cSum-SR);	Lsum+=pd[i]; Rsum+=pd[Rm1-i];
  }
  diagLsum=Lsum; diagRsum=Rsum; violation+=abs(Lsum-SR)+abs(Rsum-SR);
 	return lastViolation=violation;
} // violationMagic

void putStepRectangle(const int violation, const int ia, const int a, const int ib, const int b,
//   ----------------
                      const int ic, const int c, const int id, const int d) {
  int *p=xRectangle, *q=&stepSqs[nextStep*RC]; for (int i=0; i<RC; ++i) *q++=*p++;
  p=&stepSwap[nextStep*swapDataSize]; *p++=violation;
  *p++=ia; *p++=a; *p++=ib; *p++=b; *p++=ic; *p++=c; *p++=id; *p++=d; nextStep++;
} // putStepRectangle

int randomRestartMagic(int *X) {
//  ------------------
  if (noZero) { for (int i=0; i<RC; ++i) X[i]=specific[i]; }
  else {
    for (int i=0; i<(RCfree-1); ++i) {
      const int ri=i+random_(RCfree-i), tmp=freeNum[i]; freeNum[i]=freeNum[ri]; freeNum[ri]=tmp;
    }
    int next=0;
    for (int i=0; i<RC; ++i) if (specific[i]==0) X[i]=freeNum[next++]; else X[i]=specific[i];
  }
  const int v=rectangleType==normalMagic ? violationMagic() : violationSemi();
  putStepRectangle(v,-1,-1,-1,-1,-1,-1,-1,-1); return v;
} // randomRestartMagic

int randomRestartAssoc(int *X) {
//  ------------------
  if (noZero) { for (int i=0; i<RC; ++i) X[i]=specific[i]; }
  else {
    const int halfRCfree=RCfree/2;
    for (int i=0; i<(halfRCfree-1); ++i) {
      const int ri=i+random_(halfRCfree-i), tmp=freeNum[i]; freeNum[i]=freeNum[ri]; freeNum[ri]=tmp;
    }
    int next=0; if (odd) X[RCd2]=S1;
    for (int i=0; i<RCd2; ++i) {
      int v=specific[i]; if (v==0) { v=freeNum[next++]; if (random_(2)) v=S2-v; } X[i]=v; X[RCm1-i]=S2-v;
    }
  }
  const int v=violationSemi(); putStepRectangle(v,-1,-1,-1,-1,-1,-1,-1,-1); return v;
} // randomRestartAssoc

int doSwapMagic(int *X, const int ix1, const int ix2) {
//  -----------
  const int a=X[ix1], b=X[ix2], tmp=a; X[ix1]=b; X[ix2]=tmp;
  const int v=rectangleType==normalMagic ? updateViolationMagic(ix1,ix2)
                                         : updateViolationSemi(ix1,ix2);
  putStepRectangle(v,ix1,a,ix2,b,-1,-1,-1,-1); return v;
} // doSwapMagic

int doSwapAssoc(int *X, int ix1, int ix2) {
//  -----------
  const int ia=ix1, a=X[ia], ib=ix2, b=X[ib];
  int ic, c, id, d, tmp=a; X[ix1]=b; X[ix2]=tmp; ix1=RCm1-ix1;
  if (ix1!=ix2) { ix2=RCm1-ix2; c=X[ic=ix1]; d=X[id=ix2]; tmp=c; X[ix1]=d; X[ix2]=tmp; }
  else ic=c=id=d=-1;
  const int v=updateViolationAssoc(ia,ib);
  putStepRectangle(v,ia,a,ib,b,ic,c,id,d); return v;
} // doSwapAssoc

int swapDeltaSemi(const int *X, const int i1, const int i2) {
//  -------------
	const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
	          x1=xRectangle[i1], x2=xRectangle[i2], y=x2-x1; int x, delta=0;

  if (r1!=r2) { // Different rows
    x=rowSum[r1]-SR; delta+=abs(x+y)-abs(x); x=rowSum[r2]-SR; delta+=abs(x-y)-abs(x);
	}
  if (c1!=c2) { // Different columns
    x=colSum[c1]-SC; delta+=abs(x+y)-abs(x); x=colSum[c2]-SC; delta+=abs(x-y)-abs(x);
 	}
  return delta;
} // swapDeltaSemi

int swapDeltaMagic(const int *X, const int i1, const int i2) { // R==C
//  --------------
	const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
	          x1=xRectangle[i1], x2=xRectangle[i2], y=x2-x1; int x, delta=0;

	if (r1!=r2) { // Different rows
    x=rowSum[r1]-SR; delta+=abs(x+y)-abs(x); x=rowSum[r2]-SR; delta+=abs(x-y)-abs(x);
	}
  if (c1!=c2) { // Different columns
    x=colSum[c1]-SR; delta+=abs(x+y)-abs(x); x=colSum[c2]-SR; delta+=abs(x-y)-abs(x);
 	}
  if (r1==c1) {
    if (r2!=c2) { // First is in diag 1 but second is not.
      x=diagLsum-SR; delta+=abs(x+y)-abs(x);
    }
  } else if (r2==c2) { // First is not in diag 1 but second is.
    x=diagLsum-SR; delta+=abs(x-y)-abs(x);
  }
  if (r1==(Rm1-c1)) {
    if (r2!=(Rm1-c2)) { // First is in diag 2 but second is not.
      x=diagRsum-SR; delta+=abs(x+y)-abs(x);
    }
  } else if (r2==(Rm1-c2)) { // First is not in diag 2 but second is.
    x=diagRsum-SR; delta+=abs(x-y)-abs(x);
	 }
	 return delta;
} // swapDeltaMagic

int swapDeltaAssoc(int *X, const int i1, const int i2) {
//  --------------
  const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
	          x1=xRectangle[i1], x2=xRectangle[i2], y=x2-x1;
  int x, d, delta=0; const bool cPair=(RCm1-i1)==i2;

  if (r1!=r2) { // Different rows
    x=rowSum[r1]-SR; d=abs(x+y)-abs(x); delta+=d;
    if (cPair) delta+=d; else { x=rowSum[r2]-SR; delta+=abs(x-y)-abs(x); }
	}
	if (c1!=c2) { // Different columns
    x=colSum[c1]-SC; d=abs(x+y)-abs(x); delta+=d;
    if (cPair) delta+=d; else { x=colSum[c2]-SC; delta+=abs(x-y)-abs(x); }
 	}
 return delta;
} // swapDeltaAssoc

int tabuSumsSemi(int *X, const int iterations, int *violation) {
//  ------------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=RCm1; i>=0; --i) tabu[i]=0;
    while (it<iterations) {
		  int i1=-1, i2=-1, bestDelta=0, minR=-1, maxR=-1, minC=-1, maxC=-1,
          minRSum=INT_MAX, maxRSum=0, minCSum=INT_MAX, maxCSum=0;
      if (!fullSearch) {
        for (int i=0; i<R; ++i) {
          const int q=rowSum[i]; 
          if (q<minRSum) { minRSum=q; minR=i; } if (q>maxRSum) { maxRSum=q; maxR=i; }
        }
        for (int i=0; i<C; ++i) {
          const int q=colSum[i];
          if (q<minCSum) { minCSum=q; minC=i; } if (q>maxCSum) { maxCSum=q; maxC=i; }
        }
      }
      for (int i=RCfree_2; i>=0; --i) if (tabu[i]<=it) {
        const int fi=freeIx[i]; bool ido=fullSearch;
        if (!ido) {
          const int iR=row[fi], iC=col[fi]; ido=(iR==minR)||(iR==maxR)||(iC==minC)||(iC==maxC);
        }
        if (ido) for (int j=i+1; j<RCfree; ++j) if (tabu[j]<=it) {
          const int delta=swapDeltaSemi(X, fi, freeIx[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 ...
      ++it;
      if (i1>=0) {
		    if ((viol=doSwapMagic(X, freeIx[i1], freeIx[i2]))==0) break;
        tabu[i1]=it+tabuLen; if (Nfree>5) tabu[i2]=it+tabuLen;
      }
	   } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsSemi

int tabuSumsMagic(int *X, const int iterations, int *violation) {
//  -------------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=RCm1; i>=0; --i) tabu[i]=0;
    while (it<iterations) {
		  int i1=-1, i2=-1, bestDelta=0, minR=-1, maxR=-1, minC=-1, maxC=-1,
        minRSum=INT_MAX, maxRSum=0, minCSum=INT_MAX, maxCSum=0,
        minSum=INT_MAX, maxSum=0;  bool doL=F, doR=F;
      if (!fullSearch) {
        for (int i=0; i<N; ++i) {
          int q=rowSum[i]; 
          if (q<minRSum) { minRSum=q; minR=i; } if (q>maxRSum) { maxRSum=q; maxR=i; }
          q=colSum[i];
          if (q<minCSum) { minCSum=q; minC=i; } if (q>maxCSum) { maxCSum=q; maxC=i; }
        }
        minSum=min(minRSum,minCSum); maxSum=max(maxRSum,maxCSum);
        doL=(diagLsum<=minSum)|(diagLsum>=maxSum); doR=(diagRsum<=minSum)|(diagRsum>=maxSum);
      }
      for (int i=RCfree_2; i>=0; --i) if (tabu[i]<=it) {
        const int fi=freeIx[i]; bool ido=fullSearch;
        if (!ido) {
          const int iR=row[fi], iC=col[fi];
          ido=(iR==minR)||(iR==maxR)||(iC==minC)||(iC==maxC)||(doL&&diagL[fi])||(doR&&diagR[fi]);
        }
        if (ido) for (int j=i+1; j<RCfree; ++j) if (tabu[j]<=it) {
          const int delta=swapDeltaMagic(X, fi, freeIx[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 ...
      ++it;
      if (i1>=0) {
		    if ((viol=doSwapMagic(X, freeIx[i1], freeIx[i2]))==0) break;
        tabu[i1]=it+tabuLen; if (Nfree>5) tabu[i2]=it+tabuLen;
      }
	  } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsMagic

int tabuSumsAssoc(int *X, const int iterations, int *violation) {
//  -------------
  int it=0, viol=*violation;
  if (viol) {
    for (int i=RCm1; i>=0; --i) tabu[i]=0;
    while (it<iterations) {
		  int i1=-1, i2=-1, bestDelta=0, minR=-1, maxR=-1, minC=-1, maxC=-1,
          minRSum=INT_MAX, maxRSum=0, minCSum=INT_MAX, maxCSum=0;
      if (!fullSearch) {
        for (int i=0; i<R; ++i) {
          const int q=rowSum[i]; 
          if (q<minRSum) { minRSum=q; minR=i; } if (q>maxRSum) { maxRSum=q; maxR=i; }
        }
        for (int i=0; i<C; ++i) {
          const int q=colSum[i];
          if (q<minCSum) { minCSum=q; minC=i; } if (q>maxCSum) { maxCSum=q; maxC=i; }
        }
      }
      for (int i=RCfree_2; i>=0; --i) if (tabu[i]<=it) {
        const int fi=freeIx[i]; bool ido=fullSearch;
        if (!ido) {
          const int iR=row[fi], iC=col[fi]; ido=(iR==minR)||(iR==maxR)||(iC==minC)||(iC==maxC);
        }
        if (ido) for (int j=i+1; j<RCfree; ++j) if (tabu[j]<=it) {
          const int delta=swapDeltaAssoc(X, fi, freeIx[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 ...
      ++it;
      if (i1>=0) {
		    if ((viol=doSwapAssoc(X, freeIx[i1], freeIx[i2]))==0) break;
        tabu[i1]=it+tabuLen;
        if (Nfree>5) { tabu[i2]=it+tabuLen; tabu[RCfree_1-i1]=it+tabuLen; }
        tabu[RCfree_1-i2]=it+tabuLen;
      }
	   } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsAssoc
//================================================ main =============================================

typedef bool (*t_bno)();
typedef int  (*t_iip)(int *);
typedef int  (*t_iipciip)(int *, const int, int *);

struct t_type { const char *name; t_bno initialize; t_iip randomRestart; t_iipciip tabuSums; };
//const int numTypes=lastType-firstType+1;
#define numTypes lastType-firstType+1

struct t_type rectangleTypes[numTypes]=
{
  { "semimagic",   initializeSemi,  randomRestartMagic, tabuSumsSemi  },
  { "magic",       initializeMagic, randomRestartMagic, tabuSumsMagic },
  { "associative", initializeAssoc, randomRestartAssoc, tabuSumsAssoc }
};

bool getType() {
//   -------
  printf("Type of rectangle, (semimagic, magic, associative), as s, m, or a? ");
  char buf[bufSize], *s=getLine(buf, '\n');
  if ((s!=NULL)&&(*s!='\0')) {
    while ((*s==' ')|(*s=='\t')) ++s;
    if (strlen(s)>=1) {
      switch (s[0]) {
        case 'a': rectangleType=Associative; return T;
        case 'm': rectangleType=R==C ? normalMagic : normalSemiMagic; return T;
        case 's': rectangleType=normalSemiMagic; return T;
        default: break;
      }
    }
  }
  printf("\a\nUnknown type: \"%s\".\n", s); return F;
} // getType

const char *getName() {
//          -------
  enum magicType t=((R!=C)&&(rectangleType==normalSemiMagic)) ? normalMagic : rectangleType;
  return rectangleTypes[t].name;
} // getName

const char *magicName="normal magic", *semiName="normal semimagic",
     *assocName="associative", *ultraName="ultramagic", *errorName="\"error\"";
bool isAssoc(int *x) { for (int i=0; i<RCd2; ++i) { if ((x[i]+x[RCm1-i]!=S2)) return F; } return T; }
//   -------

bool isPan(int *x) {
//   -----
  int *pL, *pR;
	// Pandiagonal sums
	for (int i=0; i<R; ++i) {
	  int Lsum=0, Rsum=0; pL=pR=x+i;
	  for (int j=0; j<R; ++j) { Lsum+=*pL; Rsum+=*pR; pL+=(i+j)==Rm1 ? 1 : Rp1; pR+=i==j ? R+Rm1 : Rm1; }
    if ((Lsum!=SR)|(Rsum!=SR)) return F;
  }
  return T;
} // isPan

const char *checkMagic(int *x) { return isAssoc(x) ? isPan(x) ? ultraName : assocName : magicName; }
//          ----------

const char *checkSemiMagic(int *x) {
//          --------------
  if (R==C) {
    int Lsum=0, Rsum=0; for (int i=0; i<R; ++i) { Lsum+=x[Rp1*i]; Rsum+=x[Rm1*(i+1)]; }
    return ((Lsum==SR)&(Rsum==SR)) ? checkMagic(x) : semiName;
  }
  return isAssoc(x) ? assocName : magicName;
} // checkSemiMagic

const char *checkType(int *x) {
//          ---------
  if ((R==1)&(C==1)) return magicName;
  switch (rectangleType) {
    case normalSemiMagic: return checkSemiMagic(x);
    case normalMagic:     return checkMagic(x);
    case Associative:     return assocName;
    default: break;
  }
  return errorName;
} // checkType

bool checkRC() {
//   ---------
  if ((R<=0)|(C<=0)) return reportError("Invalid order"); enum magicType t=rectangleType;
  if (((R==1)|(C==1))&(R!=C)) {
    printf("\a\nThere are no %dx%d magic rectangles.\n", R, C); return F;
  }
  if ((R==2)|(C==2)) {
    if (R==C) {
      printf("\a\nThere are no 2x2 %smagic squares.", t==normalSemiMagic ? "semi" : ""); return F;
    }
    if (t==Associative) {
      printf("\a\nThere are no %dx%d associative magic rectangles.\n", R, C); return F;
    }
  }
  if ((R&1)!=(C&1)) return reportError("There are no odd, even magic rectangles");
  if ((t==Associative)&((R&3)==2)&((C&3)==2))
    return reportError("There are no singly even associative magic rectangles");
  return T;
} // checkRC

bool useAssoc() {
//   --------
  if (!(singlyEven)) {
    bool centerSym=T;
    for (int i=0; i<RCd2; ++i) if ((specific[i]+specific[RCm1-i])!=S2) centerSym=F;
    if (centerSym) { rectangleType=Associative; return T; }
  }
  return F;
} // useAssoc

void initRetry(struct t_type **sqType, FILE *wfp) {
//   ---------
  retry=F; fullSearch=T; *sqType=&rectangleTypes[Associative];
  char buf[bufSize]; strcpy(buf, "\nRetry with full search");
  if (originalRectType!=Associative) strcat(buf, " associative");
  strcat(buf, ".\n"); printf("%s",buf); fprintf(wfp, "%s", buf);
} //

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

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\n\0", local); printf("\n%s", dateTime);
  }
} // outputLocalTime

int main() {
//  ----
  bool ok=F; outputLocalTime(); seed_rand(); printf("\nOrder? ");
  if (getInts(&R, &C, 1))
  if (getType()&&checkRC()&&getExt()) { // checkRC uses rectangleType
    struct t_type *sqType=&rectangleTypes[rectangleType];
    printf("Stop at first success? Y or N: "); const bool stopAt1=getY();

    if (sqType->initialize()) {
      do { // while another
        const int runs=getInt("runs? ", 1); if (runs<=0) break;
        numSteps=getInt("Maximum steps? ", 1); if (numSteps<=0) break;
        if (numSteps>maxSteps) { reportError1("Step limit is set at %d", maxSteps); break; }
        FILE *wfpCons=NULL;
        if ((wfpCons=openOutput(NULL,-1,fileExtName[txt], getName()))!=NULL) {
          int stopRun=0; bool noRectangle=T; retry=noZero&&useAssoc();
          do { // once and retry?
            int good=0, goodIt=0, maxGoodIt=0; noRectangle=T; double elapsedT=0., goodET=0., failET=0.;
	          for (int i=0; i<runs; ++i) {
              nextStep=0; stopRun=i;
              int violation=sqType->randomRestart(xRectangle), clock1=clock();
              const int it=sqType->tabuSums(xRectangle, numSteps, &violation);
              const double t=((clock()-clock1)/((double)CLOCKS_PER_SEC)); elapsedT+=t;
              if (!outputRun(i+1, it, violation, t, wfpCons, getName(), checkType(xRectangle))) break;
              if (violation==0) { ok=T;
                //if (!checkRectangle(xRectangle)) break;
                ++good; goodIt+=it; goodET+=t; noRectangle=F;
                if (it>maxGoodIt) maxGoodIt=it; if (stopAt1) break;
                if (good==maxSuccess) { printf("Stopping at %d successes.\n", maxSuccess); break; }
              } else failET+=t;
            } // for (int i ...
            if (runs>0) {
              if (!outputSummary(++stopRun, good, goodIt, maxGoodIt,
                                 elapsedT, goodET, failET, wfpCons)) break;
            }
            if (!(noRectangle&retry)) break; initRetry(&sqType, wfpCons);
          } while (T); // once and retry?
          fclose(wfpCons);
          if (stopRun>1) printf("\nRun summary file is %s\n\n", runSummaryFile);
          else { putchar('\n'); remove(runSummaryFile); }
	       } // if ((wfpCons ...
        printf("Another? Y or N: "); if (!getY()) break;
        rectangleType=originalRectType; // possibly reset for retry
	    } while (T); // while another
      freeStore();
    } // if (sqType.initialize ...
  } //if (getType ...
  printf("\nHit return to close the console.");
  while (T) if (getchar()=='\n') break; return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main