/*
 *  File:    Pandiagonal.cpp
 *  Author:  S Harry White
 *  Created: 2011-09-04
 *  Updated: 2020-04-23:
 *    Add makeMargossianEvenKnight and makeMargossianOddKnight.
 *    Add option "how many".
 *  Updated: 2020-06-12:
 *    Sort output in Frénicle standard form removing duplicates.
 *  Updated: 2020-06-26:
 *    Add makeEvenAl_Asfizari, makeEvenPlanck.
 *    Add input selection of method.
 *  Updated: 2020-11-24:
 *    Adjust startSize and startSquaresSize.
 *  Updated 2022-01-20
 *    Replace & with && where necessary.
 *  Updated 2023-02-06
 *    Randomly select triads to reorder for makeMargossianOddKnight.
 *  Updated 2023-03-02
 *    Add makeCampbell for square number orders: 4, 9, 16, 25, ..
 *  Updated 2023-03-23
 *    Change makeCampbell to support all doubly-even and non-prime odd orders.
 *  Updated 2023-03-24
 *    Change makeCampbell to makeWithTranspose.
 */

/*
 *  Makes pandiagonal magic squares of order n.
 *  For singly-even n, the squares are near-pandiagonal.
 */

#include "stdafx.h"
#include <assert.h>
#include <conio.h>
#include <direct.h>
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <Windows.h>

typedef unsigned int Uint; Uint MagicConstant; struct tdigits { unsigned short int a, b; };
const int startSize=8, startSquaresSize=16;
int **xSquare=NULL, **tSquare=NULL, *Permuted=NULL, **Squares=NULL, *square=NULL,
    allocatedSize, allocatedSquaresSize, Sindex, sLen;
const bool F=false, T=true; bool *numberUsed=NULL, ProgramError; 
//================================================ store =================================================

char *storeAllocFail="Storage allocation failed"; bool bell=T;
bool reportError(char *msg) {
//   -----------
  printf("%sError: %s.\n", bell ? "\a\n" : "",  msg); return bell=F;
} // reportError

void freeInts(int **line) { if (*line!=NULL) { free(*line); *line=NULL; }}
//   -------

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

void freeStore() {
//   ---------
  freeSquare(&xSquare, allocatedSize); freeSquare(&tSquare, allocatedSize);
  freeSquare(&Squares, allocatedSquaresSize); freeInts(&square); freeInts(&Permuted);
  if (numberUsed!=NULL) { free(numberUsed); numberUsed=NULL; } allocatedSize=0; allocatedSquaresSize=0;
} // freeStore

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

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

bool allocateSquares(const int size) {
//   ---------------
  bool ok; int length=size*size;  // Have to enter full square to accomodate near-associative.
  Squares=(int**) malloc(startSquaresSize*sizeof(int*)); ok=(Squares!=NULL);
  if (ok) {
    int numAllocated=startSquaresSize;
    for(int i=0; i<startSquaresSize; i++) {
      int *p=(int*) malloc(length*sizeof(int)); (Squares)[i]=p; if (p==NULL) { numAllocated=i; ok=F; break; }
    }
    if (!ok) freeSquare(&Squares, numAllocated);
  }
  return ok;
} // allocateSquares

bool increaseSquares() {
//   ---------------
  bool ok; int length=allocatedSize*allocatedSize, // of square
    size=allocatedSquaresSize+allocatedSquaresSize, **tmp=(int**) malloc(size*sizeof(int*));
  if (ok=(tmp!=NULL)) {
    for (int i=0; i<allocatedSquaresSize; i++) tmp[i]=Squares[i];
    free(Squares); Squares=tmp; int numAllocated=size;
    for(int i=allocatedSquaresSize; i<size; i++) {
      int *p=(int*) malloc(length*sizeof(int)); (Squares)[i]=p; if (p==NULL) { numAllocated=i; ok=F; break; }
    }
    if (!ok) freeSquare(&Squares, numAllocated);
  }
  if (ok) allocatedSquaresSize=size; else reportError(storeAllocFail);
  return ok;
} // increaseSquares

bool allocateStore(int size) {
//   -------------
  bool ok=T; if (size<startSize) size=startSize;
  if (size>allocatedSize) { freeStore();
    if (ok=allocateSquare(&xSquare, size))
      if (ok=allocateSquare(&tSquare, size))
        if (ok=allocateSquares(size)) {
          allocatedSize=size; allocatedSquaresSize=startSquaresSize;
          if (ok=allocateInts(&Permuted, size))
            if (ok=allocateInts(&square, size*size))
              ok=(numberUsed=(bool *) malloc((size*size+1)*sizeof(bool)))!=NULL;
        }
    if (!ok) { freeStore(); reportError(storeAllocFail); }
  }
  return ok;
} //allocateStore
//============================================== input ============================================

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

