/*
 *  File:    Most-perfect.cpp
 *  Author:  S Harry White
 *  Created: 2012-01-09
 *  Updated: 2020-05-10
 *    A complete re-write to make reversible squares and convert to most-perfect. Added howMany
 *  Updated: 2022-01-11
 *    Tidy code.
 *  Updated: 2023-01-09
 *    Check if list is NULL before calling free. Fix freeSquare;
 *  Updated: 2023-03-28
 *    Fix store allocation. 
 *    Remove squareSize=0, list1Size=0, list2Size=0, listSize=0 from initGlobals().
 */

/*
 *  Makes most-perfect magic squares of order 4k.
 */

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

int N,               // The order of the square.
    Nm1,             // N-1
    Nd2,             // N/2
    Nd2m1,           // N/2-1 
    N3d4,            // 3*N/4
    N3m2d2,          // (3*N-2)/2,
    NN,              // N*N
    NNp1,            // NN+1
    magicSum,        // Nd2*NNp1
    Sum2,            // NNp1
    Sum4,            // Sum2+Sum2
    **tSquare=NULL, **xSquare=NULL,
    squareSize, list1Size, list2Size, listSize, listIndex, fieldWidth;

typedef unsigned int Uint; typedef unsigned short int Ushort;
typedef struct { Ushort x, y; } t_format; t_format *list1=NULL, *list2=NULL, *list=NULL;
const bool F=false, T=true;

int getFieldWidth(int n) { if (n<=1) return 1; int width=1; while ((n=n/10)!=0) ++width; return width; }
//  -------------

void initGlobals() {
//   -----------
 listIndex=0; Nm1=N-1; Nd2=N/2; Nd2m1=Nd2-1;
 N3d4=3*N/4; N3m2d2=(3*N-2)/2; NN=N*N; NNp1=NN+1; Sum2=NNp1; Sum4=NNp1+NNp1;
 magicSum=Nd2*NNp1; fieldWidth=getFieldWidth(NN);
} // initGlobals

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

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

void freeStore() {
//   ---------
  freeSquare(&tSquare, squareSize); freeSquare(&xSquare, squareSize); squareSize=0;
  if (list1!=NULL) { free(list1); list1=NULL; } list1Size=0;
  if (list2!=NULL) { free(list2); list2=NULL; } list2Size=0; list=NULL; listSize=0;
//  }
}

void printAllocFail() { printf("\a\nERROR: Storage allocation failed.\n"); }
//  ---------------

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

bool allocateSquares() {
//   ---------------
  bool ok=T;
  if (N>squareSize) {
    freeSquare(&tSquare,squareSize); freeSquare(&xSquare,squareSize);
    ok=allocateSquare(&tSquare,N)&&allocateSquare(&xSquare,N);
    if (ok) squareSize=N; else printAllocFail();
  }
  return ok;
} // allocateSquares

bool allocateLists() {
//   -------------
  const int size=10000; bool ok=T;
  if (size>list1Size) {
    if (list1!=NULL) { free(list1); list1=NULL; }
    ok=(list1=(t_format *) malloc(size*sizeof(t_format)))!=NULL;
    if (ok) list1Size=size;
    else { list1Size=0; if (list2!=NULL) free(list2); list2=NULL; list2Size=0; }
  }
  if (ok&(size>list2Size)) {
    if (list2!=NULL) { free(list2); list2=NULL; }
    ok=(list2=(t_format *) malloc(size*sizeof(t_format)))!=NULL;
    if (ok) list2Size=size;
    else { list2Size=0; if (list1!=NULL) free(list1); list1=NULL; list1Size=0; }
  }
  listSize=list1Size; list=list1; if (!ok) printAllocFail(); return ok;
} // allocateLists

bool increaseList(t_format **list) {
//   ------------
  bool ok=F;
  if (*list==list1) {
   const int size=list1Size+list1Size/2; t_format *tmp;
   ok=(tmp=(t_format *) malloc(size*sizeof(t_format)))!=NULL;
   if (ok) {
     for (int i=0; i<list1Size; ++i) tmp[i]=list1[i];
     if (list1!=NULL) free(list1); list1=tmp; *list=list1; list1Size=size; listSize=size;
   }
  } else {
    const int size=list2Size+list2Size/2; t_format *tmp;
    ok=(tmp=(t_format *) malloc(size*sizeof(t_format)))!=NULL;
    if (ok) {
      for (int i=0; i<list2Size; ++i) tmp[i]=list2[i];
      if (list2!=NULL) free(list2); list2=tmp; *list=list2; list2Size=size; listSize=size;
    }
  }
  if (!ok) printAllocFail(); return ok;
} // increaseList

bool allocateStore() { if (allocateLists()) return allocateSquares(); return F; }
//   -------------
//=============================================== 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

