/*
 *  File:    Associative.cpp
 *  Author:  S Harry White
 *  Created: 2011-04-04
 *  Updated: 2020-04-25
 *    Add howMany.
 *           2020-05-28
 *    Add makeOddMamzeris.
 *           2020-06-09
 *    Add makeOddWalkington.
 *    Extend makeOddMamzeris from 1 to 4.
 *    Sort output in Frénicle standard form removing duplicates.
 *           2020-06-15
 *    Add makeOddGolunski.
 *           2020-06-19
 *    Add row/col swaps of input table to increase squares made by makeOddMamzeris.
 *           2020-06-21
 *    Add makeOddFrierson.
 *           2020-06-24
 *    Add makeEvenAl_Asfizari, makeEvenPlanck, makeEvenMargossian.
 *    Add input selection of method.
 *           2020-07-22
 *    Add makeOddChan and makeOddDeLosReyes.
 *           2020-11-24
 *    Adjust startSize and startSquaresSize.
 *  Updated: 2021-12-31
 *    Tidy code.
 *  Updated: 2023-01-03
 *    Change to compile with g++ for macOS.
 */

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

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

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

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

const char *storeAllocFail="Storage allocation failed"; bool bell=T;
bool reportError(const 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(&tSquare, allocatedSize); freeSquare(&xSquare, 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, **t=(int**) malloc(size*sizeof(int*));
  if ((ok=(t!=NULL))) {
    for (int i=0; i<allocatedSquaresSize; i++) t[i]=Squares[i]; free(Squares); Squares=t;
    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(&tSquare, size)))
      if ((ok=allocateSquare(&xSquare, size)))
        if ((ok=allocateSquares(size))) {
          allocatedSize=size; allocatedSquaresSize=startSquaresSize;
          if ((ok=allocateInts(&square, size*size)))
            if ((ok=allocateInts(&Permuted, size)))
              ok=(numberUsed=(bool *) malloc((size*size+1)*sizeof(bool)))!=NULL;
        }
        if (!ok) { freeStore(); reportError(storeAllocFail); }
  }
  return ok;
} // allocateStore
//============================================== 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 isAssoc(int **x, const int n) {
//   -------
  const int g=n/2, m=n-1, sum2=n*n+1, maxCount=n&1 ? 0 : (n&3)==0 ? 0 : 4; int count=0;
  for (int r=0; r<g; ++r) for (int c=0; c<n; ++c) {
    if (x[r][c]+x[m-r][m-c]!=sum2) if (++count>maxCount) return F;
  }
  if (n&1) for (int c=0; c<g; ++c) if (x[g][c]+x[g][m-c]!=sum2) return F;
  return T;
} // isAssoc
//================================================ input ================================================

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

bool getYorOrder(int *n) {
//   -----------
  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("%d", &n); clearLine(getchar()); return n; }
//  -------  

int inputKind(const int n) {
//  ---------
  int k=0;
  if (n&1) {
    if (n<5) k=7; else {
      printf("Method, ( 1 Chan, 2 De Los Reyes, 3 Frierson, 4 Golunski,\n"
             "          5 Mamzeris I, 6 Mamzeris II, 7 Siamese, 8 Walkington)? ");
      k=getInt(); if ((k<1)|(k>8)) k=7;
    }
  } else if ((n&3)==0) {
    printf("Method, ( 1 Al_Asfizari I, 2 Al_Asfizari II, 3 Margossian, 4 Planck)? ");
    k=getInt(); if ((k<1)|(k>4)) k=1;
  }
  return k;
} // inputKind
//===================================================== output ==================================================

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