bool getYorOrder(int *n) { // 'y' or 'n' or the order
//   -----------
  bool result=F; int c; *n=0; do { c=getchar(); } while ((c==' ')|(c=='\t')|(c=='\n') );
  if ( (c=='Y')|(c=='y') ) result=T;
  else if ( (c!='N')&(c!='n') )
    if ( ('1'<=c)&(c<='9') ) {
      int i=c-'0'; while ( ('0'<=(c=getchar()))&&(c<='9') ) i=i*10+c-'0'; *n=i; result=T;
    }   
  clearLine(c); return result;
} // getYorOrder

int getInt() { int n=0; scanf_s("%d", &n); clearLine(getchar()); return n; }
//  ------

bool isPrime(const int n) {
//   -------
  if (n&1) {
    bool prime=(n!=1); if (n<9) return prime; const int r=(int)sqrt((float)n);
    for (int f=3; f<=r; f+=2) if ((n%f)==0) { prime=F; break; } return prime;
  }
  return n==2;
} // isPrime

int inputKind(const int n) {
//  ---------
  int k=0;
  if (n&1) {
    if (n==1) k=3;
    else {
      if (isPrime(n)) {
        printf("Method, ( 1 Margossian I, 2 Siamese, 3 from SODLS)? ");
        k=getInt(); if (k<1) k=1; else if (k==2) k=3; else if (k>2) k=4; 
      } else {
        if ((n%3)==0) {
          printf("Method, ( 1 Margossian I, 2 Margossian II, 3 transpose)? ");
          k=getInt(); if (k<1) k=1; else if (k>=3) k=5;
        } else {
          printf("Method, ( 1 Margossian I, 2 Siamese, 3 from SODLS, 4 transpose)? ");
          k=getInt(); if (k<1) k=1; else if (k>=2) { if (k<=4) ++k; else k=5; }
        }
      }
    }
  } else {
    if ((n&3)==0) {
      printf("Method, ( 1 Al_Asfizari I, 2 Al_Asfizari II, 3 Planck, 4 Margossian IV, "
             "5 Margossian V, 6 transpose)? ");
      k=getInt(); if (k<1) k=1; else if (k>6) k=6;
    }
  }
  return k;
} // inputKind
//================================================= output ================================================

const int bufSize=1024;
bool openDir() {
//   -------
  int sub=0; char buf[bufSize], msg[bufSize]; const char *baseName="PandiagonalSquares";
  strcpy_s(buf, bufSize, baseName);
  do {
    if (_mkdir(buf)==0) break;
    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);
  } while (T);
  if (_chdir(buf)!=0) { 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

FILE *open_File(const int n) {
//    ---------
  FILE *wfp=NULL; char buf[bufSize]; sprintf_s(buf, bufSize, "Order%d.txt", n);
  if (fopen_s(&wfp, buf, "a")!=0) {
    char msg[bufSize]; sprintf_s(msg, bufSize, "\a\nCan't open for write %s", buf); perror(msg);
  }
  return wfp; 
} // open_File

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

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;

//bool printSquare(int **x, const int n, FILE *wfp) {
////   -----------
//  const int fw0=fieldWidth(n), fw=fw0+1; if (fw>maxFieldWidth) return F;
//  for (int i=0; i<n; i++) {
//    if (!fprintFW[fw0](wfp, x[i][0])) return F;
//    for (int j=1; j<n; j++) if (!fprintFW[fw](wfp, x[i][j])) return F;
//    if (fputc('\n', wfp)==EOF) return F;
//  }
//  return fputc('\n', wfp)!=EOF;
//} // printSquare

bool outputSquares(const int n, FILE *wfp) {
//   -------------
  const int m=n-1, fw0=fieldWidth(n), fw=fw0+1; if (fw>maxFieldWidth) return F;
  for (int i=0; i<Sindex; i++) {
    int *x=Squares[i], j=0, r=n, c;
    while (r--) {
      if (!fprintFW[fw0](wfp, x[j++])) return F; c=m; while (c--) if (!fprintFW[fw](wfp, x[j++])) return F;
      if (fputc('\n', wfp)==EOF) return F;
    }
    if (fputc('\n', wfp)==EOF) return F;
  }
  return T;
} // outputSquares
//============================================= check squares ===================================================

bool isCorrect(int **x, const int n) {
//   ---------
  const Uint chkSum=MagicConstant, nn=n*n; Uint sumX, sumY, sumXY=0, sumYX=0;
  for (Uint i=1; i<=nn; i++) numberUsed[i]=F;
  for (int i=0; i<n; i++) {
    sumX=0; sumY=0; for (int j=0; j<n; j++) { sumX+=x[i][j]; sumY+=x[j][i]; numberUsed[x[j][i]]=T; }
    if ((sumX!=chkSum)|(sumY!=chkSum)) { return F; } sumXY+=x[i][n-i-1]; sumYX+=x[i][i];
  }
  for (Uint i=1; i<=nn; i++) if (!numberUsed[i]) return F; return ((sumXY==chkSum)&(sumYX==chkSum));
} // isCorrect

bool isPand(int **x, const int n) {
//   ------
  const int m=n-1, maxCount=n&1 ? 0 : (n&3)==0 ? 0 : 4; int count=0;
  // Main diagonals already checked in isCorrect.
  for (int i=1; i<n; ++i) {
    int sumYX=0, sumXY=0;
    for (int r=0, c=i; r<n; ++r, ++c) { sumYX+=x[r][c%n]; sumXY+=x[r][m-c%n]; }
    if (sumYX!=MagicConstant) ++count; if (sumXY!=MagicConstant) ++count;
  }
  return count==maxCount;
} // isPand
//================================= sort Frénicle ascending =================================

int cmpSquares(int *key, int *s) {
//  ----------
  int i=0; while(i++<sLen) if (*key<*s) return -1; else if (*key++>*s++) return 1; return 0;
} // cmpRowPerms

int findSquare(int *key, int *insertPoint) {
//  ----------
  int first=0, last=Sindex-1, middle;
  while (first<=last) {
    middle=(first+last)/2; int cmp=cmpSquares(key,Squares[middle]);
    if (cmp<0) last=middle-1; else if (cmp>0) first=middle+1; else return T;
  }
  *insertPoint=first; return F;
} // findSquare

bool insertSquare(int *sq, const int insertPoint, bool *storeFail) {
//   ------------
  if (Sindex==allocatedSquaresSize) if (!increaseSquares()) { *storeFail=T; return F; }
  int *p=Squares[Sindex]; for (int i=Sindex; i>insertPoint; --i) Squares[i]=Squares[i-1];
  for (int i=0; i<sLen; ++i) p[i]=sq[i]; Squares[insertPoint]=p; ++Sindex; return T;
} // insertSquare

void pushSquare(int **x, const int n, bool *storeFail) {
//   ----------
  int i=0, insertPoint; bool push=T;
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) square[i++]=x[r][c];
  if (Sindex==0) insertPoint=0; else if (findSquare(square, &insertPoint)) push=F;
  if (push) insertSquare(square, insertPoint, storeFail);
} // pushSquare

