/*
 *  File:    CompleteSquare.cpp
 *  Author:  S Harry White
 *  Created: 2012-08-28
 *  Updated: 2021-12-21
 *    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 RxC 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.
 *
 *  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>

int R, C,
    N,         // sqrt(R*C)
    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

    *xRectangle, *specific, *freeNum, *freeIx,
    *row, *col, *rowSum, *colSum, diagLsum, diagRsum,
    *tabu, tabuLen, Nfree, RCfree, RCfree_1, RCfree_2;
const bool F=false, T=true;
bool *diagL, *diagR, *numUsed, allZero, noZero, odd, fullSearch, singlyEven, retry; // *dupl,

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

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

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

#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
//============================================ store ==============================================

void freeStore() {
//   ---------
  if (xRectangle!=NULL) { free(xRectangle); xRectangle=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 (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 (diagL     !=NULL) { free(diagL);      diagL     =NULL; }
  if (diagR     !=NULL) { free(diagR);      diagR     =NULL; }
  //if (dupl      !=NULL) { free(dupl);       dupl       = NULL; }
}

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=allocateInts(&xRectangle, RC)))
    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("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) {
//  -------
  printf("%s",prompt); int n=0; scanf("%d", &n); clearLine(getchar());
  if (n<=0) reportError("Invalid input"); return n;
} // getInt

bool getInts(int *p, int *q) {
//   -------
 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=T;
  } 
  clearLine(c); if (!ok) return reportError("Invalid input"); return T;
} // getInts

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

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; }
    else { 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
//============================================ output ========================================

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) { // Row, column
//   -------------
  int *pr=x;
  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; 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 // Row, column, diagonal
//   --------------
   int *pr=x, Lsum=0, Rsum=0;
	 for (int i=0; i<R; ++i) {
	 	 int rSum=0, cSum=0, *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

char OutputFileName[bufSize];
FILE *openOutput(const char *typeName) {
//    ----------
  int sub=0; FILE *wfp=NULL; char baseName[bufSize], buf[bufSize];
  snprintf(baseName, bufSize, "./%dx%d%s", R, C, typeName);
  snprintf(buf, bufSize, "%s.txt", baseName);
  do {
    if (access(buf, F_OK)!=0) break;
    snprintf(buf, bufSize, "%s_%d.txt", baseName, ++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
//=========================================== 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) { const 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))) {
    const 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(); }
//   --------------

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

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

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];
  }
  return rectangleType==normalMagic ? violationMagic() : violationSemi();
} // 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;
    }
  }
  return violationSemi();
} // randomRestartAssoc

void doSwapMagic(int *X, const int i1, const int i2) {
//   -----------
	 const int tmp=xRectangle[i1]; xRectangle[i1]=xRectangle[i2]; xRectangle[i2]=tmp;
} // doSwapMagic

void doSwapAssoc(int *X, int i1, int i2) {
//   -----------
	int tmp=xRectangle[i1]; xRectangle[i1]=xRectangle[i2]; xRectangle[i2]=tmp; i1=RCm1-i1;
  if (i1!=i2) { i2=RCm1-i2; tmp=xRectangle[i1]; xRectangle[i1]=xRectangle[i2]; xRectangle[i2]=tmp; }
} // doSwapAssoc

int randomSwapAndCalcMagic(int *X, const int swaps) {
//   --------------------
  for (int swap=0; swap<swaps; ++swap) {
    const int a=freeIx[random_(RCfree)], b=freeIx[random_(RCfree)];
    doSwapMagic(X,a,b); updateViolationMagic(a, b);
  }
  return lastViolation;
} // randomSwapAndCalcMagic

int swapDeltaSemi(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(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) {
          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) {
          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) {
        int fi=freeIx[i]; bool ido=fullSearch;
        if (!ido) { 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) {
          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) {
		    doSwapMagic(X, freeIx[i1], freeIx[i2]); tabu[i1]=it+tabuLen;
        if (Nfree>5) tabu[i2]=it+tabuLen;
        if ((viol=updateViolationSemi(freeIx[i1],freeIx[i2]))==0) break;
      }
	   } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsSemi

int tabuSumsMagic(int *X, const int iterations, int *violation) { // R==C
//  -------------
  int it=0, viol=*violation;
  if (viol) {
    int zeroBD=0; const int noItGainX=Nfree<10 ? 1 : Nfree<20 ? 2 :
        Nfree<30 ? 3 : Nfree<40 ? 4 : Nfree<50 ? 5 : 6;
    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<R; ++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 (bestDelta==0) {
          if (++zeroBD>(2*tabuLen)) {
            viol=randomSwapAndCalcMagic(X,noItGainX); if (viol==0) break; zeroBD=0;
          } else {
            const int a=freeIx[i1], b=freeIx[i2]; doSwapMagic(X,a,b); viol=updateViolationMagic(a,b);
          }
        } else {
          const int a=freeIx[i1], b=freeIx[i2];
          doSwapMagic(X,a,b); viol=updateViolationMagic(a,b); if (viol==0) break; zeroBD=0;
        }   
        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) {
		    doSwapAssoc(X, freeIx[i1], freeIx[i2]); tabu[i1]=it+tabuLen; tabu[RCfree_1-i2]=it+tabuLen;
        if (Nfree>5) { tabu[i2]=it+tabuLen; tabu[RCfree_1-i1]=it+tabuLen; }
        if ((viol=updateViolationAssoc(freeIx[i1],freeIx[i2]))==0) break;
      }
	   } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsAssoc
//==================================================== main ================================================

typedef bool (*t_bNoParm)(); typedef int  (*t_iP)(int *); typedef int  (*t_iPcIiP)(int *, const int, int *);

struct t_type { const char *name; t_bNoParm initialize; t_iP randomRestart; t_iPcIiP  tabuSums; };

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

bool isPan(int *x) {
//   -----
  int *pL, *pR;
  for (int i=0; i<R; ++i) { // Pandiagonal sums
	  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: 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 (noZero&!(singlyEven|(rectangleType==Associative))) {
    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) {
//   ---------
  retry=F; fullSearch=T; printf("Retry with full search");
  if (useAssoc()) { *sqType=&rectangleTypes[Associative]; printf(" associative"); } printf(".\n");
} //

//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))
  if (getType()&&checkRC()) { // checkRC uses rectangleType
    struct t_type *sqType=&rectangleTypes[rectangleType];

    if (sqType->initialize()) {
      do { // while another
        const int runs=getInt("Runs? "); if (runs<=0) break;
        const int iterations=getInt("Iterations? "); if (iterations<=0) break;
        
        FILE *wfp=openOutput(getName());
        if (wfp!=NULL) {
          bool noRectangle=T; fullSearch=(RCfree<(75*RC/100));
          retry=noZero|(!fullSearch&(iterations>(400*N)));
          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) {
              int violation=sqType->randomRestart(xRectangle), clock1=clock();
              const int it=sqType->tabuSums(xRectangle, iterations, &violation);
              const double t=((clock()-clock1)/((double)CLOCKS_PER_SEC)); elapsedT+=t;

              if (violation==0) {
                //if (!checkRectangle(xRectangle)) break;
                ++good; goodIt+=it; goodET+=t; if (it>maxGoodIt) maxGoodIt=it;
                noRectangle=F; if (i==0) putchar('\n');
                printf("Run: %d, iterations: %d, time: %.1fs.\n", i+1, it, t);
                printf("%s is %s.\n\n", (R==C) ? "Square" : "Rectangle", checkType(xRectangle));
                if (!(consolePrint(xRectangle)&&filePrint(xRectangle, wfp))) {
                  reportError("Output failed"); break;
                } fflush(wfp);
              } else failET+=t;       
	          } // for(int i ...
            if (runs>0) {
              double avgT=elapsedT / runs;
              if (good==0) {
                printf("\nRuns: %d, %ss: 0, Average Time: %.1fs",
                  runs, (R==C) ? "Square" : "Rectangle", avgT);
              } else {
                printf("\nRuns: %d, %ss: %d, Success Iterations: Avg %d Max %d\n",
                        runs, (R==C) ? "Square" : "Rectangle", good, goodIt/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");
            }
            if (!(noRectangle&retry)) break; initRetry(&sqType);
          } while (T); // once and retry?
          fclose(wfp); if (noRectangle) remove(OutputFileName); else ok=T;
        } // if (wfp ...
        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