FILE *open_File(const int n) {
//    ---------
  char buf[bufSize]; FILE *wfp=NULL; snprintf(buf, bufSize, "./Order%d.txt", n);
  if ((wfp=fopen(buf, "a"))==NULL) {
    char msg[bufSize+50]; snprintf(msg, bufSize+50, "\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
//================================= sort Frenicle 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

bool pushSquare(int **x, const int n, bool *storeFail) {
//   ----------
  int i=0, insertPoint; 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)) return F;
  return 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 ============================================

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

bool isPrime(const int n) { int i=2; if (n<=1) return F; while ((i*i)<=n) if ((n%i++)==0) return F; return T; }
//   -------

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 **x, const int u, const int v) { int *t=x[u]; x[u]=x[v]; x[v]=t; }
//   --------

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

void swapRowsRand(int **x, const int n) {
//   ------------
  const int v=n/2-1, w=random_(n/2-1), y=n-1-v, z=n-1-w;
  int *t=x[v]; x[v]=x[w]; x[w]=t; t=x[y]; x[y]=x[z]; x[z]=t;
} // swapRowsRand

void swapRowsCols1(int **x, const int n, const int u) {
//   -------------
  const int y=n-1-u; const bool se=((n&3)==2);
  if (se||random_(2)) {
    int *t=x[u]; x[u]=x[y]; x[y]=t;
    if (se||random_(2)) { for (int r=0; r<n; ++r) { const int t=x[r][u]; x[r][u]=x[r][y]; x[r][y]=t; }}
  } else {
    if (random_(2)) { int *t=x[u]; x[u]=x[y]; x[y]=t; }
    for (int r=0; r<n; ++r) { const int t=x[r][u]; x[r][u]=x[r][y]; x[r][y]=t;
    }
  }
} // swapRowsCols1

void swapRowsCols2(int **x, const int n, const int u, const int v) {
//   -------------
  const int y=n-1-u, z=n-1-v; int *t; const bool se=((n&3)==2);
  if (se||random_(2)) {
    t=x[u]; x[u]=x[v]; x[v]=t; t=x[y]; x[y]=x[z]; x[z]=t;
    if (se||random_(2)) {
      for (int r=0; r<n; ++r) {
        int t=x[r][u]; x[r][u]=x[r][v]; x[r][v]=t; t=x[r][y]; x[r][y]=x[r][z]; x[r][z]=t;
      }
    }
  } else {
    if (random_(2)) { t=x[u]; x[u]=x[v]; x[v]=t; t=x[y]; x[y]=x[z]; x[z]=t; }
    for (int r=0; r<n; ++r) {
      int t=x[r][u]; x[r][u]=x[r][v]; x[r][v]=t; t=x[r][y]; x[r][y]=x[r][z]; x[r][z]=t;
    }
  }
} // swapRowsCols2
//
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 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

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

void xForm(int **x, const int n) {
//   -----
  if (n>3) {
    const int m=n-1, g=n/2; int loop=random_(g+1);
    while (loop--) {
      const int row1=random_(n); if ((n&1)&(row1==g)) { ++loop; continue; } swapRowsCols1(x, n, row1);
    }
    loop=random_(g+1);
    while (loop--) {
      const int row1=random_(n); if ((n&1)&(row1==g)) { ++loop; continue; }
      int row2=random_(n); while ((row1==row2)|(row1==(m-row2))) row2=random_(n);
      if ((n&1)&(row2==g)) { ++loop; continue; } if ((row1>=g)!=(row2>=g)) { ++loop; continue; }
      swapRowsCols2(x, n, row1, row2);
    }
  }
} // xForm

typedef void (*t_Rotate)(const int, const int, const int);
void rotate_0  (const int n, const int i, const int j) { xSquare[i][j]    =tSquare[i][j]; };
void rotate_90 (const int n, const int i, const int j) { xSquare[j][n-i]  =tSquare[i][j]; };
void rotate_180(const int n, const int i, const int j) { xSquare[n-i][n-j]=tSquare[i][j]; };
void rotate_270(const int n, const int i, const int j) { xSquare[n-j][i]  =tSquare[i][j]; };
void rotate_Y  (const int n, const int i, const int j) { xSquare[i][n-j]  =tSquare[i][j]; };
void rotate_XY (const int n, const int i, const int j) { xSquare[n-j][n-i]=tSquare[i][j]; };
void rotate_X  (const int n, const int i, const int j) { xSquare[n-i][j]  =tSquare[i][j]; };
void rotate_YX (const int n, const int i, const int j) { xSquare[j][i]    =tSquare[i][j]; };
static t_Rotate rotateCell[]=
{ rotate_0,    // no rotation
  rotate_90,   
  rotate_180,
  rotate_270,
  rotate_Y,    // about vertical axis,    (mirror)
  rotate_XY,   // about forward diagonal, (mirror + rotate_90)
  rotate_X,    // about horizontal axis,  (mirror + rotate_180)
  rotate_YX    // about back diagonal,    (mirror + rotate_270)
};

void rotateSquare(const int rotation, const int n) {
//   ------------
  const int m=n-1; t_Rotate rotate=rotateCell[rotation];
  for (int i=0; i<n; i++) for (int j=0; j<n; j++) rotate(m, i, j);
} // rotateSquare
//===================================== singly even ====================================

// [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(const int m, int r, int c) {
//  -----------
  const int f=m/2; int x=-1; 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, const int n) {
//   ---------
  const int g=n/2, which=random_(2);
  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==0 ? 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];
  }
  swapRowsRand(x, n);
}  // makeSEven
//================================================ doubly even ==============================================

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;
  }
} // 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--; }}
  }
} // 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 makeEvenMargossian(int **x, const int n) {
//   -------------------
  const int g=n/2, m=n-1, j=n+g-1; 
  for (int i=0; i<g; ++i) Permuted[i]=i; for (int i=g; i<n; ++i) Permuted[i]=j-i;
  struct tdigits **q=(struct tdigits **)x;
  for (int i=0; i<n; ++i) {
    const int beta=(g*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;
  }
  xFormA_D(x, n);
} // makeEvenMargossian

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) swapCols(x, n, col, m-col);
  }
} // 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: makeEvenMargossian(x, n);     break;
    case 4: makeEvenPlanck(x, n);         break;
    default: break;
  }
} // makeDEven