bool getYorOrder(int *n) { // if user inputs size instead of 'y' or 'n', use it
//   -----------
  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 getSize() { int n=0; scanf_s("%d", &n); clearLine(getchar()); return n; }
//  -------
//================================================== output ==============================================

const int bufSize=1024;
bool openDir() {
//   -------
  char buf[bufSize], msg[bufSize]; int sub=0;
  const char *baseName="Most-perfectSquares"; 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() {
//    ---------
  char buf[bufSize]; FILE *wfp=NULL; 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

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(FILE *wfp) {
//   -----------
  const int fw0=fieldWidth, fw=fw0+1; if (fw>maxFieldWidth) return F;
  for (int i=0; i<N; i++) {
    if (!fprintFW[fw0](wfp, xSquare[i][0])) return F;
    for (int j=1; j<N; j++) if (!fprintFW[fw](wfp, xSquare[i][j])) return F;
    if (fputc('\n', wfp)==EOF) return F;
  }
  return fputc('\n', wfp)!=EOF;
} // printSquare
//============================================= check =================================================

bool isMagic(int **x) {
//   -------
  Uint sumXY=0, sumYX=0;
  for (int i=0; i<N; ++i) {
    Uint sumR=0, sumC=0; for (int j=0; j<N; ++j) { sumR+=x[i][j]; sumC+=x[j][i]; }
    if ((sumR!=magicSum)|(sumC!=magicSum)) return F; sumXY+=x[i][N-i-1]; sumYX+=x[i][i];
  }
  return (sumXY==magicSum)&(sumYX==magicSum);
} // isMagic

bool isComplete(int **x) {
//   ----------
  for (int r=0; r<Nd2m1; ++r) { // no need to check last pair in each diag
    const int i=r+Nd2;
    for (int c=0; c<N; ++c) { int j=c+Nd2; j=j<N ? j : j-N; if ((x[r][c]+x[i][j])!=Sum2) return F; }
  }
  return T;
} // isComplete

bool isPandiagonal(int **x) {
//   -------------
  bool panBack=T, panForward=T;
  // Main diagonals already checked in isMagic.
  for (int i=1; i<N; ++i) {
    int sumYX=0, sumXY=0;
    for (int r=0, c=i; r<N; ++r, ++c) {
      const int cModn=c<N ? c : c-N; sumYX+=x[r][cModn]; sumXY+=x[r][Nm1-cModn];
    }
    if ((sumYX!=magicSum)|(sumXY!=magicSum)) return F;
  }
  return T;
} // isPandiagonal

bool isCompact(int **x) {
//   ---------
  for (int r=0; r<Nm1; ++r) {
    for (int c=0; c<Nm1; ++c) { if ((x[r][c]+x[r][c+1]+x[r+1][c]+x[r+1][c+1])!=Sum4) return F; }
    if ((x[r][Nm1]+x[r][0]+x[r+1][Nm1]+x[r+1][0])!=Sum4) return F;
    if ((x[Nm1][r]+x[0][r]+x[Nm1][r+1]+x[0][r+1])!=Sum4) return F;
  }
  return T;
} // isCompact

const char *notMostPerfect="\a\nERROR: Square is NOT most-perfect! Please report by email\n";

bool checkMostPerfect(int **x) {
//   ----------------
  if (isMagic(x)&&isPandiagonal(x)&&isComplete(x)&&isCompact(x)) return T;
  else { printf(notMostPerfect); return F; }
} // checkMostPerfect
//============================================= make =================================================

bool putFormat(const int x, const int y) {
//   ---------
  if ((listIndex<listSize)||increaseList(&list)) {
    t_format *q=&list[listIndex++]; q->x=x; q->y=y; return T;
  }
  return F;
} // putFormat

bool putFirstLevelFormats(bool *moreFormats) {
//   --------------------
  bool more=F;
  //printf("\nput unnested formats\n");
  if (!putFormat(0, 0)) return F; // separator
  if (!putFormat(N, N)) return F; /* full square */ if (!putFormat(0, 0)) return F; // separator

  for (int x=Nd2; x>=2; --x) if ((N%x)==0) { // col
     for (int m=2; m<(N/x); ++m) { const int mx=m*x; if ((mx<N)&((N%mx)==0)) { more=T; break; } }
    for (int y=2; y<=N; ++y) if ((N%y)==0) { // row
      if (!putFormat(x, y)) return F; // sub-square rectangle
      if (!putFormat(0, 0)) return F; // separator
    }
  }
  *moreFormats=more; return T;
} // putFirstLevelFormats

bool putNestedFormats(const int l, bool *moreFormats) {
//   ----------------
  //printf("put nested %d formats\n", l);
  t_format *oldList=list; bool more=F;
  if (list==list1) { list=list2; listSize=list2Size; } else { list=list1; listSize=list1Size; }
  const int oldIndex=listIndex; int count=0; listIndex=0;
  for (int i=0; i<oldIndex; ++i) {
    t_format *p=&oldList[i]; if (!putFormat(p->x, p->y)) return F; ++count;
    if (p->x==0) { 
      if (count==(l+1)) {
        p=&oldList[i-1]; const int x1=p->x;
        if ((x1+x1)<N) {
          const int i1=i-l, i2=i1+l, y1=p->y; int x=N-x1, y=y1+y1;
          while (x>x1) {
            if (N%x==0) {
              for (int m=2; m<(N/x); ++m) { const int mx=m*x; if ((mx<N)&((N%mx)==0)) { more=T; break; }}
              while (y<=N) {
                if (N%y==0) {
                  for (int j=i1; j<i2; ++j) { // copy the format
                    p=&oldList[j]; if (!putFormat(p->x, p->y)) return F;
                  }
                  if (!putFormat(x, y)) return F; // nest rectangle
                  if (!putFormat(0, 0)) return F; // separator
                }
                y+=y1;
              } // while (y<=N)
            }
            x-=x1; y=y1+y1;
          } // while (x<N)
        }
      }
      count=0;
    } // if (p->x==0)
  } // for (int i
  *moreFormats=more; return T;
} // putNestedFormats

bool putFormats() {
//   ----------
  bool moreFormats=F, ok=putFirstLevelFormats(&moreFormats);
  if (ok) { int l=1; while (moreFormats) { if (!(ok=putNestedFormats(l++, &moreFormats))) break; }}
  return ok;
} // putFormats

void fillRectangle(const int ra, const int ca, const int rb, const int cb, int *pv, int i) {
//   -------------
  int v=*pv; t_format *q=&list[--i];
  if (q->x==0) { int **p=xSquare; for (int r=ra; r<rb; ++r) for (int c=ca; c<cb; ++c) p[r][c]=v++; }
  else {
    const int x=q->x, y=q->y;
    for (int r=ra; r<rb; r+=y) for (int c=ca; c<cb; c+=x) { fillRectangle(r, c, r+y, c+x, &v, i); }
  }
  *pv=v;
} // fillRectangle

void reverseRowsRightHalf(int **x) {
//    -------------------
  for (int r=0; r<N; ++r) for (int c=Nd2; c<N3d4; ++c) {
    const int t=x[r][c]; x[r][c]=x[r][N3m2d2-c]; x[r][N3m2d2-c]=t;
  }
} // reverseRowsRightHalf

void  reverseColumnsBottomHalf(int **x) {
//    ------------------------
  for (int c=0; c<N; ++c) for (int r=Nd2; r<N3d4; ++r) {
    const int t=x[r][c]; x[r][c]=x[N3m2d2-r][c]; x[N3m2d2-r][c]=t;
  }
} // reverseColumnsBottomHalf

#define xForm(iX, jX) for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) t[iX][jX]=x[i][j]
#define xFormT(iX, jX) for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) t[iX][jX]=x[i][j]; found=T
#define xCopy() for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) x[i][j]=t[i][j]

void applyTransform(int **x) { int **t=tSquare; const int k=Nd2, l=k+1; xForm((i+k*j)%N, (k*i+l*j)%N); xCopy(); }
//   --------------

// In the square below:
//
//   a, b, c, d are the square corners, e is (N/2-1, N/2), f is (N/2-1, N/2-1)
//
//            a . . b
//            . f e .
//            . . . .
//            c . . d
//
// Orients the square corners as above where a < b < c < d.
// Then, if a > e and b > f inverts the square, else, if a > e and b < f mirrors the square.
//
void orientRSquare(int **x) {
//   -------------
  bool change=T; int **t=tSquare, minC=min(min(x[0][0], x[0][Nm1]), min(x[Nm1][0], x[Nm1][Nm1]));
  if (x[0][0]==minC) {
    if (x[0][Nm1]<x[Nm1][0]) change=F; else xForm(j,i); // YX
  } else if (x[0][Nm1]==minC) {
    if (x[Nm1][Nm1]<x[0][0]) xForm(Nm1-j,i); /* -90 */ else xForm(i,Nm1-j); // Y   
  } else if (x[Nm1][Nm1]==minC) {
    if (x[Nm1][0]<x[0][Nm1]) xForm(Nm1-i,Nm1-j); /* 180 */ else xForm(Nm1-j,Nm1-i); // XY
  } else { // x[Nm1][0]==minC
    if (x[0][0]<x[Nm1][Nm1]) xForm(j,Nm1-i); /* 90 */ else xForm(Nm1-i,j); // X
  }
  if (change) xCopy();
  if (x[0][0]>x[Nd2m1][Nd2]) {
    if (x[0][Nm1]>x[Nd2m1][Nd2m1]) xForm(Nm1-i,j); /* X */ else xForm(i,Nm1-j); /* Y */ xCopy();
  }
}  // orientRSquare

void swapRowsCols1(const int w) {
//   -------------
  const int y=Nm1-w; int **X=xSquare;
  if (random(2)) {
    int *tmp=X[w]; X[w]=X[y]; X[y]=tmp;
    if (random(2)) { for (int r=0; r<N; ++r) { const int tmp=X[r][w]; X[r][w]=X[r][y]; X[r][y]=tmp; }}
  } else {
    if (random(2)) { int *tmp=X[w]; X[w]=X[y]; X[y]=tmp; }
    for (int r=0; r<N; ++r) { const int tmp=X[r][w]; X[r][w]=X[r][y]; X[r][y]=tmp; }
  }
} // swapRowsCols1

void swapRowsCols2(const int w, const int x) {
//   -------------
  if (w!=x) {
    const int y=Nm1-w, z=Nm1-x; int **X=xSquare;
    if (random(2)) {
      int *tmp=X[w]; X[w]=X[x]; X[x]=tmp; tmp=X[y]; X[y]=X[z]; X[z]=tmp;
      if (random(2)) {
        for (int r=0; r<N; ++r) {
          int tmp=X[r][w]; X[r][w]=X[r][x]; X[r][x]=tmp; tmp=X[r][y]; X[r][y]=X[r][z]; X[r][z]=tmp;
        }
      }
    } else {
      if (random(2)) { int *tmp=X[w]; X[w]=X[x]; X[x]=tmp;  tmp=X[y]; X[y]=X[z]; X[z]=tmp; }
      for (int r=0; r<N; ++r) {
        int tmp=X[r][w]; X[r][w]=X[r][x]; X[r][x]=tmp; tmp=X[r][y]; X[r][y]=X[r][z]; X[r][z]=tmp;
      }
    }
  }
} // swapRowsCols2

void xFormReversible() {
//   ---------------
  if (random(N)) {
    int a=random(Nd2), b=random(Nd2); while (a--) swapRowsCols1(random(Nd2));
    while (b--) swapRowsCols2(random(Nd2), random(Nd2));
  }
} // xFormReversible

bool getMostPerfect() {
//   --------------
  int **x=xSquare; orientRSquare(x); reverseRowsRightHalf(x); reverseColumnsBottomHalf(x);
  applyTransform(x); return checkMostPerfect(x);
} // getMostPerfect

bool makeSquares(FILE *wfp, int *howMany, bool *writeError) {
//   -----------
  for (int i=1; i<listIndex; ++i) {
    t_format *q=&list[i];
    if (q->x==0) {
      q=&list[i-1]; const int x=q->x, y=q->y; int v=1; 
      for (int r=0; r<N; r+=y) for (int c=0; c<N; c+=x) fillRectangle(r, c, r+y, c+x, &v, i-1);
      xFormReversible(); if (!getMostPerfect()) return F;
      if (!printSquare(wfp)) { *writeError=T; return F; } if (--*howMany==0) return T;
    }
  }
  return T;
} // makeSquares
//================================================== main ==============================================

//void outputLocalTime(FILE *wfp) {
////   ---------------
//  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\n\0", &local);
//    fprintf(wfp, "\n%s", dateTime);
//  }
//} // outputLocalTime
//
//void printElapsedTime(FILE *wfp, 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

bool validOrder() {
//   ----------
  if ((N<=0)|((N&3)!=0)) {
    printf("\aERROR: Order %d is not a positive integral multiple of 4.\n", N); return F;  
  }
  return T;
} // validOrder

int main() {
//  ----
  bool ok=F;
  if (openDir()) { 
    bool another=T, inputSize=T, writeError=F; squareSize=0; list1Size=0; list2Size=0; listSize=0;
    seed_rand();
    do {     
      if (inputSize) { printf("\nInput the order of the square: "); N=getSize(); }
      if (validOrder()) {
        initGlobals();
	      if (allocateStore()) {
          if (putFormats()) {
	          FILE *wfp=NULL;
	          if ( (wfp=open_File())!=NULL) {
              printf("\nHow many? "); int howMany=getSize();
              while (howMany>0) { if (!makeSquares(wfp, &howMany, &writeError)) break; } fclose(wfp);
	          }
          }
	      } else printf("\a\nERROR: Storage allocation failed.\n");
      }
      if (writeError) { perror("\a\nError writing file"); another=F; }
      else { ok=T;
        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; 
      }
    } while (another);
    freeStore(); 
  } // if (openDir()
  printf("\nPress a key to close the console.");
  while (!_kbhit()) Sleep(250); return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main