#define fForm(iX, jX) for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) t[iX][jX]=x[i][j]
#define xCopy() for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) x[i][j]=t[i][j]

void putInFrenicleForm(int **x, const int n) {
//   -----------------
  const int m=n-1, l=n-2; bool inFrenicleForm=F;
  int **t=tSquare, minC=min(min(x[0][0], x[0][m]), min(x[m][0], x[m][m]));
  if (x[0][0]==minC) {
    if (x[0][1]<x[1][0]) inFrenicleForm=T;
    else fForm(j,i);                      // YX
  } else if (x[0][m]==minC) {
    if (x[1][m]<x[0][l]) fForm(m-j,i);    // -90
    else fForm(i,m-j);                    // Y   
  } else if (x[m][m]==minC) {
    if (x[m][l]<x[l][m]) fForm(m-i,m-j);  // 180
    else fForm(m-j,m-i);                  // XY
  } else { // x[m][0]==minC
    if (x[l][0]<x[m][1]) fForm(j,m-i);    // 90
    else fForm(m-i,j);                    // X
  }
  if (!inFrenicleForm) xCopy();
} // putInFrenicleForm

bool putSquare(int **x, const int n) {
//   ---------
  bool storeFail=F; putInFrenicleForm(x,n); pushSquare(x, n, &storeFail); return !storeFail;
} // putSquare
//===================================== common procedures =========================================

void seed_rand() { srand((unsigned int)time(NULL)); }
int random(const int x) { return rand()%x; } bool randomBool() { return random(2)==0; }

void swapRows(int **Q, const int r1, const int r2) { int *t=Q[r1]; Q[r1]=Q[r2]; Q[r2]=t; }
//   --------

void swapCols(int **Q, const int n, const int c1, const int c2) {
//   --------
  for (int r=0; r<n; ++r) { const int t=Q[r][c1]; Q[r][c1]=Q[r][c2]; Q[r][c2]=t; }
} // swapCols

void reverseRows(int **x, const int n, const int a) {
//   -----------
  const int m=n-1, b=m-a, g=n/2;
  for (int c=0; c<g; ++c) {
    int t=x[a][c]; x[a][c]=x[a][m-c]; x[a][m-c]=t; t=x[b][c]; x[b][c]=x[b][m-c]; x[b][m-c]=t;
  }
} // reverseRows

void reverseCols(int **x, const int n, const int a) {
//   -----------
  const int m=n-1, b=m-a, g=n/2;
  for (int r=0; r<g; ++r) {
    int t=x[r][a]; x[r][a]=x[m-r][a]; x[m-r][a]=t; t=x[r][b]; x[r][b]=x[m-r][b]; x[m-r][b]=t;
  }
} // reverseCols

bool notMainDiags(int **x, const int n) {
//   ------------
  const Uint chkSum=MagicConstant; Uint sumXY=0, sumYX=0;
  for (int i=0; i<n; i++) { sumXY+=x[i][n-i-1]; sumYX+=x[i][i]; } return (sumXY!=chkSum)|(sumYX!=chkSum);
} // notMainDiags