void makeEven(int **x, const int n, const int kind) {
//   --------
  if ((n&3)==0) makeDEven(x, n, kind); else makeSEven(x, n);
}  // makeEven
//==================================================== odd ====================================================

void makeOddChan(int **x, const int n) { // Chan, Mainkar, Narayan, Webster
//   -----------
  const int s1=(n*n+1)/2, m=n-1, g=n/2; int **t=tSquare; bool *u=numberUsed;

  t[0][0]=0; for (int i=1; i<=g; ++i) u[i]=F;
  for (int c=1; c<=g; ++c) {
    int a=random_(g)+1; while (u[a]) a=random_(g)+1; u[a]=T;
    if (randomBool()) a=-a; t[0][c]=a; t[0][n-c]=-a;
  }
  for (int r=1; r<n; ++r) for (int c=0; c<n; ++c) { int b=c-1; if (b<0) b=m; t[r][c]=t[r-1][b]; }
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]=s1+n*t[r][c]+t[r][m-c];
} // makeOddChan

int getFactorDLR(const int n) {
//  ------------
  const int m=n-1; int f;
  if (isPrime(n)) f=random_(n-2)+2;
  else {
    do {
      int i=3; bool good=T; f=random_(n-2)+2;
      while ((i*i)<=n) {
        if ((n%i)==0) {
          int j=n/i; 
          if ( ((f%i)==0)|(((f-1)%i)==0)) { good=false; break; }
          else if (j!=i) if ( ((f%j)==0)|(((f-1)%j)==0)) { good=false; break; }
        }
        if (!good) break; i+=2;
      }
      if (good) break;
    } while (T);
  }
  return f;
} // getFactorDLR

void makeOddDeLosReyes(int **x, const int n) { // De Los Reyes, Pourdarvish, Midha, Das
//   -----------------
  const int m=n-1, g=n/2, h=g+1, f=getFactorDLR(n); int i=h, j=n-i, **t=tSquare; t[0][0]=h;
  for (int c=1; c<n; ++c) { ++i; if (i==n) i-=n; t[0][c]=i; }
  for (int r=1; r<n; ++r) for (int c=0; c<n; ++c) t[r][c]=t[r-1][c==m?0:c+1];
  for (int c=0; c<n; ++c) x[0][c]=(f*t[0][c])%n;
  for (int r=1; r<n; ++r) for (int c=0; c<n; ++c) { int z=x[r-1][c]; z=(z==m)?0:z+1; x[r][c]=z; }
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]+=n*t[r][c]+1;
} // makeOddDeLosReyes

void makeOddFrierson(int **x, const int n) {
//   ---------------
  const int m=n-1, o=n+1, g=n/2, f=g-1, h=g+1, G=g*n; int **t=tSquare; bool *u=numberUsed;

  for (int c=1; c<=n; ++c) u[c]=F; x[g][g]=h; u[h]=T;
  for (int c=0; c<g; ++c) {
    int j=random_(n)+1; while (u[j]) j=random_(n)+1; u[j]=T; u[o-j]=T; x[g][c]=j; x[g][m-c]=o-j;
  }
  for (int c=0; c<n; ++c) u[c]=F; t[g][g]=G; u[g]=T;
  for (int c=0; c<g; ++c) {
    int j=random_(n); while (u[j]) j=random_(n); u[j]=T; u[m-j]=T; t[g][c]=j*n; t[g][m-c]=(m-j)*n;
  }
  for (int c=0; c<n; ++c) u[c]=F; u[g]=T;
  for (int r=0; r<g; ++r) {
    int C=g, c;
    if ((r==(f-1))&(!u[f])) c=random_(2)?f:m-f; // Use here to leave a choice for row f.
    else { c=random_(n); while (u[c]|(c==r)|(c==(m-r))) c=random_(n); } u[c]=T; u[m-c]=T;
    x[r][c]=h; x[m-r][m-c]=h; t[r][m-c]=G; t[m-r][c]=G; int d=m-c, loop=m;
    while (loop--) {
      ++c; if (c>=n) c-=n; ++d; if (d>=n) d-=n; ++C; if (C>=n) C-=n;
      x[r][c]=x[g][C]; x[m-r][m-c]=x[g][m-C]; t[r][d]=t[g][C]; t[m-r][m-d]=t[g][m-C];
    }
  }
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]+=t[r][c];
} // makeOddFrierson

