/*
*  File:    MagicRectangles.cpp
*  Author:  S Harry White
*  Created: 2013-12-02
*  Updated: 2022-01-11
*    Tidy code.
*  Updated: 2022-06-20
*    Add concentric rectangles.
*  Updated: 2022-06-24
*    Fix even order concentric rectangles.
*  Updated: 2022-08-11
*    Change name from MagicSquares.
*    Modify the user interface:
*      - input "how many" instead of "runs" and "iterations"
*      - change "another" to select a new order instead of a re-run of the same
*  Updated: 2022-08-14
*    Add pandiagonal for non-square rectangles.
*/

/*
*  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 magic rectangles.
*
*  The type of rectangle can be specified as:
*    semimagic, magic, associative, pandiagonal, or concentric
*
*  Note: If the rectangle is not square, magic equates to semimagic.
*/

#include "stdafx.h"
#include <conio.h>
#include <direct.h>
#include <errno.h>
#include <io.h>
#include <math.h>
#include <time.h>
#include <Windows.h>

int R, C,
B,             // bigger of R or C
N,             // sqrt(R*C)
Rm1, Cm1,      // R-1, C-1
Rp1, Cp1,      // R+1, C+1
Rd2, Cd2,      // R/2, C/2
RC,            // R *C
RCm1,          // RC-1
RCm2,          // RC-2
RCd2,          // RC/2
S2,            // RC+1
S1,            // S2/2 for odd R, C
SR,            // C*S2/2
SC,            // R*S2/2
SD;            // lesser of SR or SC

int *xRectangle, *numbers, *row, *col, *opp, *diagsL, *diagsR,
    *rowSum, *colSum, *diagsLsum, *diagsRsum, *tabu, diagLsum, diagRsum, tabuLen;
const bool F=false, T=true;
bool odd, *diagL, *diagR, allocOpp, allocDiags, allocPanDiags, sevn; // singly even

enum magicType { normalSemiMagic, normalMagic, Associative, nearAssociative, Concentric, Pandiagonal };
char *typeName[]={"Semimagic","Magic","Associative","NearAssociative","Concentric","Pandiagonal" };
const magicType firstType=normalSemiMagic, lastType=Pandiagonal; magicType rectangleType;

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(char *msg) { printf("%sError: %s.\n", bell ? "\a\n" : "", msg); return F; }
//   -----------
//=================================================== store ===================================================