void panForm(int **x, int **t, const int n) {
//   -------
  const int m=n-1, numRows=random(n); int loop=2, numCols=random(n), num=numRows;
  while(loop--) {
    if (num!=0) {
      for (int r=0; r<n; ++r) Permuted[r]=(r+num)%n; int *saveP=x[0]; int saveR=0, next=0;
      for (int r=0; r<m; ++r) {
        int tmp=Permuted[next];
        if (tmp==saveR) { x[next]=saveP; saveP=x[next=++saveR]; } else { x[next]=x[tmp]; next=tmp; }
      } x[next]=saveP;
    }
    if (numCols!=0) { // rotate 90
      for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) { t[c][m-r]=x[r][c]; }
      for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]=t[r][c]; } num=numCols;
    if ((loop==0)&((n&3)==2)&&notMainDiags(x,n)) { loop=1; numCols=0; num=1; }
  } // while(loop--)
} // panForm

void xFormA_D(int **x, const int n) {
//   --------
  const int g=n/2, m=n-1, k=g+n/4, t=m+g;
  for (int r=g; r<k; ++r) swapRows(x, r, t-r); for (int c=g; c<k; ++c) swapCols(x, n, c, t-c);
} // xFormA_D

bool getEvenFactor(const int n, int *nA) {
//   -------------
  const int r=(int)sqrt((float)n); int num=0, *fac=NULL; bool ok=F;
  if (allocateInts(&fac,n)) {
    int f, num=2; fac[0]=1; fac[1]=n;
    for (int a=2; a<=r; ++a)
      if ((n%a)==0) { const int b=n/a; if (((a&3)==0)|((b&3)==0)) { fac[num++]=a; fac[num++]=b; }}
    f=fac[random(num)]; if ((f&3)!=0) f=n/f; *nA=f; ok=T; freeInts(&fac);
  }
  return ok; 
} // getEvenFactor

//============================================== transpose =================================================
//------------------------------------------ get sub-rectangle --------------------------------------------

int R, C,
    RC,   // R*C
    N,    // sqrt(RC)
    RCm1, // RC-1
    RCm2, // RC-2
    S2,   // RC-1
    SR,   // C*S2/2
    SC,   // R*S2/2
   *xRectangle=NULL, *numbers=NULL, *row=NULL, *rowSum=NULL, *col=NULL, *colSum=NULL, *tabu=NULL, tabuLen;

void freeStoreRect() {
//   -------------
  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 (rowSum!=NULL)     { free(rowSum);     rowSum    =NULL; }
  if (colSum!=NULL)     { free(colSum);     colSum    =NULL; }
  if (tabu!=NULL)       { free(tabu);       tabu      =NULL; }
} // freeStoreRect

bool allocateStoreRect() {
//   -----------------
  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) { freeStoreRect(); reportError("Storage allocation failed"); }
  return ok;
} // allocateStoreRect

bool initializeRect(const int r, const int c) {
//   --------------
  R=r; C=c; RC=R*C; RCm1=RC-1; RCm2=RCm1-1; S2=RCm1; SR=C*S2/2; SC=R*S2/2;
  N=(R==C) ? R : (int)sqrt((double)RC); tabuLen=(N==4) ? 5 : (N+N)/3;
  if (!allocateStoreRect()) return F;
  for (int i=0; i<RC; ++i) { const int idC=i/C, imC=i%C; numbers[i]=i; row[i]=idC; col[i]=imC; }
  return T;
} // initializeRect

int getFactor(const int n) {
//  ---------
  int num=0, *fac=Permuted;
  if (n&1) {
    int a=3, b=n/a; while (a<=b) { if (a*b==n) fac[num++]=a; a+=2; b=n/a; }
  } else {
    int a=2, b=n/a; while (a<=b) { if ((a*b==n)&&(!(b&1))) fac[num++]=a; a+=2; b=n/a; }
  }
  //return (num>0) ? fac[random(num)] : 1;
  return (num>0) ? (n>256)? fac[num-1]: fac[random(num)] : 1;
} // getFactor

int lastViolation; bool chkC=F;
int updateViolation(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 (chkC&(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);
} // updateViolation

int violation() {
//  ---------
  int violation=0; int *pr=xRectangle;
  for (int i=0; i<R; ++i) { // Row constraints
    int rSum=0;
    for (int j=0; j<C; ++j) rSum+=*pr++; rowSum[i]=rSum; violation+=abs(rSum-SR);
  }
  if (chkC) for (int i=0; i<C; ++i) { // Col constraints
    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;
} // violation

void randomFill(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];
} // randomFill

int restart(int *X) { randomFill(X); return violation(); }
//  -------

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

int swapDelta(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 (chkC&(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;
} // swapDelta

int tabuSums(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; }
      }
      if (chkC) 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=swapDelta(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) {
        doSwap(X, i1, i2); tabu[i1]=it+tabuLen; if (N>5) tabu[i2]=it+tabuLen;
        if ((viol=updateViolation(i1, i2))==0) break;
      }
    } // while (it ...
  } // if (viol ...
  *violation=viol; return it;
} // tabuSums