void makeOddGolunski(int **x, const int n) {
//   ---------------
  const int m=n-1, g=n/2, h=g+1; int **t=tSquare;

  for (int r=0; r<n; r++) for (int c=0; c<n; c++) x[r][(c+h*r)%n]=(c<=g) ? 1+2*c : 2*(n-c);
  for (int r=0; r<n; r++) for (int c=0; c<n; c++) t[r][c]=x[(m-r+2)%n][(c+1)%n];
  for (int r=0; r<n; r++) for (int c=0; c<n; c++) x[r][c]=t[r][c]*n+1-x[r][c];
} // makeOddGolunski

void makeOddMamzeris(int **x, const int n, const bool xFormInput) {
//   ---------------
  const int m=n-1, g=n/2, aspect[4]={0,1,4,5}, rot=aspect[random_(4)]; int i=1, j=n-i, **t=tSquare;

  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) t[r][c]=i++;
  if (xFormInput) xForm(t, n); rotateSquare(rot, n); i=1;
  for (int r=0; r<n; ++r) {
    if (r==g) for (int c=0; c<n; ++c) t[r][c]=x[r][c];
    else {
      for (int c=0; c<i; ++c) t[r][c+j]=x[r][c]; for (int c=0; c<j; ++c) t[r][c]=x[r][c+i]; ++i; --j;
    }
  }
  i=1, j=n-i;
  for (int c=0; c<n; ++c) {
    if (c==g) for (int r=0; r<n; ++r) x[r][c]=t[r][c];
    else {
      for (int r=0; r<i; ++r) x[r+j][c]=t[r][c]; for (int r=0; r<j; ++r) x[r][c]=t[r+i][c]; ++i; --j;
    }
  }
} // makeOddMamzeris

void makeOddSiamese(int **x, const int n) {
//   --------------
  const int which=random_(2);
  for (int r=0; r<n; ++r) {
    const int rr=r+r;
    for (int c=0; c<n; ++c) {
      int i=rr+c+1, j=which==0 ? rr-c : r+c-(n-1)/2;
      if (i<0) i+=n; if (j<0) j+=n; x[r][c]=i%n+(j%n)*n+1;
    }
  }
} // makeOddSiamese

void makeOddWalkington(int **x, const int n) {
//   -----------------
  const int which=random_(2), incrR0=-2, incrC0=(which==0)?n-2:1;
  int R0=n-2, C0=(which==0)?n-1:n-n/2, i=1, loop1=n, loop2=n;
  while (loop1--) {
    int r=R0, c=C0; while (loop2--) { x[r--][c++]=i++; if (r<0) r+=n; if (c>=n) c-=n; }
    R0+=incrR0; if (R0<0) R0+=n; C0+=incrC0; if (C0>=n) C0-=n; loop2=n;
  }
} // makeOddWalkington

void makeOdd(int **x, const int n, const int kind) {
//   -------
  switch (kind) {
  case 1: makeOddChan(x, n);        break;
  case 2: makeOddDeLosReyes(x,n);   break;
  case 3: makeOddFrierson(x, n);    break;
  case 4: makeOddGolunski(x, n);    break;
  case 5: makeOddMamzeris(x, n, F); break; // no input transform
  case 6: makeOddMamzeris(x, n, T); break; // transform input
  case 7: makeOddSiamese(x, n);     break;
  case 8: makeOddWalkington(x, n);  break;
  default: break;
  }
} // makeOdd

bool makeSquare(const int n, const int kind, FILE *wfp) {
//   ----------
  int **x=xSquare;
  if (n&1) makeOdd(x, n, kind); else makeEven(x, n, kind); if (Sindex>0) xForm(x, n);
  if (isCorrect(x, n)) {
    if (!isAssoc(x, n)) return printError((n&3)==2 ? "near-associative" : "associative");
  } else return printError("magic");
  return putSquare(x, n);
} // makeSquare
//================================================== main =================================================

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

int main() {
//  ----
  bool another=T, inputSize=T, writeError=F, ok=F; seed_rand();
  if (openDir()) {
    int n=0;
    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) {
            int kind=inputKind(n), howMany=1; if (n>3) { printf("\nHow many? "); howMany=getInt(); }
            int loop=max(howMany,1000000); 
            while (loop--) { ok=makeSquare(n, kind, wfp); if (!ok) break; if (Sindex==howMany) break; }
            if (ok) writeError=!outputSquares(n, wfp); fclose(wfp);
	        }
        }
      } // 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);
  } // if (openDir())

  freeStore(); printf("\nHit return to close the console.");
  while (T) if (getchar()=='\n') break; return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main