void freeStore() {
//   ---------
  if (xRectangle!=NULL) { free(xRectangle); xRectangle=NULL; }
  if (numbers!=NULL)    { free(numbers);    numbers   =NULL; }
  if (row!=NULL)        { free(row);        row       =NULL; }
  if (col!=NULL)        { free(col);        col       =NULL; }
  if (opp!=NULL)        { free(opp);        opp       =NULL; }
  if (diagsL!=NULL)     { free(diagsL);     diagsL    =NULL; }
  if (diagsR!=NULL)     { free(diagsR);     diagsR    =NULL; }
  if (rowSum!=NULL)     { free(rowSum);     rowSum    =NULL; }
  if (colSum!=NULL)     { free(colSum);     colSum    =NULL; }
  if (tabu!=NULL)       { free(tabu);       tabu      =NULL; }
  if (diagL!=NULL)      { free(diagL);      diagL     =NULL; }
  if (diagR!=NULL)      { free(diagR);      diagR     =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; }

bool allocateStore() {
//   -------------
  bool ok=T;
  if (ok=allocateInts(&xRectangle, RC))
    if (ok=allocateInts(&numbers, RC))
      if (ok=allocateInts(&row, RC))
        if (ok=allocateInts(&col, RC))
          if (ok=allocateInts(&rowSum, R))
            if (ok=allocateInts(&colSum, C))
              ok=allocateInts(&tabu, RC);
  if ((ok)&(allocOpp)) ok=allocateInts(&opp, RC);
  if ((ok)&(allocDiags)) {
    if (ok=allocateBools(&diagL, RC))
      ok=allocateBools(&diagR, RC);
  }
  if ((ok)&(allocPanDiags)) {
    if (ok=allocateInts(&diagsL, RC))
      if (ok=allocateInts(&diagsR, RC))
        if (ok=allocateInts(&diagsLsum, B))
          ok=allocateInts(&diagsRsum, B);
  }
  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(char *prompt) {
//  ------
  printf(prompt); int n=0; scanf_s("%d", &n); clearLine(getchar());
  if (n<=0) reportError("Invalid input"); return n;
} // getInt

bool getInts(int *p, int *q, int c) {
//   -------
  bool ok=F; *p=0; *q=0; if (c<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

bool getYorInts(int *p, int *q) {
//   ----------
  bool ok=F; int c; *p=-1; do { c=getchar(); } while ((c==' ')|(c=='\t')|(c=='\n') );
  if ( (c=='Y')|(c=='y') ) ok=T; else if ( (c!='N')&(c!='n') ) return getInts(p, q, c);  
  clearLine(c); return ok;
} // getYorInts

char *getLine(char *s, 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 rectangle, (semimagic, magic, associative, "
    "concentric, pandiagonal), as s, m, a, c, or p? ");
  char buf[bufSize], *s=getLine(buf, '\n');
  if ((s!=NULL)&&(*s!='\0')) {
    while ((*s==' ')|(*s=='\t')) ++s;
    if (strlen(s)>=1) {
      sevn=((R&3)==2)&((C&3)==2); // singly even
      switch (s[0]) {
        case 'a': rectangleType=sevn?nearAssociative:Associative; return T;
        case 'c': rectangleType=Concentric; return T;
        case 'm': rectangleType=(R==C) ? normalMagic : normalSemiMagic; return T;
        case 'p': rectangleType=Pandiagonal; return T;
        case 's': rectangleType=normalSemiMagic; return T;
        default: break;
      }
    }
  }
  printf("\a\nUnknown type: \"%s\".\n", s); return F;
} // getType
//=============================================== output =================================================

bool openDir() {
//   -------
  char buf[bufSize], msg[bufSize]; const char *baseName="MagicRectangles";
  strcpy_s(buf, bufSize, baseName); int sub=0;

  do {
    if (_mkdir(buf)==-1) {
      if (errno!=EEXIST) {
        sprintf_s(msg, bufSize, "Can't make folder %s", buf); perror(msg); return F;
      }
      sprintf_s(buf, bufSize, "%s_%d", baseName, ++sub);
    } else break;
  } while (T);

  if (_chdir(buf)==-1) {
    sprintf_s(msg, bufSize, "Can't open folder %s", buf); perror(msg); return F;
  }
  printf("Output files are in folder %s\n\n", buf); return T;
} // openDir

char OutputFileName[bufSize];
FILE *openOutput() {
//    ----------
  char buf[bufSize], defaultName[bufSize]; int sub=0; FILE *wfp=NULL;
  magicType t=((R!=C)&(rectangleType==normalSemiMagic))?normalMagic:rectangleType;
  sprintf_s(buf, bufSize, "%s", typeName[t]);
  sprintf_s(defaultName, bufSize, "%dx%d%s", R, C, buf);
  sprintf_s(buf, bufSize, "%s.txt", defaultName);
  do {
    if (_access_s(buf, 00)==ENOENT) break;
    else sprintf_s(buf, bufSize, "%s_%d.txt", defaultName, ++sub);
  } while (T);
  if (fopen_s(&wfp, buf, "w")==0) {
    printf("\nOutput file is %s\n", buf); strcpy_s(OutputFileName, bufSize, buf);
  }
  else {
    char msg[bufSize]; sprintf_s(msg, bufSize, "\a\nCan't open for write %s", buf); perror(msg);
  }
  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); }
//===================================================== make ==================================================

void initializeDiags() {
//   ---------------
  if (R<=C)
    for (int i=0; i<C; ++i) {
      int k=i, l=i, r=0;
      while (r<R) { diagsL[k]=i; k+=((i+r)==Cm1) ? 1 : Cp1; diagsR[l]=i; l+=(i==r) ? C+Cm1 : Cm1; ++r; };
    }
  else
    for (int i=0; i<R; ++i) {
      int k=i*C, l=k, c=0;
      while (c<C) { diagsL[k]=i; k+=(i+c)==Rm1 ? Cp1-RC : Cp1; diagsR[l]=i; l+=(i==c) ? RC-Cm1 : -Cm1; ++c; };
    }
  //for (int i=0; i<RC; ++i) printf("%d %d\n", i, diagsL[i]); putchar('\n');
  //for (int i=0; i<RC; ++i) printf("%d %d\n", i, diagsR[i]);
} // initializeDiags

void initializeOpp() {
//   -------------
  for (int i=0; i<RC; ++i) opp[i]=RCm1-i;
  if (rectangleType==Concentric) {
    const int s=(R<C)?R:C, sub=odd?1:(R==C)?2:1, lim=s/2-sub;
    for (int i=0; i<RC; ++i) {
      const int r=row[i], c=col[i], ro=Rm1-r, co=Cm1-c;
      if ((r<lim)&(c>r)&(c<(Cm1-r))) { const int t=ro*C+c; opp[i]=t; opp[t]=i; }
      if ((c<lim)&(r>c)&(r<(Rm1-c))) { const int t=r*C+co; opp[i]=t; opp[t]=i; }
      if (!odd&(R!=C)) {
        if (R<C) { if ((r==lim)&(c>=r)&(c<=(Cm1-r))) { const int t=ro*C+c; opp[i]=t; opp[t]=i; }
        } else { if ((c==lim)&(r>=c)&(r<=(Rm1-c))) { const int t=r*C+co; opp[i]=t; opp[t]=i; } }
      }
    }
  } else { // Associative, nearAssociative
    if (sevn) { int i=(Rd2-1)*C; opp[i]=i+C; opp[i+C]=i; i=Rd2*C-1; opp[i]=i+C; opp[i+C]=i; }
  }
} // initializeOpp

bool initialize() {
//   ----------
  Rm1=R-1; Rp1=R+1; Rd2=R/2; Cm1=C-1; Cp1=C+1; Cd2=C/2; RC=R*C; RCd2=RC/2; RCm1=RC-1; RCm2=RCm1-1;
  B=R>C?R:C; S2=RC+1; S1=S2/2; SR=C*S2/2; SC=R*S2/2; SD=(R<C)?SC:SR; odd=R&1;
  N=(R==C) ? R : (int)sqrt((double)RC); tabuLen=(N==4) ? 5 : (N+N)/3;
  FW0=fieldWidth(); FW=FW0+1; if (FW>maxFieldWidth) return reportError("Output format overflow");
  magicType t=rectangleType; allocOpp=(t==Associative)|(t==nearAssociative)|(t==Concentric);
  allocDiags=(t==normalMagic); allocPanDiags=(t==Pandiagonal); if (!allocateStore()) return F;
  for (int i=0; i<RC; ++i) {
    const int idC=i/C, imC=i%C; numbers[i]=i+1; row[i]=idC; col[i]=imC;
    if (rectangleType==normalMagic) { diagL[i]=(idC==imC); diagR[i]=((idC+imC)==Cm1); }
  }
  if (allocOpp) { initializeOpp(); if (odd) xRectangle[RCd2]=S1; } if (allocPanDiags) initializeDiags();
  seed_rand(); return T;
} // initialize

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

bool pairSwapped;
int updateViolationSym(int i1, int i2) {
//  ------------------
  if (pairSwapped) {
    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=opp[i1];
    if (i1!=i2) {
      i2=opp[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);
  }
  return lastViolation;
} // updateViolationSym

int updateViolationPan(const int i1, const int i2) { // R==C
//  ------------------
  const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
            L1=diagsL[i1], R1=diagsR[i1], L2=diagsL[i2], R2=diagsR[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);
  }
  if (L1!=L2) {
    oldv+=abs(diagsLsum[L1]-SD)+abs(diagsLsum[L2]-SD); diagsLsum[L1]+=d; diagsLsum[L2]-=d;
    newv+=abs(diagsLsum[L1]-SD)+abs(diagsLsum[L2]-SD);
  }
  if (R1!=R2) {
    oldv+=abs(diagsRsum[R1]-SD)+abs(diagsRsum[R2]-SD); diagsRsum[R1]+=d; diagsRsum[R2]-=d;
    newv+=abs(diagsRsum[R1]-SD)+abs(diagsRsum[R2]-SD);
  }
  return lastViolation+=(newv-oldv);
} // updateViolationPan

int violationSemi() {
//  -------------
  int violation=0; int *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; int *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 violationPan() {
//  ------------
  int violation=0; int *pr, *pc, *pL, *pR;
  // Row, column, pandiagonal constraints
  if (R<=C) {
    pr=xRectangle;
    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, Lsum=0, Rsum=0; pc=pL=pR=xRectangle+i;
      for (int j=0; j<R; ++j) {
        cSum+=*pc; Lsum+=*pL; Rsum+=*pR;
        pc+=C; pL+=((i+j)==Cm1) ? 1 : Cp1; pR+=(i==j) ? C+Cm1 : Cm1;
      }
      colSum[i]=cSum; diagsLsum[i]=Lsum; diagsRsum[i]=Rsum;
      violation+=abs(cSum-SC)+abs(Lsum-SD)+abs(Rsum-SD);
    }
  } else {
     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);
    }
    for (int i=0; i<R; ++i) {
      int rSum=0, Lsum=0, Rsum=0; pr=pL=pR=xRectangle+i*C;
      for (int j=0; j<C; ++j) {
        rSum+=*pr++; Lsum+=*pL; Rsum+=*pR;
        pL+=((i+j)==Rm1) ? Cp1-RC : Cp1; pR+=(i==j) ? RC-Cm1 : -Cm1;
      }
      rowSum[i]=rSum; diagsLsum[i]=Lsum; diagsRsum[i]=Rsum;
      violation+=abs(rSum-SR)+abs(Lsum-SD)+abs(Rsum-SD);
    }
  }
  return lastViolation=violation;
} // violationPan

void randomFillMagic(int *X) {
//   ---------------
  for (int i=0; i<RCm1; ++i) {
    const int ri=i+random(RC-i);
    if (i!=ri) { const int v=numbers[i]; numbers[i]=numbers[ri]; numbers[ri]=v; }
  }
  for (int i=0; i<RC; ++i) X[i]=numbers[i];
} // randomFillMagic

void randomFillSym(int *X) {
//   -------------
  for (int i=0; i<(RCd2-1); ++i) {
    const int ri=i+random(RCd2-i), v=numbers[i]; numbers[i]=numbers[ri]; numbers[ri]=v;
  }
  int next=0;
  for (int i=0; i<RC; ++i) if (i<opp[i]) {
    int v=numbers[next++]; if (random(2)) v=S2-v; X[i]=v; X[opp[i]]=S2-v;
  }
} // randomFillSym

int restartSemi(int *X) { randomFillMagic(X); return violationSemi(); }
//  -----------

int restartMagic(int *X) { randomFillMagic(X); return violationMagic(); }
//  ------------

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

int restartPan(int *X) { randomFillMagic(X); return violationPan(); }
//  ----------

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

void doSwapSym(int *X, int i1, int i2) {
//   ---------
  pairSwapped=!(odd&((i1==RCd2)|(i2==RCd2)));
  if (pairSwapped) {
    int tmp=X[i1]; X[i1]=X[i2]; X[i2]=tmp;
    if ((i1=opp[i1])!=i2) { i2=opp[i2]; tmp=X[i1]; X[i1]=X[i2]; X[i2]=tmp; }
  }
} // doSwapSym

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=X[i1], x2=X[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=X[i1], x2=X[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 swapDeltaSym(const int *X, const int i1, const int i2) {
//  ------------
  const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2],
            x1=X[i1], x2=X[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;
} // swapDeltaSym

int swapDeltaPan(const int *X, const int i1, const int i2) { // R==C
//  ------------
  //printf("swapDeltaPan in\n");
  const int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2], L1=diagsL[i1], R1=diagsR[i1],
            L2=diagsL[i2], R2=diagsR[i2], x1=X[i1], x2=X[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);
  }
  if (L1!=L2) { // Different left diagonals
    x=diagsLsum[L1]-SD; delta+=abs(x+y)-abs(x); x=diagsLsum[L2]-SD; delta+=abs(x-y)-abs(x);
  }
  if (R1!=R2) { // Different right diagonals
    x=diagsRsum[R1]-SD; delta+=abs(x+y)-abs(x); x=diagsRsum[R2]-SD; delta+=abs(x-y)-abs(x);
  }
  //printf("swapDeltaPan out\n");
  return delta;
} // swapDeltaPan

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;
      for (int i=0; i<R; ++i) {
        const int q=rowSum[i]; if (q<minRSum) { minRSum=q; minR=i; } else 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; } else if (q>maxCSum) { maxCSum=q; maxC=i; }
      }
      for (int i=RCm2; i>=0; --i) if (tabu[i]<=it) {
        if ((row[i]==minR)||(row[i]==maxR)||(col[i]==minC)||(col[i]==maxC))
          for (int j=i+1; j<RC; ++j) if (tabu[j]<=it) {
            const int delta=swapDeltaSemi(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 ...
      ++it;
      if (i1>=0) {
        doSwapMagic(X, i1, i2); tabu[i1]=it+tabuLen; if (N>5) tabu[i2]=it+tabuLen;
        if ((viol=updateViolationSemi(i1, i2))==0) break;
      }
    } // 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;
      for (int i=0; i<N; ++i) {
        int q=rowSum[i]; if (q<minRSum) { minRSum=q; minR=i; } else if (q>maxRSum) { maxRSum=q; maxR=i; }
        q=colSum[i]; if (q<minCSum) { minCSum=q; minC=i; } else if (q>maxCSum) { maxCSum=q; maxC=i; }
      }
      const int mxnSum=max(minRSum, minCSum), mnxSum=min(maxRSum, maxCSum);
      const bool doL=(diagLsum<=mxnSum)|(diagLsum>=mnxSum), doR=(diagRsum<=mxnSum)|(diagRsum>=mnxSum);
      for (int i=RCm2; i>=0; --i) if (tabu[i]<=it) {
        if ((row[i]==minR)||(row[i]==maxR)||(col[i]==minC)||(col[i]==maxC)||(doL&&diagL[i])||(doR&&diagR[i]))
          for (int j=i+1; j<RC; ++j) if (tabu[j]<=it) {
            const int delta=swapDeltaMagic(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 ...
      ++it;
      if (i1>=0) {
        doSwapMagic(X, i1, i2); tabu[i1]=it+tabuLen; if (N>5) tabu[i2]=it+tabuLen;
        if ((viol=updateViolationMagic(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsMagic

int tabuSumsSym(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;
      for (int i=0; i<R; ++i) {
        const int q=rowSum[i]; if (q<minRSum) { minRSum=q; minR=i; } else 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; } else if (q>maxCSum) { maxCSum=q; maxC=i; }
      }
      for (int i=RCm2; i>=0; --i) if (tabu[i]<=it) {
        if ((row[i]==minR)||(row[i]==maxR)||(col[i]==minC)||(col[i]==maxC))
          for (int j=i+1; j<RC; ++j) if (tabu[j]<=it) {
            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 ...
      ++it;
      if (i1>=0) {
        doSwapSym(X, i1, i2); tabu[i1]=it+tabuLen; tabu[opp[i2]]=it+tabuLen;
        if (N>5) { tabu[opp[i1]]=it+tabuLen; tabu[i2]=it+tabuLen; }
        if ((viol=updateViolationSym(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsSym

int tabuSumsPan(int *X, const int iterations, int *violation) { // R==C
//  -----------
  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,
        minLd=-1, maxLd=-1, minRd=-1, maxRd=-1,
        minRSum =INT_MAX, maxRSum =0, minCSum =INT_MAX, maxCSum =0,
        minLdSum=INT_MAX, maxLdSum=0, minRdSum=INT_MAX, maxRdSum=0;
      for (int i=0; i<R; ++i) {
        const int q=rowSum[i]; if (q<minRSum) { minRSum=q; minR=i; } else 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; } else if (q>maxCSum) { maxCSum=q; maxC=i; }
      }
      for (int i=0; i<B; ++i) {
        int q=diagsLsum[i]; if (q<minLdSum) { minLdSum=q; minLd=i; } else if (q>maxLdSum) { maxLdSum=q; maxLd=i; }
            q=diagsRsum[i]; if (q<minRdSum) { minRdSum=q; minRd=i; } else if (q>maxRdSum) { maxRdSum=q; maxRd=i; }
      }
      for (int i=RCm2; i>=0; --i) if (tabu[i]<=it) {
        if ((row[i]==minR)||(row[i]==maxR)||(col[i]==minC)||(col[i]==maxC)||
            (diagsL[i]==minLd)||(diagsL[i]==maxLd)||(diagsR[i]==minRd)||(diagsR[i]==maxRd))
          for (int j=i+1; j<RC; ++j) if (tabu[j]<=it) {
            const int delta=swapDeltaPan(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 ...
      ++it;
      if (i1>=0) {
        doSwapMagic(X, i1, i2); tabu[i1]=it+tabuLen; if (N>5) tabu[i2]=it+tabuLen;
        if ((viol=updateViolationPan(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSumsPan
//================================================= main ================================================

typedef int(*t_iP)(int *); typedef int(*t_iPcIiP)(int *, const int, int *);
struct t_type { char *name; t_iP restart; t_iPcIiP tabuSums; }; const int numTypes=lastType-firstType+1;
t_type rectangleTypes[numTypes]=
{
  { typeName[normalSemiMagic], restartSemi,  tabuSumsSemi  },
  { typeName[normalMagic],     restartMagic, tabuSumsMagic },
  { typeName[Associative],     restartSym,   tabuSumsSym   },
  { typeName[nearAssociative], restartSym,   tabuSumsSym   },
  { typeName[Concentric],      restartSym,   tabuSumsSym   },
  { typeName[Pandiagonal],     restartPan,   tabuSumsPan   }
};

bool checkRC() {
//   -------
  if ((R<=0)|(C<=0)) return reportError("Invalid order"); 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)|(t==nearAssociative)|(t==Concentric)|(t==Pandiagonal)) {
      printf("\a\nThere are no %dx%d %s magic rectangles.\n",
        R, C, t==Concentric ? "concentric" : t==Pandiagonal ? "pandiagonal" : "associative"); return F;
    }
  }
  if ((R&1)!=(C&1)) return reportError("There are no odd, even magic rectangles");
  if (t==Pandiagonal) {
    if (sevn) return reportError("There are no singly even pandiagonal magic squares");
    if ((R==3)&(C==3)) return reportError("There are no order 3 pandiagonal magic squares");
  }
  if (t==Concentric) {
    if (R==C) {
      if (R==4) { printf("\a\nThere are no order 4 concentric magic squares.\n"); return F; }
    } else {
      if ((R==3)|(C==3)) { printf("\a\nThere are no %dx%d concentric magic rectangles.\n", R, C); return F; }
    }
  }
  return T;
} // checkRC

//bool checkRectangle(const 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;
  if (localtime_s(&local, &startTime)==0) {
    char dateTime[100];
    size_t slen=strftime(dateTime, 100, "%A %Y-%m-%d %X %Z\n\0", &local); printf("%s\n", dateTime);
  }
} // outputLocalTime

void printElapsedTime(FILE *wfp, const time_t startTime) {
//   ----------------
  const int et=(int)difftime(time(NULL), startTime);
  if (et>0) fprintf(wfp, "\nelapsed time %d:%02d:%02d\n", et/3600, et%3600/60, et%60);
} // printElapsedTime

int main() {
//  ----
  bool another=T, inputOrder=T, ok=F; //outputLocalTime();
  if (openDir()) {
    do { 
      if (inputOrder) { printf("\nOrder? "); if (!getInts(&R, &C, -1)) break; }
      if (getType()&&checkRC()&&initialize()) { // checkTC uses rectangleType
        const int howMany=getInt("How many? "); if (howMany<=0) break;
        FILE *wfp=openOutput();
        if (wfp!=NULL) {
          const bool p=rectangleType==Pandiagonal;
          const int iterations=p|(R&1)?10000:1000,
                runs=p?1000000:(howMany>1000)?howMany:1000;
          int good=0; time_t startTime=time(NULL);
          t_type *sqType=&rectangleTypes[rectangleType];
	        for (int i=0; i<runs; ++i) {
            int violation=sqType->restart(xRectangle);
            sqType->tabuSums(xRectangle, iterations, &violation);
            if (violation==0) {
              ok=T; printf("%d\n", ++good);
              if (!filePrint(xRectangle, wfp)) { reportError("Output failed"); break; } fflush(wfp);
              if (good==howMany) break;
            } // if (violation==0)
	        } // for (int i ...
          printf("Number of %ss: %d\n", (R==C) ? "square" : "rectangle", good);
          printElapsedTime(stdout, startTime);
          fclose(wfp); if (good==0) remove(OutputFileName);
        } // if (wfp ...
        freeStore();
      } // if (checkRC ....
      printf("\nAnother order? input y (yes), n (no) or the order: ");
      if (getYorInts(&R, &C)) inputOrder=(R<0); else break;
    } while (another);
  }
  printf("\nPress a key to close the console.");
  while (!_kbhit()) Sleep(250);	return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main