bool makeRectangle(const int a, const int b) {
//   -------------
  bool ok=F;
  int loop=1000, iterations=a*b; 
  while (loop--) {
    int violation=restart(xRectangle);
    tabuSums(xRectangle, iterations, &violation); if (violation==0) { ok=T; break; }  
  }
  return ok;
} // makeRectangle

bool getSub(int **x, const int a, const int b) {
//   ------
  if (!makeRectangle(a, b)) { ProgramError=T; return F; }
  int i=0; for (int r=0; r<a; ++r) for (int c=0; c<b; ++c) x[r][c]=xRectangle[i++];
  return T; 
} // getSub

bool chkDiags(int **x, const int n) {
//   --------
  const int m=n-1; bool *u=numberUsed;
  for (int i=0; i<n; ++i) u[i]=F;
  for (int r=0; r<n; ++r) { int v=x[r][r]; if (u[v]) return F; u[v]=T; }
  for (int i=0; i<n; ++i) u[i]=F;
  for (int r=0; r<n; ++r) { int v=x[r][m-r]; if (u[v]) return F; u[v]=T; }
  return T;
} // chkDiags

bool fillWithShifts(int **x, const int n, const int a, const int b) {
//   --------------
  const int m=n-1; int s=1, loop=(b>=9)?9:b; bool ok=F;
  while (loop--) {
    for (int r=0; r<a; ++r) for (int c=b; c<n; ++c) x[r][c]=x[r][c-b];
    for (int r=a; r<n; ++r) {
      for (int c=0; c<s; ++c) x[r][c]=x[r-a][n+c-s]; for (int c=s; c<n; ++c) x[r][c]=x[r-a][c-s];
    }
    if (chkDiags(x,n)) { ok=T; break; } s+=2;
  }
  //if (ok) printf("factor %d shift %d\n", a, s);
  return ok;
} // fillWithShifts

void transpose(int **x, int **t, const int n) {
//   ---------
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) t[c][r]=x[r][c];
} // transpose

bool getWithTranspose(int **x, int **t, const int n, const int a, const int b) {
//   ----------------
  bool ok=F;
  if (getSub(x,a,b)) { if (fillWithShifts(x,n,a,b)) { transpose(x,t,n); ok=T; } }
  if (!ok) ProgramError=T; return ok;
} // getWithTranspose

void makeWithTranspose(int **x, const int n) {
//   -----------------
  const int a=getFactor(n), b=n/a; bool ok=F;
  if (a==1) { ProgramError=T; }
  else {
    if (initializeRect(a,b)) {
      int loop=2; chkC=F;
      while (loop--) {
        if (getWithTranspose(x,tSquare,n,a,b)) {
          for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]+=(n*tSquare[r][c]+1);
          if (isCorrect(x, n)) { ok=T; break; }
        }
        chkC=T;
      }
    }
    freeStoreRect(); ProgramError=!ok;
  }
} // makeWithTranspose
//============================================ doubly evan ================================================

void makeAl_Asfizari(int **x, const int n) {
//   ---------------
  const int k=n/4, nmk=n-k; int i=1, j=n*n;
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) {
    const bool cmid=(c>=k)&(c<nmk); if ((r<k)|(r>=nmk)) x[r][c]=cmid?j:i; else x[r][c]=cmid?i:j; ++i; --j;
  }
  xFormA_D(x, n);
} // makeAl_Asfizari

void makeAl_Asfizari_of_Subsquares(int **x, int **t, const int n, const int a) {
//   -----------------------------
  int i=1, j=n*n;
  for (int r=0; r<n; ++r) {
    const int rma=r%a;
    for (int c=0; c<n; ++c) { if (t[rma][c%a]==0) { x[r][c]=j--; i++; } else { x[r][c]=i++; j--; }}
  }
  xFormA_D(x, n);
} // makeAl_Asfizari_of_Subsquares

void makeEvenAl_Asfizari(int **x, const int n, const bool nSquare) {
//   -------------------
  if ((n==4)|nSquare) makeAl_Asfizari(x, n);
  else {
    int f;
    if (getEvenFactor(n, &f)) {
      const int k=f/4, fmk=f-k; int **t=tSquare;
      for (int r=0; r<f; ++r) for (int c=0; c<f; ++c) {
        const bool cmid=(c>=k)&(c<fmk); if ((r<k)|(r>=fmk)) t[r][c]=cmid?1:0; else t[r][c]=cmid?0:1;
      }
      makeAl_Asfizari_of_Subsquares(x, t, n, f);
    }
  }
} // makeEvenAl_Asfizari

void makeMargossian4(int **x, const int n) {
//   ---------------
  const int g=n/2, m=n-1, half=n/2, j=n+half-1; 
  for (int i=0; i<half; ++i) Permuted[i]=i; for (int i=half; i<n; ++i) Permuted[i]=j-i;
  struct tdigits **q=(struct tdigits **)x;
  for (int i=0; i<n; ++i) {
    const int beta=(half*i)%n; q[m-i][0].a=beta; q[m-i][0].b=i; q[m][i].a=i; q[m][i].b=beta;
  }

  for (int r=0; r<m; ++r) {
    struct tdigits t=q[r][0]; const int ra=t.a, rb=t.b;
    for (int c=1; c<n; ++c) { t=q[m][c]; q[r][c].a=(ra+t.a)%n; q[r][c].b=(rb+t.b)%n; }
  }
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) {
    struct tdigits t=q[r][c]; x[r][c]=Permuted[t.a]*n+Permuted[t.b]+1;
  }
} // makeMargossian4

void makeGroupEven(const int n) {
//   -------------
  int **t=tSquare; const int m=n-1, g=n/2, n3d4=3*n/4;
  for (int r=0, i=0; r<n; r+=2) {
    for (int c=0; c<g; ++c) t[r][c]=++i; for (int c=0; c<g; ++c) t[r+1][c]=++i;
    for (int c=m; c>=g; --c) t[r][c]=++i; for (int c=m; c>=g; --c) t[r+1][c]=++i;
  }
  for (int r=g, i=0; r<n3d4; ++r,++i) swapRows(t, r, m-i);
} // makeGroupEven

void makeMargossianEvenKnight(int **x, const int n) {
//   ------------------------
  makeGroupEven(n); int **t=tSquare; const int m=n-1;
  do {
    const int a=random(m)+1, b=random(m)+1, a1=random(m)+1, b1=random(m)+1;
    int sum=0, R=0, C=0, R0=0, C0=0; if ((a==a1)|(b==b1)|((a+a1)==n)|((b+b1)==n)) continue;
    // Check row 0.
    for (int c=0; c<n; ++c) { x[0][c]=t[R][C]; sum+=t[R][C]; R=(R+b)%n; C=(C+a)%n; }
    if (sum!=MagicConstant) continue; R=(R0+b1)%n; C=(C0+a1)%n;
    for (int r=1; r<n; ++r) {
      R0=R, C0=C; for (int c=0; c<n; ++c) { x[r][c]=t[R][C]; R=(R+b)%n; C=(C+a)%n; } R=(R0+b1)%n; C=(C0+a1)%n;
    }
    if (isCorrect(x, n)&&isPand(x, n)) break;
  } while (T);
} // makeMargossianEvenKnight

void makeEvenPlanck(int **x, const int n) {
//   --------------
  const int m=n-1, k=n/4, g=n/2; int i=1, loop=k; bool *u=numberUsed, swapR=randomBool(), swapC=randomBool();
  for (int r=0; r<g; ++r) u[r]=F; for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]=i++;
  while (loop--) {
    int row=random(g); while (u[row]) row=random(g); u[row]=T; reverseRows(x, n, row);
    if (swapR) swapRows(x, row, m-row);
  }
  loop=k; for (int c=0; c<g; ++c) u[c]=F;
  while (loop--) {
    int col=random(g); while (u[col]) col=random(g); u[col]=T;
    reverseCols(x, n, col); if (swapC) swapRows(x, col, m-col);
  }
  xFormA_D(x, n);
} // makeEvenPlanck

void makeDEven(int **x, const int n, const int kind) {
//   ---------
  switch (kind) {
    case 1: makeEvenAl_Asfizari(x, n, T);   break; // n square pattern
    case 2: makeEvenAl_Asfizari(x, n, F);   break; // n square and subsquare patterns
    case 3: makeEvenPlanck(x, n);           break;
    case 4: makeMargossian4(x, n);          break;
    case 5: makeMargossianEvenKnight(x, n); break;
    case 6: makeWithTranspose(x, n);        break;
    default: break;
  }
}  // makeDEven
//=================================== singly evan =====================================

// [r][c], [r][c+1], [r+1][c], [r+1][c+1]
enum orders { Lo, Uo, Xo, Lp, Up, Xp, Uq };
int Order2x2[7][4] =
{
  { 3, 0, 1, 2 },  // 0  L
  { 0, 3, 1, 2 },  // 1  U
  { 0, 3, 2, 1 },  // 2  X
  { 1, 2, 3, 0 },  // 3  L'
  { 1, 2, 0, 3 },  // 4  U'
  { 2, 1, 0, 3 },  // 5  X'
  { 2, 1, 3, 0 }   // 6  U", U` mirror
};

int get2x2Order(int m, int r, int c) {
//  -----------
  int x=-1; const int f=m/2; r/=2; c/=2;
  if (r<f) { // L.LLL.L
    if (c==f) { x=r==0 ? ((m&3)==1) ? Uo : Lo : (r&1)==0 ? Lo : r==(f-1) ? Lp : Xp; }
    else { if (c>f) ++c; x=((r+c)&1)==0 ? r==(f-1) ? Up : Xp : Lo; }
  } else if (r==f) { // LL.U.LL
    x=c==f ? (m&3)==1 ? Uq : Uo : c==(m-1) ? Lp : (c<f)^((c&1)==0) ? Lp : Lo;
  } else { // L.LUL.L or X.XXX.X
    if (c==f) { x=r==(m-1) ? ((m&3)==1) ? Up : Lp : (r&1)==0 ? Lp : r==(f+1) ? Lo : Xo; }
    else { if (c>f) ++c; x=((r+c)&1)==0 ? Lp : r==(f+1) ? Uo : Xo; }
  }
  return x;
} // get2x2Order

void makeSEven(int **x, int **t, const int n) {
//   ---------
  const int g=n/2, which=randomBool();
  for (int r=0; r<g; ++r) {
    const int rr=r+r;
    for (int c=0; c<g; ++c) {
      int i=rr+c+1, j=which ? rr-c : r+c-(g-1)/2; if (i<0) i+=g; if (j<0) j+=g; x[rr][c+c]=(i%g+(j%g)*g)*4+1;
    }
  }
  for (int r=0; r<n; r+=2) for (int c=0; c<n; c+=2) {
    const int lo=x[r][c]; int *p=Order2x2[get2x2Order(g,r,c)];
    x[r][c]=lo+p[0]; x[r][c+1]=lo+p[1]; x[r+1][c]=lo+p[2]; x[r+1][c+1]=lo+p[3];
  }
  xFormA_D(x, n);
}  // makeSEven
//================================================== odd ======================================================

void makeSiamese(int **x, const int n) {
//   -----------
  const int k=random(n), l=random(n);
  for (int r=0; r<n; ++r) {
    const int rr=r+r;
    for (int c=0; c<n; ++c) { int i=rr+c+k, j=rr-c+l; if (i<0) i+=n; if (j<0) j+=n; x[r][c]=i%n+(j%n)*n+1; }
  }
} // makeSiamese

void makeMargossian3(int **x, const int n) {
//   ---------------
  const int m=n-1, e=n/3, numNoChange=(e-1)/2; int i=0, j=0, count=0, length=e;
  /* Use x[0] temporarily for 0, 1, ..., e-1 */ int *p=x[0]; for (i=0; i<e; ++i) p[i]=i;
  // Leave numNoChange triads unchanged.
  while (count++<numNoChange) {
    i=random(length); j=3*p[i];
    Permuted[j]=j; Permuted[j+1]=j+1; Permuted[j+2]=j+2; for (j=i; j<length; ++j) p[j]=p[j+1]; --length;
  }
  // Permute 2 triads.
  i=random(length); j=3*p[i]; Permuted[j]=j+1; Permuted[j+1]=j+2; Permuted[j+2]=j;
  for (j=i; j<length; ++j) p[j]=p[j+1]; --length; i=random(length); j=3*p[i];
  Permuted[j]=j+2; Permuted[j+1]=j; Permuted[j+2]=j+1; for (j=i; j<length; ++j) p[j]=p[j+1]; --length;
  // Reverse remaining triads.
  for (i=0; i<length; ++i) { j=3*p[i]; Permuted[j]=j+2; Permuted[j+1]=j+1; Permuted[j+2]=j; }
  struct tdigits **q=(struct tdigits **)x;
  for (int i=0; i<n; ++i) {
    const int iii=(3*i)%n; q[m-i][0].a=iii; q[m-i][0].b=i; q[m][i].a=i; q[m][i].b=iii;
  }
  for (int r=0; r<m; ++r) {
    struct tdigits t=q[r][0]; int ra=t.a, rb=t.b;
    for (int c=1; c<n; ++c) { t=q[m][c]; q[r][c].a=(ra+t.a)%n; q[r][c].b=(rb+t.b)%n; }
  }
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) {
    struct tdigits t=q[r][c]; x[r][c]=Permuted[t.a]*n+Permuted[t.b]+1;
  }
} // makeMargossian3

void makeGroupOdd(int **x, const int n) {
//   ------------
  int **t=tSquare; for (int r=0, i=0; r<n; ++r) for (int c=0; c<n; ++c) t[r][c]=++i;
  if ((n%3)==0) { 
    const int e=n/3, mu=(e-1)/2; int mr=(e-3)/2, i=0, j=0, k=0, *per=x[0], *rev=x[1];
    bool *u=numberUsed; for (int i=0; i<e; ++i) u[i]=F;
    // Mark mu triads to remain unchanged.
    while (T) { j=random(e); if (!u[j]) { u[j]=T; if (++i==mu) break; }}
    if (mr>0) { // Reverse mr triads.
      i=0; while (T) { j=random(e); if (!u[j]) { rev[i]=j; u[j]=T; if (++i==mr) break; }}
      int h=0; i=3*rev[h]; j=i+1; k=j+1;
      while (mr--) { swapRows(t,i,k); swapCols(t,n,i,k); i=3*rev[++h]; j=i+1; k=j+1; }
    }
    // Permute the remaining 2 triads.
    i=0; for (j=0; j<e; ++j) if (!u[j]) per[i++]=j; i=3*per[0]; j=i+1; k=j+1;
    swapRows(t,i,j); swapRows(t,j,k); swapCols(t,n,i,j); swapCols(t,n,j,k); i=3*per[1]; j=i+1; k=j+1;
    swapRows(t,i,k); swapRows(t,j,k); swapCols(t,n,i,k); swapCols(t,n,j,k);
  }
} // makeGroupOdd

void makeMargossianOddKnight(int **x, const int n) {
//   -----------------------
  makeGroupOdd(x, n); int **t=tSquare; const int m=n-1;
  do {
    const int a=random(m)+1, b=random(m)+1, a1=random(m)+1, b1=random(m)+1;
    int sum=0, R=0, C=0, R0=0, C0=0; if ((a==a1)|(b==b1)|((a+a1)==n)|((b+b1)==n)) continue;
    // Check row 0.
    for (int c=0; c<n; ++c) { x[0][c]=t[R][C]; sum+=t[R][C]; R=(R+b)%n; C=(C+a)%n; }
    if (sum!=MagicConstant) continue; R=(R0+b1)%n; C=(C0+a1)%n;
    for (int r=1; r<n; ++r) {
      R0=R, C0=C; for (int c=0; c<n; ++c) { x[r][c]=t[R][C]; R=(R+b)%n; C=(C+a)%n; } R=(R0+b1)%n; C=(C0+a1)%n;
    }
    if (isCorrect(x, n)&&isPand(x, n)) break;
  } while (T);
} // makeMargossianOddKnight

void makeFromSODLS(int **x, const int n) {
//   -------------
  int **t=tSquare; const int m=n-1;
  if (random(2)) for (int c=0; c<n; ++c) t[0][c]=c;
  else { int i=0; t[0][0]=i++; for (int c=m; c>0; --c) t[0][c]=i++; }
  for (int r=1; r<n; ++r) {
    t[r][0]=t[r-1][m-1]; t[r][1]=t[r-1][m]; for (int c=2; c<n; ++c) t[r][c]=t[r-1][c-2];
  }
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]=n*t[r][c]+t[c][r]+1;
} // makeFromSODLS

void makeOdd(int **x, const int n, const int kind) {
//   -------
  switch (kind) {
    case 1: makeMargossianOddKnight(x,n); break;  // all
    case 2: makeMargossian3(x, n);        break;  // mult 3 only
    case 3: makeSiamese(x, n);            break;  // not mult 3
    case 4: makeFromSODLS(x, n);          break;  // not mult 3
    case 5: makeWithTranspose(x, n);      break;  // not prime
    default: break;
  }
} // makeOdd

void printError(char *s) { printf("\a\nERROR: Square is NOT %s! Please report by email.\n", s); }
//   ----------

bool makeSquare(const int n, const int kind, const bool put) {
//   ----------
  int **x=xSquare, **t=tSquare; const bool se=(n&3)==2; ProgramError=F;
  if (n&1) makeOdd(x, n, kind); else if (se) makeSEven(x, t, n); else makeDEven(x, n, kind);
  if (ProgramError) return T; // One make failed, but try others..
  if (se|(Sindex>0)) panForm(x, t, n);
  if (isCorrect(x, n)) {
    if (!isPand(x, n)) { printError(se ? "near-pandiagonal" : "pandiagonal"); return F; }
  } else { printError("magic"); return F; }
  return put ? putSquare(x, n) : T;
} // makeSquare
//============================================== main =================================================

bool validOrder(const int n) {
//   ----------
  switch (n) {
  case 2: printf("\aERROR: There is no magic square of order 2.\n"); return F;
  case 3: printf("\aERROR: There is no pandiagonal magic square of order 3.\n"); return F;
  default:
    if (n<=0) { printf("\aERROR: N must be a positive integer.\n"); return F; }
  }
  return T;
}// validOrder

int main() {
//  ----
  bool ok=F;
  if (openDir()) {
    int n=0; bool another=T, inputSize=T, writeError=F; seed_rand();
    do {     
      if (inputSize) { printf("\nInput the order of the square: "); n=getInt(); }
      if (validOrder(n)) {
        MagicConstant=n&1 ? n*((n*n+1)/2) : (n/2)*(n*n+1); Sindex=0; sLen=n*n;
	      if (allocateStore(n)) {
	        FILE *wfp=NULL;
	        if ( (wfp=open_File(n))!=NULL) {
            const int kind=inputKind(n);
            int howMany=1; if (n>3) { printf("\nHow many? "); howMany=getInt(); }
            int loop=max(howMany,1000000); 
            while (loop--) { ok=makeSquare(n, kind, T); if (!ok) break; if (Sindex==howMany) break; }
            if (ok) writeError=!outputSquares(n, wfp); fclose(wfp);
            if (writeError) { perror("\a\nError writing file"); break; }
            else if (Sindex!=howMany) printf("Number of squares: %d\n", Sindex);
	        }
        }
      } // if (validOrder(n))
      if (writeError) { perror("\a\nError writing file"); ok=F; another=F; }
      else if (ok) {
        printf("\nMake another square? input y (yes) or n (no) or the order of the square: ");
        if (getYorOrder(&n)) inputSize=(n==0); else another=F;
      } else another=F;
    } while (another);
    freeStore();
  } // if (openDir())
  printf("\nPress a key to close the console.");
  while (!_kbhit()) Sleep(250); return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main