/*
  * File:    BlockSquares.cpp
  * Author:  S Harry White
  * Created: 2019-01-20
  * Updated: 2019-04-19
  * Updated: 2021-02-02
  *   Add bordered blocks.
  * Updated: 2021-12-06
  *   Tidy code.
  *  Updated: 2021-12-29
  *    Replace & with && where required.
  *  Updated: 2022-05-26
  *    Increase maxOrder to 1000. Fix freeRectangle size for stack.
  *  Updated: 2023-01-07
  *    Fix getFactors.
  *    Remove unformatted data from printf.
  */

/*
 *  Makes block magic squares of order N.
 *  These are based on some of the squares described by Inder J.Taneja.
 *  ijtaneja@gmail.com">ijtaneja@gmail.com
 *  Formerly, Professor of Mathematics, Federal University of Santa Catarina,
 *  88040-400 Florianópolis, SC, Brazil.
 *  https://zenodo.org/search?page=1&size=20&q=%22inder%20j%20taneja%22&keywords=magic&keywords=squares&sort=mostviewed"
 */

#include "stdafx.h"
#include <assert.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>

const int maxOrder=1000, maxSquares=100, startSize=32;
int R, C,
    N,             // sqrt(R*C)
    Rm1, Cm1,      // 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

    **xSquare, **bSquare, *xRectangle, *numbers, *row, *col, *opp,
    *rowSum, *colSum, *tabu, tabuLen, squaresMade, allocated, **stack, diagLsum, diagRsum;
const bool F=false, T=true; bool odd, *numberUsed, writeError, *diagL, *diagR;

enum magicType { normalSemiMagic, normalMagic, Associative };
const magicType firstType=normalSemiMagic, lastType=Associative;
magicType rectangleType;

enum sqKind { none, semi, magic, assoc, pan, ultra };
char *kindName[]={ "Not ", "Semi", "", "Associative ", "Pandiagonal ", "Ultra" };

typedef int(*t_iP)(int *); // restart function
typedef int(*t_iPcIiP)(int *, const int, int *); // tabuSums function
//======================================= common routines ========================================

const int bufSize=1024; bool bell=T; char fmt[bufSize];
bool reportError(char *msg) { printf("%sError: %s.\n", bell?"\a\n":"", msg); return bell=F; }
//   -----------

void seed_rand() { srand((unsigned int)time(NULL)); }
//   ---------

int random(int x) { return rand()%x; }
//  ------

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

bool getFactors(const int n, int *nA, int *nB) {
//   ----------
  const int r=n/3; int num=0, *fac=tabu; bool ok=F;
  if (n&1) {
    for (int a=3; a<=r; a+=2) if ((n%a)==0) { fac[num++]=a; fac[num++]=n/a; }
  } else {
    for (int a=4; a<=r; ++a) if ((n%a)==0) { fac[num++]=a; fac[num++]=n/a; }
  }
  if (num>0) { const int a=fac[random(num)], b=n/a; ok=T;
    if (random(2)) { *nA=a; *nB=b; } else { *nA=b; *nB=a; }
  }
  //printf("factors %d %d\n", *nA, *nB);
  return ok;
} // getFactors
//======================================= allocate store =========================================

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

bool allocateRectangle(int*** rectangle, const int nR, const int nC) {
//   -----------------
  bool ok; *rectangle=(int**) malloc(nR*sizeof(int*)); ok=(*rectangle!=NULL);
  if (ok) {
    int numAllocated=0;
    for(int i=0; i<nR; i++) {
      int* p=(int*) malloc(nC*sizeof(int)); (*rectangle)[i]=p;
      if (p==NULL) { numAllocated=i; ok=F; break; }
    }
    if (!ok) freeRectangle(rectangle, numAllocated);
  }
  return ok;
} // allocateRectangle

void freeStore(const int n) {
//   ---------
  if (allocated==0) {
    stack=NULL; xSquare=NULL; bSquare=NULL; xRectangle=NULL; numbers=NULL; row=NULL; col=NULL;
    opp=NULL; rowSum=NULL; colSum=NULL; tabu=NULL; numberUsed=NULL;
  } else {
    freeRectangle(&stack,maxSquares); freeRectangle(&xSquare,n); freeRectangle(&bSquare,n);
    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 (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; }
    if (numberUsed!=NULL) { free(numberUsed); numberUsed=NULL; }
    allocated=0;
  }
} // 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(int size) {
//   -------------
  if (size<startSize) size=startSize; bool ok=T;
  if (size>allocated) {
    freeStore(allocated);
    const int sqSize=size*size;
    if (ok=allocateRectangle(&stack, maxSquares, sqSize))
      if (ok=allocateRectangle(&xSquare, size, size))
        if (ok=allocateRectangle(&bSquare, size, size))
          if (ok=allocateInts(&xRectangle, sqSize))
            if (ok=allocateInts(&numbers, sqSize))
              if (ok=allocateInts(&row, sqSize))
                if (ok=allocateInts(&col, sqSize))
                  if (ok=allocateInts(&rowSum, size))
                    if (ok=allocateInts(&colSum, size))
                      if (ok=allocateInts(&tabu, sqSize))
                        if (ok=allocateInts(&opp, sqSize))
                          if (ok=allocateBools(&diagL, sqSize))
                            if (ok=allocateBools(&diagR, sqSize))
                              ok=allocateBools(&numberUsed, sqSize+1);
                      
    if (ok) allocated=size; else { freeStore(size); reportError("Storage allocation failed"); }
  }
  return ok;
} // allocateStore
//============================================== store ===========================================

int cmpSquare(const int nn, int *x, int *xx) {
//  ---------
  int i=0; while (i++<nn) if (*x<*xx) return -1; else if (*x++>*xx++) return 1; return 0;
} // cmpSquare

int findSquare(int *x, const int nn, int *insertPoint, int **stack) {
//  ----------
  int first=0, last=squaresMade-1, middle;
  while (first<=last) {
    middle=(first+last)/2; const int cmp=cmpSquare(nn, x, stack[middle]);
    if (cmp<0) last=middle-1; else if (cmp>0) first=middle+1; else return T;
  }
  *insertPoint=first; return F;
} // findSquare

void insertSquare(int *x, const int nn, const int insertPoint, int **stack) {
//   ------------
  int *p=stack[squaresMade];
  for (int i=squaresMade; i>insertPoint; --i) stack[i]=stack[i-1];
  stack[insertPoint]=p; for (int i=0; i<nn; ++i) { p[i]=x[i]; } ++squaresMade;
} // insertSquare

bool pushSquare(int **x, const int n, int **stack) {
//   ----------
  const int nn=n*n; int insertPoint=0; int *xx=tabu, i=0;
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) xx[i++]=x[r][c];
  if (squaresMade>0) if (findSquare(xx, nn, &insertPoint, stack)) return F;
  insertSquare(xx, nn, insertPoint, stack); return T;
} // pushSquare
//============================================ input ==================================================

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

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 ( ('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, c; scanf_s("%d", &n); c=getchar(); clearLIne(c); return n; }
//  -------

//=================================================== type ====================================================

bool isAssociative(int **x, const int r0, const int c0, const int n, const int S2) {
//   -------------
  if ((n&3)==2) return F; const int nd2=n/2, mr=r0+r0+n-1, mc=c0+c0+n-1;
  for (int r=r0; r<(nd2+r0); ++r) for (int c=c0; c<(n+c0); ++c) if (x[r][c]+x[mr-r][mc-c]!=S2) return F;
  if (n&1) for (int c=c0; c<(nd2+c0); ++c) if (x[r0+nd2][c]+x[r0+nd2][mc-c]!=S2) return F;
  return T;
} // isAssociative

bool isPandiagonal(int **x, const int r0, const int c0, const int n, const int magicSum) {
//   -------------
  if ((n&3)==2) return F; const int mc=c0+c0+n-1;
  for (int i=r0; i<(n+r0-1); ++i) {
    int sumYX=0, sumXY=0;
    for (int r=r0, c=i-r0+c0+1; r<(n+r0); ++r, ++c) { // i -1, c +1 because main diags checked in isMagic
      const int cModn=(c<(n+c0))?c:c-n; sumYX+=x[r][cModn]; sumXY+=x[r][mc-cModn];
    }
    if ((sumYX!=magicSum)|(sumXY!=magicSum)) return F;
  }
  return T;
} // isPandiagonal

sqKind isMagic(int **M, const int r0, const int c0, const int n, const int nn, const int magicSum, const bool u) {
//     -------
  const int m=n+r0+c0+c0-1; int sumX, sumY, sumXY=0, sumYX=0;

  if (u) for (int i=0; i<=nn; i++) numberUsed[i]=F;
  for (int i=r0; i<(n+r0); i++) {
    sumX=0; sumY=0;
    for (int j=c0; j<(n+c0); j++) {
      const int t=M[i][j]; if(u) numberUsed[t]=T; sumX+=t; sumY+=M[j+r0-c0][i-r0+c0];
    }
    if ((sumX!=magicSum)|(sumY!=magicSum)) return none;
    sumXY+=M[i][m-i-c0]; sumYX+=M[i][i+c0-r0];
  }
  if (u) for (int i=1; i<=nn; i++) if (!numberUsed[i]) return none;
  return ((sumXY==magicSum)&(sumYX==magicSum))?magic:semi;
} // isMagic

sqKind squareKind(int **M, const int n) {
//     ----------
  const int nn=n*n, S2=nn+1, magicSum=n&1?n*(S2/2):(n/2)*S2;
  sqKind kind=isMagic(M,0,0,n,nn,magicSum,T);
  if (kind==magic) {
    if (n==1) return magic;
    if (isAssociative(M,0,0,n,S2)) kind=isPandiagonal(M,0,0,n,magicSum)?ultra:assoc;
    else if (isPandiagonal(M,0,0,n,magicSum)) kind=pan;
  }
  return kind;
} // squareKind

sqKind subSquareKind(int **M, const int n, const int r0, const int c0, int *sum) {
//     -------------
  const int nn=n*n, m=n-1, S2=M[r0][c0]+M[r0+m][c0+m]; int magicSum=0;
  for (int c=c0; c<(c0+n); ++c) magicSum+=M[r0][c];
  sqKind kind=isMagic(M,r0,c0,n,nn,magicSum,F);
  if (kind==magic) {
    if (n==1) return magic;
    if (isAssociative(M,r0,c0,n,S2)) kind=isPandiagonal(M,r0,c0,n,magicSum)?ultra:assoc;
    else if (isPandiagonal(M,r0,c0,n,magicSum)) kind=pan;
  }
  *sum=magicSum; return kind;
} // subSquareKind

sqKind resolveKind(sqKind k0, const sqKind k, const int s0, const int s, bool *eq) {
//     -----------
  if (k0!=k) {
    if ((k0==semi)&(k!=magic)&(k!=assoc)&(k!=pan)&(k!=ultra)) k0=none;
    else if ((k0==magic)&(k!=assoc)&(k!=pan)&(k!=ultra)) {
      if (k==semi) k0=semi; else k0=none;
    } else if (k0==assoc) {
      if (k==semi) k0=semi; else if (k==magic) k0=magic; else if (k!=ultra) k0=none;
    } if (k0==pan) {
      if (k==semi) k0=semi; else if (k==magic) k0=magic; else if (k!=ultra) k0=none;
    } if (k0==ultra) {
      if (k==semi) k0=semi; else if (k==magic) k0=magic; else if (k==assoc) k0=assoc;
      else if (k==pan) k0=pan; else k0=none;
    }
  }
  *eq=(s==s0); return k0;
} // resolveKind

sqKind blockKind(int **x, const int n, const int b, const int o, int *bmsum) {  // equal sums
//     ---------
  const int no=n-o; int sum=0, c0=o+b; bool equal=T, rbreak=F;
  sqKind kind=subSquareKind(x,b,o,o,&sum);
  for (int r=o; r<no; r+=b) {
    for (int c=c0; c<no; c+=b) {
      int sum_=0; bool eq=F; sqKind kind_=subSquareKind(x,b,r,c,&sum_);
      kind=resolveKind(kind,kind_,sum,sum_,&eq);
      if (!eq) equal=F; if (kind==none) { rbreak=T; break; }
    }
    c0=o; if (rbreak) break;
  }
  if (!equal)  { printf("Program error: unequal block sums.\n"); sum=0; } *bmsum=sum; return kind;
} // blockKind

sqKind blockKindBB(int **x, const int n, int b, const int bb, int *bmsum) {  // equal sums
//     -----------
  int sum=0; sqKind kind=none; bool equal=T, bbreak=F;
  if (b==0) {
    b=(n-bb-bb)/bb; int c0=b+3; kind=subSquareKind(x,b,1,1,&sum);
    for (int r=1; r<n; r+=b+2) {
      for (int c=c0; c<n; c+=b+2) {
        int sum_=0; bool eq=F; sqKind kind_=subSquareKind(x,b,r,c,&sum_);
        kind=resolveKind(kind,kind_,sum,sum_,&eq);
        if (!eq) equal=F; if (kind==none) { bbreak=T; break; }
      }
      c0=1; if (bbreak) break;
    } 
  } else {
    int bbn=(n-bb-bb)/bb, c0=b; kind=subSquareKind(x,b,1,1,&sum);
    for (int rb=1; rb<n; rb+=bbn+2) {
      for (int cb=1; cb<n; cb+=bbn+2) {
        for (int r=0; r<bbn; r+=b) {
          for (int c=c0; c<bbn; c+=b) {
            int sum_=0; bool eq=F; sqKind kind_=subSquareKind(x,b,r+rb,c+cb,&sum_);
            kind=resolveKind(kind,kind_,sum,sum_,&eq);
            if (!eq) equal=F; if (kind==none) { bbreak=T; break; }
          }
          c0=0; if (bbreak) break;
        }
        if (bbreak) break;
      }
      if (bbreak) break;
    }
  }
  if (!equal) { printf("Program error: unequal block sums.\n"); sum=0; } *bmsum=sum; return kind;
} // blockKindBB

sqKind blockKindU(int **x, const int n, const int b) { // unequal sums
//     ----------
  int sum=0, c0=b; bool equal=T, rbreak=F;
  sqKind kind=subSquareKind(x,b,0,0,&sum);
  for (int r=0; r<n; r+=b) {
    for (int c=c0; c<n; c+=b) {
      int sum_=0; bool eq=F; sqKind kind_=subSquareKind(x,b,r,c,&sum_);
      kind=resolveKind(kind,kind_,sum,sum_,&eq);
      if (!eq) equal=F; if (kind==none) { rbreak=T; break; }
    }
    c0=0; if (rbreak) break;
  }
  if (equal) printf("Program error: equal block sums.\n"); return kind;
} // blockKindU
//============================================ output =================================================

bool openDir() {
//   -------
  int sub=0; char buf[bufSize], msg[bufSize];
  const char *baseName="BlockSquares";
  strcpy_s(buf, bufSize, baseName);
  
  do {
    if (_mkdir(buf)!=-1) 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)==-1) {
    sprintf_s(msg, bufSize, "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 *ext) {
//    ---------
 FILE *wfp=NULL; int sub=0; char buf[bufSize], defaultName[bufSize];
  sprintf_s(defaultName, bufSize, "%d", n, ext);
  sprintf_s(buf, bufSize, "%s.%s", defaultName, ext);
  do {
    if (_access_s(buf, 00)==ENOENT) break;
    else sprintf_s(buf, bufSize, "%s_%d.%s", defaultName, ++sub, ext);
  } while (T);
  if (fopen_s(&wfp, buf, "w")==0) {
    printf("%s %s%s", ext[0]=='t'?"\nOutput files are":" and", buf, ext[0]=='t'?"":".\n");
  } else {
    char msg[bufSize+50];
    sprintf_s(msg, bufSize+50, "\a\nCan't open for write %s", buf); perror(msg);
  }
  return wfp;
} // open_File

int fWidthQ(const int n) { if (n==1) return 1; int i=n, width=1; while ((i=i/10)!=0) ++width; return width; }
//  -------
int fWidth(const int n) { if (n==1) return 1; int i=n*n, width=1; while ((i=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; }

static t_fprintFW fprintFW[]={ NULL, fprintFW1, fprintFW2, fprintFW3, fprintFW4,
                                     fprintFW5, fprintFW6, fprintFW7, fprintFW8 };
const int maxFieldWidth=8;

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

bool txtPrint(int **x, const int n, FILE *wfp) {
//   --------
  //if (txtPrintQ(x,n,wfp)&&txtPrintR(x,n,wfp)) {
    int fw0=fWidth(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;
  //}
  return F;
} // txtPrint

bool writeHead(const int n, FILE *wfp) {
//   ---------
  if (fprintf(wfp,
    "<!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML 1.0 Strict//EN\""
    " \"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd\">\n"
    "<html lang=\"en\" xml:lang=\"en\" "
    "xmlns=\"http://www.w3.org/1999/xhtml\">\n\n"
    "<head>\n"
    "<title>Order %d</title>\n"
    "<style type=\"text/css\" id=\"internalStyle\">\n"
    "  body { font-family: Verdana, Arial, Helvetica, sans-serif; color : #000080 }\n"
    "  h3 { text-align: center }\n"
    "  span.m { font-size: medium } span.s { font-size: small } \n"
    "  table { margin-left: auto; margin-right: auto; border: 4px solid #aaaaaa }\n"
    "  td { %sfont-size: small; font-weight: bold; text-align: right;\n"
    "       color: #000080; border-right: none; border-bottom: none }\n"
    "    .no, .lt, .to, .tl { background-color : #eeeeee }\n", n, n==10?"width: 26px; ":"")<0) return F;
  if (fprintf(wfp,"    .nog, .ltg, .tog, .tlg { background-color : #77ff66 }\n")<0) return F;
  if (fprintf(wfp,"    .nob, .ltb, .tob, .tlb { background-color : #66ffff }\n")<0) return F;
  if (fprintf(wfp,"    .no, .nob, .nog { border-left: none; border-top: none }\n")<0) return F;
  if (fprintf(wfp,"    .lt, .ltb, .ltg { border-left: 2px solid #999999; border-top: none }\n")<0) return F;
  if (fprintf(wfp,"    .to, .tob, .tog { border-top: 2px solid #999999; border-left: none }\n")<0) return F;
  return fprintf(wfp,
    "    .tl, .tlb, .tlg { border-top: 2px solid #999999; border-left: 2px solid #999999 }\n"
    "</style>\n</head>\n\n<body>\n")>0;
} // writeHead

bool writeCaption(int **x, const int n, const int b, const int o, const int bb, FILE *wfp) {
//  -------------
  const int noo=n-o-o-bb-bb, a=noo/b; int bmsum=0;
  sqKind sKind=squareKind(x,n);
  if (fprintf(wfp, "<h3>Order %d %sMagic Square\n<span class=\"m\">", n,
       kindName[sKind])<0) return F;
  if (o>0) {
    sqKind iKind=blockKind(x,n,noo,o,&bmsum);
    if (fprintf(wfp, "<br />Order %d %sMagic Square with Border%s",
         noo, kindName[iKind], o==2?"s":"")<0) return F;
    if ((a==1)&(b>1))
      if (fprintf(wfp, "<br /><span class=\"s\">order %d magic sum %d</span>", noo, bmsum)<0) return F;
  }
  if ((a>1)&(b>1)) {
    if (bb>0) {
      const int t=noo/bb; sqKind iKind=blockKindBB(x,n,0,bb,&bmsum);
      if (fprintf(wfp, 
        "<br />%d %sMagic %dx%d Blocks With Border<br /><span class=\"s\">",bb*bb,kindName[iKind], t, t)<0) return F;
      if (t>b) {
        sqKind bKind=iKind=blockKindBB(x,n,b,bb,&bmsum);
        if (fprintf(wfp, "each having %d %sMagic %dx%d blocks with magic sum %d</span>",
          t*t/(b*b), kindName[bKind], b, b, bmsum)<0) return F;
      } else if (fprintf(wfp, "each %dx%d magic sum %d</span>", b, b, bmsum)<0) return F;
    } else {
      sqKind bKind=(o>=0)?blockKind(x,n,b,o,&bmsum):blockKindU(x,n,b);
      if (fprintf(wfp, 
        "<br />%d %sMagic %dx%d Blocks<br /><span class=\"s\">",a*a,kindName[bKind], b, b)<0) return F;
      if (o>=0) {
        if (fprintf(wfp, "each with magic sum %d</span>", bmsum)<0) return F;
      } else {
        if (fprintf(wfp, "with unequal magic sums</span>")<0) return F;
      }
    }
  }
  if (fprintf(wfp, "</span></h3>\n")<0) return F;
  return T;
} // writeCaption

bool writeCell(FILE *wfp, char *_class, int value) {
  return fprintf(wfp, "<td class=\"%s\">%d\n", _class, value)>0; }

bool writeSquareHTML2(int **x, const int n, const int b, FILE *wfp) {
//   ----------------
  const int m=n-1, l=m-1;
  if (fprintf(wfp, "<table frame=\"box\" summary=\"magic square\">\n")<0) return F;
  if (fprintf(wfp, "<tr>\n")<0) return F;
  for (int c=0; c<n; ++c) if (!writeCell(wfp,"nog",x[0][c])) return F;
  if (fprintf(wfp, "<tr>\n")<0) return F;
  if (!writeCell(wfp,"nog",x[1][0])) return F;
  if (!writeCell(wfp,"tlb",x[1][1])) return F;
  for (int c=2; c<m; ++c) if (!writeCell(wfp,"tob",x[1][c])) return F;
  if (!writeCell(wfp,"ltg",x[1][m])) return F;

  for (int r=2; r<l; ++r) {
    if (fprintf(wfp, "<tr>\n")<0) return F;
    if (!writeCell(wfp,"nog",x[r][0])) return F;
    if (!writeCell(wfp,"ltb",x[r][1])) return F;
    for (int c=2; c<l; ++c) {
      char *cl=(c%b==2)?(r%b==2)?"tl":"lt":(r%b==2)?"to":"no";  if (b==1) cl="tl";
      if (!writeCell(wfp,cl,x[r][c])) return F;
    }
    if (!writeCell(wfp,"ltb",x[r][l])) return F;
    if (!writeCell(wfp,"ltg",x[r][m])) return F;
  }
  if (fprintf(wfp, "<tr>\n")<0) return F;
  if (!writeCell(wfp,"nog",x[l][0])) return F;
  if (!writeCell(wfp,"ltb",x[l][1])) return F;
  for (int c=2; c<l; ++c) if (!writeCell(wfp,"tob",x[l][c])) return F;
  if (!writeCell(wfp,"nob",x[l][l])) return F;
  if (!writeCell(wfp,"ltg",x[l][m])) return F;

  if (fprintf(wfp, "<tr>\n")<0) return F;
  if (!writeCell(wfp,"nog",x[m][0])) return F;
  if (!writeCell(wfp,"tog",x[m][1])) return F;
  for (int c=2; c<m; ++c) if (!writeCell(wfp,"tog",x[m][c])) return F;
  if (!writeCell(wfp,"nog",x[m][m])) return F;
  return fprintf(wfp, "</table>\n")>0;
} // writeSquareHTML2

bool writeSquareHTML1(int **x, const int n, const int b, FILE *wfp) {
//   ----------------
  const int m=n-1;
  if (fprintf(wfp, "<table frame=\"box\" summary=\"magic square\">\n")<0) return F;
  if (fprintf(wfp, "<tr>\n")<0) return F;
  for (int c=0; c<n; ++c) if (!writeCell(wfp,"nog",x[0][c])) return F;

  for (int r=1; r<m; ++r) {
    if (fprintf(wfp, "<tr>\n")<0) return F;
    if (!writeCell(wfp,"nog",x[r][0])) return F;
    for (int c=1; c<m; ++c) {
      char *cl=(c%b==1)?(r%b==1)?"tl":"lt":(r%b==1)?"to":"no"; if (b==1) cl="tl";
      if (!writeCell(wfp,cl,x[r][c])) return F;
    }
    if (!writeCell(wfp,"ltg",x[r][m])) return F;
  }
  if (fprintf(wfp, "<tr>\n")<0) return F;
  if (!writeCell(wfp,"nog",x[m][0])) return F;
  for (int c=1; c<m; ++c) if (!writeCell(wfp,"tog",x[m][c])) return F;
  if (!writeCell(wfp,"nog",x[m][m])) return F;
  return fprintf(wfp, "</table>\n")>0;
} // writeSquareHTML1

bool writeSquareHTMLbb(int **x, const int n, const int b, const int bb, FILE *wfp) {
//   -----------------
  const int m=n-1, bbn=(n-bb-bb)/bb;
  if (fprintf(wfp, "<table frame=\"box\" summary=\"magic square\">\n")<0) return F;

  int r0=1, rm=bbn+1;
  for (int i=0; i<bb; ++i) { 
    if (fprintf(wfp, "<tr>\n")<0) return F;
    int C0=1, Cm=bbn+1;
    for (int j=0; j<bb; ++j) {
      if (!writeCell(wfp,C0==1?r0==1?"nog":"tog":r0==1?"ltg":"tlg",x[r0-1][C0-1])) return F;
      for (int c=C0; c<=Cm; ++c) if (!writeCell(wfp,r0==1?"nog":"tog",x[r0-1][c])) return F;
      C0+=bbn+2; Cm+=bbn+2;
    }

    for (int r=r0; r<rm; ++r) {
      if (fprintf(wfp, "<tr>\n")<0) return F;
      int c0=1, cm=bbn+1;
      for (int j=0; j<bb; ++j) {
        if (!writeCell(wfp,c0==1?"nog":"ltg",x[r][c0-1])) return F;
        for (int c=c0; c<cm; ++c) {
          char *cl=((c-c0)%b==0)?((r-r0)%b==0)?"tl":"lt":((r-r0)%b==0)?"to":"no"; if (b==1) cl="tl";
          if (!writeCell(wfp,cl,x[r][c])) return F;
        }
        if (!writeCell(wfp,"ltg",x[r][cm])) return F;
        c0+=bbn+2; cm+=bbn+2;
      }
    }
    r0+=bbn+2; rm+=bbn+2; C0=1; Cm=bbn+1;
    if (fprintf(wfp, "<tr>\n")<0) return F;
    for (int j=0; j<bb; ++j) {
      if (!writeCell(wfp,C0==1?"nog":"ltg",x[r0-2][C0-1])) return F;
      for (int c=C0; c<Cm; ++c) if (!writeCell(wfp,"tog",x[r0-2][c])) return F;
      if (!writeCell(wfp,"nog",x[r0-2][Cm])) return F;
      C0+=bbn+2; Cm+=bbn+2;
    }
  }
  return fprintf(wfp, "</table>\n")>0;
 } // writeSquareHTMLbb

bool writeSquareHTML(int **x, const int n, const int b, FILE *wfp) {
//   ---------------
  if (fprintf(wfp, "<table frame=\"box\" summary=\"magic square\">\n")<0) return F;
  for (int r=0; r<n; ++r) {
    if (fprintf(wfp, "<tr>\n")<0) return F;
    for (int c=0; c<n; ++c) {
      char *cl=(c>0)&(c%b==0)?(r>0)&(r%b==0)?"tl":"lt":(r>0)&(r%b==0)?"to":"no";
      if (!writeCell(wfp,cl,x[r][c])) return F;
    }
  }
  return fprintf(wfp, "</table>\n")>0;
} // writeSquareHTML

bool htmlPrint(int **x, const int n, const int b, int borders, int bb, FILE *wfp) {
//   ---------
  if (bb<0) bb=0;
  bool ok=writeCaption(x,n,b,borders,bb,wfp); if (borders<0) borders=0;
  if (ok) ok=(borders==0)?(bb==0)?writeSquareHTML(x,n,b,wfp):writeSquareHTMLbb(x,n,b,bb,wfp):
             (borders==1)?writeSquareHTML1(x,n,b,wfp):writeSquareHTML2(x,n,b,wfp);
  return ok;
} // htmlPrint

bool writeFoot(FILE *wfp) { return fprintf(wfp, "</body>\n</html>")>0; }
//=============================================== make =================================================
//----------------------------------------- start local search -----------------------------------------

bool initialize(int *x, const int r, const int c, const magicType t) {
//   ----------
  R=r; C=c; Rm1=R-1; Rd2=R/2; Cm1=C-1; Cd2=C/2; RC=R*C; RCd2=RC/2; RCm1=RC-1; RCm2=RCm1-1; 
  S2=RC+1; S1=S2/2; SR=C*S2/2; SC=R*S2/2; odd=R&1; N=(R==C)?R:(int)sqrt((double)RC);
  tabuLen=(N==4)?5:(N+N)/3; rectangleType=t;
  for (int i=0; i<RC; ++i) {
    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 (t==Associative) { for (int i=0; i<RC; ++i) opp[i]=RCm1-i; if (odd) x[RCd2]=S1; }
  return T;
} // initialize

int lastViolation;
int updateViolationSemi(int i1, int i2) {
//  -------------------
  int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2];
  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);
  }

  return lastViolation+=(newv-oldv);
} // updateViolationSemi

int updateViolationMagic(int i1, int i2) { // R==C
//  --------------------
  int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2];
  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]-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];
    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 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

void randomFillMagic(int *X) {
//   ---------------
  for (int i=0; i<RCm1; ++i) {
    int ri=i+random(RC-i);
    if (i!=ri) { 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) {
    int ri=i+random(RCd2-i), v=numbers[i]; numbers[i]=numbers[ri]; numbers[ri]=v;
  }
  int next=0;
  if (R<=C) {
    for (int r=0; r<Rm1; ++r) {
      int cLim=(r<Rd2)?C-r:(odd&(r==Rd2))?Cd2:Rm1-r;
      for (int c=0; c<cLim; ++c) {
        int i=r*C+c, v=numbers[next++]; if (random(2)) v=S2-v; X[i]=v; X[opp[i]]=S2-v;
      }
    }
  } else {
    for (int c=0; c<Cm1; ++c) {
      int rLim=(c<Cd2)?R-c:(odd&(c==Cd2))?Rd2:Cm1-c;
      for (int r=0; r<rLim; ++r) {
        int i=c*R+r, 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(); }
//  ----------

void doSwapMagic(int *X, const int i1, const int i2) { 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) {
//  -------------
  int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2];
  int x1=X[i1], x2=X[i2], x, y=x2-x1, 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
//  --------------
  int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2];
  int x1=X[i1], x2=X[i2], x, y=x2-x1, 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) {
//  ------------
  int r1=row[i1], c1=col[i1], r2=row[i2], c2=col[i2];
  int x1=X[i1], x2=X[i2], x, y=x2-x1, 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 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) {
        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) {
        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) {
            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; }
      }
      int mxnSum=max(minRSum, minCSum), mnxSum=min(maxRSum, maxCSum);
      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) {
            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) {
        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) {
        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) {
            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
//---------------------------------------- end local search -------------------------------------------

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

sqKind getRowShift2(int **Q, const int n) {
//     ------------
  const int k=random(3); sqKind favor=(k==0)?assoc:(k==1)?pan:ultra;
  if (favor==pan) { for (int c=0; c<n; ++c) Q[0][c]=c; }
  else { Q[0][0]=0; for (int c=1; c<n; ++c) Q[0][c]=n-c; } // assoc, ultra
  for (int r=1; r<n; ++r) for (int c=0; c<n; ++c) {
    const int k=c>1?c-2:c+n-2; Q[r][c]=Q[r-1][k];
  }
  if (favor==assoc) { swapRows(Q,1,n-2); swapCols(Q,n,1,n-2); } return favor;
} // getRowShift2

void reverseRowsRightHalf(int **Q, const int n) {
//    -------------------
  const int n3=3*n, n3d4=n3/4, n3m2d2=(n3-2)/2, nd2=n/2;
  for (int r=0; r<n; ++r) for (int c=nd2; c<n3d4; ++c) {
    const int t=Q[r][c]; Q[r][c]=Q[r][n3m2d2-c]; Q[r][n3m2d2-c]=t;
  }
} // reverseRowsRightHalf

void reverseColumnsBottomHalf(int **Q, const int n) {
//   ------------------------
  const int n3=3*n, n3d4=n3/4, n3m2d2=(n3-2)/2, nd2=n/2;
  for (int r=nd2; r<n3d4; ++r) { int *t=Q[r]; Q[r]=Q[n3m2d2-r]; Q[n3m2d2-r]=t; }
} // reverseColumnsBottomHalf

void XformA_D(int **Q, const int n) { reverseColumnsBottomHalf(Q,n); reverseRowsRightHalf(Q,n); }
//   --------
//---------------------------------------------- start border ----------------------------------------

typedef void (*t_Turn)(int**, int**, int, int, int);
void no_Turn   (int **x, int **b, int n, int i, int j) { x[i][j]    =b[i][j]; };
void rotate_90 (int **x, int **b, int n, int i, int j) { x[j][n-i]  =b[i][j]; };
void rotate_180(int **x, int **b, int n, int i, int j) { x[n-i][n-j]=b[i][j]; };
void rotate_270(int **x, int **b, int n, int i, int j) { x[n-j][i]  =b[i][j]; };
void rotate_Y  (int **x, int **b, int n, int i, int j) { x[i][n-j]  =b[i][j]; };
void rotate_XY (int **x, int **b, int n, int i, int j) { x[n-j][n-i]=b[i][j]; };
void rotate_X  (int **x, int **b, int n, int i, int j) { x[n-i][j]  =b[i][j]; };
void rotate_YX (int **x, int **b, int n, int i, int j) { x[j][i]    =b[i][j]; };
static t_Turn turnCell[]={ no_Turn,  rotate_90, rotate_180, rotate_270,
                           rotate_Y, rotate_XY, rotate_X,   rotate_YX };
const int numRotations=8;
void Turn(int **x, int **b, int o, int n) {
//   ----
  int nmo=n-o, npo=n+o, k=random(numRotations);

  if ((n-o==0)|(n-o==3)) { // turn all 1x1, 4x4
    for (int i=o; i<=n; i++) for (int j=o; j<=n; j++) turnCell[k](x, b, npo, i, j);
  } else { // turn only border rows and columns  
    for (int i=o; i<=n; i+=nmo) for (int j=o; j<=n; j++) turnCell[k](x, b, npo, i, j);
    for (int j=o; j<=n; j+=nmo) for (int i=o+1; i<n; i++) turnCell[k](x, b, npo, i, j);
  }
} // Turn

void makeActualBorder(int **x, const int n) {
//   ----------------
  const int nn=n*n, m=n-1;
  if (n&1) {
    const int midc=(nn+1)/2;
    for (int c=0; c<n; c++) { x[0][c]+=midc; x[m][c]+=midc; }
    for (int r=1; r<m; r++) { x[r][0]+=midc; x[r][m]+=midc; }
  } else {
    const int pPlus=nn/2, mPlus=pPlus+1;
    for (int c=0; c<n; c++) { x[0][c]+=x[0][c]>0?pPlus:mPlus; x[m][c]+=x[m][c]>0?pPlus:mPlus; }
    for (int r=1; r<m; r++) { x[r][0]+=x[r][0]>0?pPlus:mPlus; x[r][m]+=x[r][m]>0?pPlus:mPlus; }
  }
} // makeActualBorder

void putBorder(int **x, int **b, int o, int n) { // random choices from row and col
//   ---------
  int length=n-o-1; // remaining choices in row, col 
  for (int i=o+1; i<(n-1); ++i) {
    int k=random(length); b[o][i]=row[k]; b[n][i]=-row[k];
    for (int j=k; j<length; ++j) row[j]=row[j+1]; k=random(length);
    b[i][o]=col[k]; b[i][n]=-col[k]; for (int j=k; j<length; ++j) col[j]=col[j+1]; --length;
  }
  b[o][n-1]=row[0]; b[n][n-1]=-row[0]; b[n-1][o]=col[0]; b[n-1][n]=-col[0]; Turn(x, b, o, n);
} // putBorder

void makeOddBorder(int **x, int **b, const int size) {
//   -------------
  int v=((size-2)*(size-2)+1)/2, m=(size-1)/2, c=0, r=0, n=size-1, k=(size-3)/2, j=k;
  while (j--) { row[c++]=v++; col[r++]=-(v++); }
  b[n][0]=-v; b[0][n]=v++; col[r++]=-(v++); b[n][n]=-v; b[0][0]=v++;
  row[c++]=-(v++); j=k; while (j--) { row[c++]=-(v++); col[r++]=v++; }
  putBorder(x, b, 0, n); makeActualBorder(x, size);
} // makeOddBorder

void makeEvenBorder(int **x, int **b, const int size) { // singly even
//   --------------
  int v=(size-2)*(size-2)/2+1, i=size, n=size-1, c=0, r=0;
  if ((i&3)==2) {
    int k=(i-2)/4, j=k;
    while (j--) { row[c++]=-(v++); row[c++]=v++; } row[c++]=-(v++); col[r++]=-(v++);
    j=k-1; while (j--) { col[r++]=v++; col[r++]=-(v++); }
    b[n][n]  =-v; b[0][0]=v++; b[n][0]=-v; b[0][n]=v++; col[r++]=v++;
    j=k-1; while (j--) { col[r++]=-(v++); col[r++]=v++; } col[r++]=v++; row[c++]=-(v++);
    j=k-1; while (j--) { row[c++]=v++; row[c++]=-(v++); } col[r]=-v;
  } else {
    int k=(i-4)/4, j=k;
    row[c++]=-(v++); while (j--) { row[c++]=v++; row[c++]=-(v++); }
    if (i==8) {
      col[r++]=v++;
    } else {
      col[r++]=-(v++); col[r++]=v++;
      j=k-2; while (j--) { col[r++]=-(v++); col[r++]=v++; } col[r++]=v++;
    }
    col[r++]=-(v++); col[r++]=-(v++);
    b[n][n]=-v; b[0][0]=v++; b[n][0]=-v; b[0][n]=v++; col[r++]=v++; col[r++]=v++;
    j=k-1; while (j--) { col[r++]=-(v++); col[r++]=v++; } col[r++]=-(v++);
    j=k; while (j--) { row[c++]=-(v++); row[c++]=v++; } row[c]=-v;
  }
  putBorder(x, b, 0, n); makeActualBorder(x, size);
}  // makeEvenBorder

//----------------------------------------------------- end border ---------------------------------------------------

char u9[][81]={
{1,50,69,56,21,43,27,76,17,44,54,22,15,28,77,70,2,48,75,16,29,49,71,0,23,42,55,13,35,72,68,6,46,39,61,20,47,66,7,18,40,
62,73,14,33,60,19,41,34,74,12,8,45,67,25,38,57,80,9,31,51,64,5,32,78,10,3,52,65,58,26,36,63,4,53,37,59,24,11,30,79},
{3,70,47,24,55,41,9,76,35,44,18,58,29,12,79,50,6,64,73,32,15,67,53,0,61,38,21,31,17,72,52,2,66,37,23,60,69,46,5,54,40,
26,75,34,11,20,57,43,14,78,28,8,63,49,59,42,19,80,27,13,65,48,7,16,74,30,1,68,51,22,62,36,45,4,71,39,25,56,33,10,77},
{5,36,79,31,65,24,57,10,53,64,26,30,9,52,59,38,78,4,51,58,11,80,3,37,25,32,63,20,33,67,46,62,12,72,7,41,61,14,45,6,40,
74,35,66,19,39,73,8,68,18,34,13,47,60,17,48,55,43,77,0,69,22,29,76,2,42,21,28,71,50,54,16,27,70,23,56,15,49,1,44,75},
{7,36,77,13,51,56,19,30,71,48,62,10,27,68,25,42,74,4,65,22,33,80,1,39,59,16,45, 60,11,49,66,26,28,72,5,43,23,34,63,2,
40,78,17,46,57,37,75,8,52,54,14,31,69,20,35,64,21,41,79,0,47,58,15,76,6,38,55,12,53,70,18,32,9,50,61,24,29,67,3,44,73},
{11,30,79,8,45,67,23,42,55,58,26,36,73,14,33,70,2,48,51,64,5,39,61,20,27,76,17,37,59,24,34,74,12,49,71,0,3,52,65,18,40,
62,15,28,77,80,9,31,68,6,46,56,21,43,63,4,53,60,19,41,75,16,29,32,78,10,47,66,7,44,54,22,25,38,57,13,35,72,1,50,69},
{57,10,53,72,7,41,69,22,29,38,78,4,35,66,19,50,54,16,25,32,63,13,47,60,1,44,75,31,65,24,46,62,12,43,77,0,9,52,59,6,40,
74,21,28,71,80,3,37,68,18,34,56,15,49,5,36,79,20,33,67,17,48,55,64,26,30,61,14,45,76,2,42,51,58,11,39,73,8,27,70,23},
{61,50,9,8,75,37,33,22,65,38,6,76,63,34,23,10,62,48,21,64,35,49,11,60,77,36,7,67,29,24,14,54,52,39,1,80,53,12,55,78,40,
2,25,68,27,0,79,41,28,26,66,56,51,13,73,44,3,20,69,31,45,16,59,32,18,70,57,46,17,4,74,42,15,58,47,43,5,72,71,30,19}
};

void swap9(int **x) {
//   -----
  int j=random(12);
  switch (j) {
    case 0:
      swapRows(x,0,1); swapRows(x,7,8); swapCols(x,9,0,1); swapCols(x,9,7,8); break;
    case 1:
      swapRows(x,0,1); swapRows(x,7,8); swapCols(x,9,0,1); swapCols(x,9,7,8);
      swapRows(x,3,5); swapCols(x,9,3,5); break;
    case 2:
      swapRows(x,0,1); swapRows(x,7,8); swapCols(x,9,0,1); swapCols(x,9,7,8);
      swapRows(x,0,2); swapRows(x,6,8); swapCols(x,9,0,2); swapCols(x,9,6,8); break;
    case 3:
      swapRows(x,0,1); swapRows(x,7,8); swapCols(x,9,0,1); swapCols(x,9,7,8);
      swapRows(x,0,2); swapRows(x,6,8); swapCols(x,9,0,2); swapCols(x,9,6,8);
      swapRows(x,3,5); swapCols(x,9,3,5); break;
    case 4:
      swapRows(x,0,1); swapRows(x,7,8); swapCols(x,9,0,1); swapCols(x,9,7,8);
      swapRows(x,1,2); swapRows(x,6,7); swapCols(x,9,1,2); swapCols(x,9,6,7); break;
    case 5:
      swapRows(x,0,1); swapRows(x,7,8); swapCols(x,9,0,1); swapCols(x,9,7,8);
      swapRows(x,1,2); swapRows(x,6,7); swapCols(x,9,1,2); swapCols(x,9,6,7);
      swapRows(x,3,5); swapCols(x,9,3,5); break;
    case 6:
      swapRows(x,0,2); swapRows(x,6,8); swapCols(x,9,0,2); swapCols(x,9,6,8); break;
    case 7: // ultra
      swapRows(x,0,2); swapRows(x,6,8); swapCols(x,9,0,2); swapCols(x,9,6,8);
      swapRows(x,3,5); swapCols(x,9,3,5); break;
    case 8:
      swapRows(x,1,2); swapRows(x,6,7); swapCols(x,9,1,2); swapCols(x,9,6,7); break;
    case 9:
      swapRows(x,1,2); swapRows(x,6,7); swapCols(x,9,1,2); swapCols(x,9,6,7);
      swapRows(x,3,5); swapCols(x,9,3,5); break;
    case 10:
      swapRows(x,3,5); swapCols(x,9,3,5); break;
    default: break; // ultra
  }
} // swap9

bool makeSODLS(int **x, const int a) {
//   ---------
  if (a==9) {
    const int j=random(7), k=random(2); int i=k?0:80;
    for (int r=0; r<9; ++r) for (int c=0; c<9; ++c) x[r][c]=k?u9[j][i++]/9:u9[j][i--]/9; swap9(x);
  } else getRowShift2(x, a);
  return T;
} // makeSODLS

bool make9(int **x) {
//   -----
  const int j=random(7), k=random(2); int i=k?0:80;
  for (int r=0; r<9; ++r) for (int c=0; c<9; ++c) x[r][c]=k?u9[j][i++]+1:u9[j][i--]+1;
  swap9(x); return T;
} // make9

bool make3kq(int **x, const int n, int **xQ) {
//   -------
  const int bc=n/3; int **blocks=NULL;
  if (allocateRectangle(&blocks, 3, bc)) {
    for (int r=0; r<3; ++r) { blocks[r][0]=r; blocks[r][1]=5-r;
      for (int c=2; c<bc; ++c) blocks[r][c]=blocks[r][c-2]+6;
    }
    blocks[0][2]=7; blocks[1][0]=2; blocks[1][2]=6; blocks[2][0]=1;
    for (int R=0,xr=0; R<3; ++R,xr+=bc) for (int C=0,xc=0; C<3; ++C,xc+=bc)
      for (int r=0; r<bc; ++r) for (int c=0; c<bc; ++c)
        x[xr+r][xc+c]=blocks[R][xQ[r][c]]*n+blocks[C][xQ[c][r]]+1;
    freeRectangle(&blocks, 3); return T;
  }
  return F;
} // make3kq

// Any odd multiple of 3. Pandiagonal with 3x3 semimagic blocks.
bool make3kq3(int **x, const int n, int **rect) {
//   --------
  const int a=n/3; int **q=NULL;
  if (allocateRectangle(&q, n, n)) {
    for (int r=0, s=0; r<a; ++r, s+=3) {
      q[s][0]=rect[0][r]; q[s][1]=rect[1][r]; q[s][2]=rect[2][r];
      for (int c=3; c<n; ++c) q[s][c]=q[s][c-3];
      for (int t=s+1; t<(s+3); ++t) {
        q[t][0]=q[t-1][n-1]; for (int c=1; c<n; ++c) q[t][c]=q[t-1][c-1];
      }
    }
    for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) x[r][c]=q[r][c]*n+q[c][n-r-1]+1;
    freeRectangle(&q, 3); return T;
  }
  return F;
} // make3kq3

// Only for 27. Makes magic or pandiagonal with 9x9 magic blocks.
bool make3kqr(int **x, const int n, int **xQ, int **xR) {
//   --------
  const int bc=n/3; int **blocks=NULL;

  if (allocateRectangle(&blocks, 3, bc)) {
    for (int r=0; r<3; ++r) { blocks[r][0]=r; blocks[r][1]=5-r;
      for (int c=2; c<bc; ++c) blocks[r][c]=blocks[r][c-2]+6;
    }
    if (bc&1) { blocks[0][2]=7; blocks[1][0]=2; blocks[1][2]=6; blocks[2][0]=1; }
    for (int R=0,xr=0; R<3; ++R,xr+=bc) for (int C=0,xc=0; C<3; ++C,xc+=bc)
      for (int r=0; r<bc; ++r) for (int c=0; c<bc; ++c)
        x[xr+r][xc+c]=blocks[R][xQ[r][c]]*n+blocks[C][xR[r][c]]+1;
    freeRectangle(&blocks, 3); return T;
  }
  return F;
} // make3kqr

bool make6k(int **x, const int n, int **m6) { // m6: special magic
//   ------
  const int nd6=n/6; int **blocks=NULL, br0=nd6*nd6; const int bc=36;
  if (allocateRectangle(&blocks, br0, bc)) {
    for (int r=0; r<br0; ++r) { const int brpbr=br0+br0;
      blocks[r][0]=r+1; blocks[r][1]=brpbr-r;
      for (int c=2; c<bc; ++c) blocks[r][c]=blocks[r][c-2]+brpbr;
    }

    for (int R=0,br=0,xr=0; R<nd6; ++R,br+=nd6,xr+=6) {
      for (int C=0,xc=0; C<nd6; ++C,xc+=6) {
        for (int r=0; r<6; ++r) {
          for (int c=0; c<6; ++c) { x[xr+r][xc+c]=blocks[br+C][m6[r][c]]; }
        }
      }
    }
    freeRectangle(&blocks, br0); return T;
  }
  return F;
} // make6k

bool make4k(int **x, const int n, int **ap4) { // ap4: assoc or pan
//   ------
  const int nd4=n/4;
  int **blocks=NULL, br0=nd4*nd4; const int bc=16;
  if (allocateRectangle(&blocks, br0, bc)) {
    for (int r=0; r<br0; ++r) { const int brpbr=br0+br0;
      blocks[r][0]=r+1; blocks[r][1]=brpbr-r;
      for (int c=2; c<bc; ++c) blocks[r][c]=blocks[r][c-2]+brpbr;
    }

    for (int R=0,br=0,xr=0; R<nd4; ++R,br+=nd4,xr+=4) {
      for (int C=0,xc=0; C<nd4; ++C,xc+=4) {
        for (int r=0; r<4; ++r) {
          for (int c=0; c<4; ++c) { x[xr+r][xc+c]=blocks[br+C][ap4[r][c]]; }
        }
      }
    }
    freeRectangle(&blocks, br0); return T;
  }
  return F;
} // make4k

// 1 associative, 1 pandiagonal, 1 ultramagic, except 100 of 81.
bool makeAA(int **x, const int a, int **sodls) { // printf("makeAA\n");
//   ------
  const int n=a*a; int **blocks=NULL;
  if (allocateRectangle(&blocks, a, a)) {
    for (int r=0; r<a; ++r) for (int c=0; c<a; ++c) blocks[r][c]=sodls[r][c]*a+sodls[c][r];
    for (int r=0; r<a; ++r) {
      for (int c=0; c<(a-1); ++c) {
        bool action=T;
        while(action) {
          action=F;
          for (int d=c; d<a; ++d) {
            if (blocks[r][d]<blocks[r][c]) {
              action=T; int t=blocks[r][c]; blocks[r][c]=blocks[r][d]; blocks[r][d]=t;
            }
          }
        }
      }
    }
    for (int R=0,xr=0; R<a; ++R,xr+=a) for (int C=0,xc=0; C<a; ++C,xc+=a)
      for (int r=0; r<a; ++r) for (int c=0; c<a; ++c)
        x[xr+r][xc+c]=blocks[R][sodls[r][c]]*n+blocks[C][sodls[c][r]]+1;
    freeRectangle(&blocks, a); return T;
  }
  return F;
} // makeAA

// Pandiagonal, 100.
bool makeAB(int **x, const int a, const int b, int **rect) { // printf("makeAB\n");
//   ------
  const int n=a*b;
  for (int R=0,xr=0; R<a; ++R,xr+=b) for (int C=0,xc=0; C<a; ++C,xc+=b)
    for (int r=0; r<b; ++r) {
      int i=b-r*2; if (i<0) i+=b; int j=b-r*3; while (j<0) j+=b;
      for (int c=0; c<b; ++c) {
        x[xr+r][xc+c]=rect[R][(c+i)%b]*n+rect[C][(c+j)%b]+1;
      }
    }
  return T;
} // makeAB

bool getRectangle(int**x, const int a, const int b, magicType t) {
//   ------------
  t_iP restart; t_iPcIiP tabuSums;
  if (t==normalSemiMagic) {restart=restartSemi; tabuSums=tabuSumsSemi;
  } else if (t==Associative) { restart=restartSym; tabuSums=tabuSumsSym;
  } else if (t==normalMagic) { restart=restartMagic; tabuSums=tabuSumsMagic; }
  initialize(xRectangle, a, b, t);
  int violation=restart(xRectangle); tabuSums(xRectangle, 10000, &violation);
  if (violation==0) { int i=0;
    for (int r=0; r<a; ++r) for (int c=0; c<b; ++c) x[r][c]=xRectangle[i++]-1; return T;
  }
  return F;
} // getRectangle

bool getRectangleSA(int **x, const int a, const int b) { // semi, assoc
//   --------------
  return getRectangle(x, a, b, random(2)?Associative:normalSemiMagic);
} // getRectangleSA

int x6[][36]={
{0,34,33,32,1,5,29,7,27,8,10,24,23,22,14,15,19,12,17,13,20,21,16,18,6,25,9,26,28,11,30,4,2,3,31,35},
{0,22,27,33,16,7,28,6,34,13,20,4,11,5,12,26,30,21,31,15,3,23,9,24,18,32,10,2,29,14,17,25,19,8,1,35}
};
bool get6x6(int **x, const int j) {
//   ------
  const int k=random(2), l=random(3); int i=k?0:35;
  for (int r=0; r<6; ++r) for (int c=0; c<6; ++c) x[r][c]=k?x6[j][i++]:x6[j][i--];
  if (l) { swapRows(x,0,l); swapRows(x,5,5-l); swapCols(x,6,0,l); swapCols(x,6,5,5-l); }
  switch (random(8)) {
    case 1: swapRows(x,1,4); swapCols(x,6,1,4); break;
    case 2: swapRows(x,2,3); swapCols(x,6,2,3); break;
    case 3: swapRows(x,1,4); swapRows(x,2,3); swapCols(x,6,1,4); swapCols(x,6,2,3); break;
    case 4: swapRows(x,1,2); swapRows(x,3,4); swapCols(x,6,1,2); swapCols(x,6,3,4); break;
    case 5: swapRows(x,1,2); swapRows(x,3,4); swapCols(x,6,1,2); swapCols(x,6,3,4);
            swapRows(x,1,4); swapCols(x,6,1,4); break;
    case 6: swapRows(x,1,2); swapRows(x,3,4); swapCols(x,6,1,2); swapCols(x,6,3,4);
            swapRows(x,2,3); swapCols(x,6,2,3); break;
    case 7: swapRows(x,1,2); swapRows(x,3,4); swapCols(x,6,1,2); swapCols(x,6,3,4);
            swapRows(x,1,4); swapCols(x,6,1,4); swapRows(x,2,3); swapCols(x,6,2,3); break;
    default: break;
  }
  return T;
} // get6x6

bool get4x4(int **x) {
//   ------
  int loop=1000; initialize(xRectangle, 4, 4, Associative);
  while (loop--) {
    int violation=restartSym(xRectangle); tabuSumsSym(xRectangle, 1000, &violation);
    if (violation==0) { int i=0;
      for (int r=0; r<4; ++r) for (int c=0; c<4; ++c) x[r][c]=xRectangle[i++]-1;
      if (random(2)) XformA_D(x,4); return T;
    }
  }
  return F;
} // get4x4

bool makeSquare(int **, const int, int *, int *, int *, FILE *, FILE *);

bool makePrime(int **x, const int n, int *bn, int *borders) {
//   ---------
  const int w=n-2; int **y=NULL; bool ok=F; int bb=0;
  if (allocateRectangle(&y, w, w)) {
    if (makeSquare(y,w,bn,borders,&bb,NULL,NULL)) { const int k=n+n-2; ok=T;
      for (int r=0; r<w; ++r) for (int c=0; c<w; ++c) x[r+1][c+1]=y[r][c]+k;
      makeOddBorder(x, bSquare, n); if (n==7) { *bn=3; *borders=2; } else ++*borders;
    }
    freeRectangle(&y, w);
  }
  return ok;
} // makePrime

// Magic or pandiagonal with 9x9 magic blocks.
bool make27(int **x, const int n, int *bn) {
//   -------
  int b=n/3, **X; bool ok=F;
  if (allocateRectangle(&X,b,b)) {
    int **R; make9(X);
    if (allocateRectangle(&R,b,b)) {
      for (int r=0; r<b; ++r) for (int c=0; c<b; ++c) {
        R[r][c]=(X[r][c]-1)%b; X[r][c]=(X[r][c]-1)/b;
      }
      ok=make3kqr(x,n,X,R); *bn=b; freeRectangle(&R, b);
    }
    freeRectangle(&X, b);
  }
  return ok;
} // make27

bool makeUnequalA(int **x, const int n, const int a, int **xa) {
//   ------------
  int **xb; bool ok=F; const int aa=a*a, b=n/a;
  if (allocateRectangle(&xb,b,b)) {
    if (b==3) ok=getRectangle(xb, b, b, Associative);
    else {
      int **t;
      if (allocateRectangle(&t,b,b)) {
        ok=makeSODLS(t,b);
        for (int r=0; r<b; ++r) for (int c=0; c<b; ++c) xb[r][c]=b*t[r][c]+t[c][r];
        freeRectangle(&t, b);
      }
    }
    if (ok) {
      for (int R=0,xr=0; R<a; ++R,xr+=b) {
        for (int C=0,xc=0; C<a; ++C,xc+=b) {
          for (int r=0; r<b; ++r) {
            for (int c=0; c<b; ++c) { x[xr+r][xc+c]=aa*xb[r][c]+xa[R][C]+1; }
          }
        }
      }
    }
    freeRectangle(&xb, b);
  }
  return ok;
} // makeUnequalA

bool makeUnequalB(int **x, const int n, int **xa) { 
//   ------------
  int **xb; bool ok=F; const int a=n/3, aa=a*a, b=3, ab=a*b;
  if (allocateRectangle(&xb,b,b)) {
    ok=getRectangle(xb, b, b, Associative); int **dist;

    if (ok&&allocateRectangle(&dist,ab,ab)) {
      int Qb[b][b], Rb[b][b];
      for (int r=0; r<b; ++r) for (int c=0; c<b; ++c) {
        Qb[r][c]=xb[r][c]/b; Rb[r][c]=xb[r][c]%b;
      }
      int i=0;
      for (int r=0; r<ab; ++r) for (int c=0; c<ab; ++c) dist[r][c]=++i;
      for (int i=3; i<ab; i+=6) { swapRows(dist,i,i+2); swapCols(dist,ab,i,i+2); }

      for (int R=0,xr=0; R<a; ++R,xr+=b) {
        for (int C=0,xc=0; C<a; ++C,xc+=b) {
          int d=xa[R][C], dr=b*(d/a), dc=b*(d%a);
          for (int r=0; r<b; ++r) {
            for (int c=0; c<b; ++c) { x[xr+r][xc+c]=dist[dr+Qb[r][c]][dc+Rb[r][c]]; }
          }
        }
      }
      freeRectangle(&dist, ab);
    }
    freeRectangle(&xb, b);
  }
  return ok;
} // makeUnequalB

bool makeOdd(int **x, const int n, int *bn, int *borders) { // website
//   -------
  if (n==1) { x[0][0]=1; *bn=1; return T; }
  if (isPrime(n)) return makePrime(x, n, bn, borders);
  if ((n==9)&&random(2)) { *bn=3; return make9(x); }
  if ((n==27)&&random(2)) return make27(x, n, bn);
  int a, b, **y; bool ok=F;
  if (getFactors(n,&a,&b)&&allocateRectangle(&y,n,n)) {
    if (b==3) {
      if (((a%3)==0)||!random(5)) { *bn=3;
        ok=getRectangleSA(y,b,a)&&make3kq3(x,n,y);
      } else {
        getRowShift2(y, a); ok=make3kq(x,n,y); *bn=a; 
      }
    } else if ((a==b)&&((a==9)||random(2))) { makeSODLS(y,a); ok=makeAA(x,a,y); *bn=b;
    } else if ((b%3)!=0) if (getRectangleSA(y,a,b)) { ok=makeAB(x,a,b,y); *bn=b; }
    freeRectangle(&y,n);
  }
  return ok;
} // makeOdd

void fillBBcorners(int **x, const int g, const int pair, const int q, int rs, int re, const int cs0, const int ce0) {
//   -------------
  const int qq=q*q; int cs=cs0, ce=ce0, t=x[rs][cs]; t=qq*(t-1)+1;
  for (int i=0; i<q; ++i) {
    for (int k=0; k<q; ++k) { x[rs][cs]=t; x[re][ce]=pair-t++; cs+=g; ce+=g; } rs+=g; cs=cs0; re+=g; ce=ce0;
  }
} // fillBBcorners

void fillBBrows(int **x, const int g, const int pair, const int q, int r, const int cs0, const int ce0) {
//   ----------
  const int qq=q*q; int cs=cs0, ce=ce0, t=x[r][cs]; t=qq*(t-1)+1;
  for (int i=0; i<q; ++i) {
    for (int k=0; k<q; ++k) { x[r][cs]=t; x[r][ce]=pair-t++; cs+=g; ce+=g; } r+=g; cs=cs0; ce=ce0;
  }
} // fillBBrows

void fillBBcols(int **x, const int g, const int pair, const int q, const int c0, int rs, int re) {
//   ----------
  const int qq=q*q; int c=c0, t=x[rs][c]; t=qq*(t-1)+1;
  for (int i=0; i<q; ++i) {
    for (int k=0; k<q; ++k) { x[rs][c]=t; x[re][c]=pair-t++; c+=g; } c=c0; rs+=g; re+=g;
  }
} // fillBBcols

bool makeBB(int **x, const int n, int *bn, int *bb) {
//   ------
  if ((*bn!=4)&(*bn!=6)) return F; int q=2, q2, nbb, loop=(n<12)?0:n/6-1, qs[20], nq=0;
  while (loop--) {
    q2=q+q; nbb=n-q2; const int bs=nbb/q; if (((nbb%q)==0)&((bs%*bn)==0)) qs[nq++]=q; ++q;
  }
  if (nq==0) return F; q=qs[random(nq)]; q2=q+q;
  const int m=n-1, g=n/q, f=g-1, j=n-q2, e=j/q, qq=q*q, add=q2*(n-q), pair=n*n+1;
  if (*bn==4) { if (get4x4(bSquare)) make4k(x,j,bSquare); else return F; }
  else { if (get6x6(bSquare,0)) make6k(x,j,bSquare); else return F; }
  for (int r=0; r<j; ++r) for (int c=0; c<j; ++c) x[r][c]+=add;

  // Distribute sub-squares leaving space for borders.
  int d=q2-1, ds=j-1, de=ds-e;
  for (int i=0; i<q; ++i) {
    for (int r=ds; r>de; --r) for (int c=0; c<j; ++c) x[r+d][c]=x[r][c]; d-=2; ds-=e; de-=e;
  }
  d=q2-1, ds=j-1, de=ds-e;
  for (int i=0; i<q; ++i) {
    for (int r=1; r<m; ++r) for (int c=ds; c>de; --c) x[r][c+d]=x[r][c]; d-=2; ds-=e; de-=e;
  }
  makeEvenBorder(x,bSquare,g);
  // Corners -TL BR
  if (x[0][0]<x[f][f]) fillBBcorners(x,g,pair,q,0,f,0,f); else fillBBcorners(x,g,pair,q,f,0,f,0);
  // Corners -TR BL
  if (x[0][f]<x[f][0]) fillBBcorners(x,g,pair,q,0,f,f,0); else fillBBcorners(x,g,pair,q,f,0,0,f);
  for (int r=1; r<f; ++r) // Borders left, right
    if (x[r][0]<x[r][f]) fillBBrows(x,g,pair,q,r,0,f); else fillBBrows(x,g,pair,q,r,f,0);
  for (int c=1; c<f; ++c) // Borders top, bottom
    if (x[0][c]<x[f][c]) fillBBcols(x,g,pair,q,c,0,f); else fillBBcols(x,g,pair,q,c,f,0);
  *bb=q; return T;
} // makeBB

bool makeEven(int **x, const int n, int *bn, int *borders, int *bb) {
//   --------
  bool ok=F; const int f=*bn;
  if ((*bb==0)&&(*borders==0)&&random(2)) ok=makeBB(x,n,bn,bb);
  if (!ok) {
    *bb=-1;
    if ((n%f)==0) {
      if (f==4) { if (get4x4(bSquare)) { make4k(x,n,bSquare); ok=T; } }
      else { if (get6x6(bSquare,0)) { make6k(x,n,bSquare); ok=T; } }
    } else {
      int w=n-2, **y; 
      if (allocateRectangle(&y,w,w)) {
        if (makeSquare(y,w,bn,borders,bb,NULL,NULL)) { const int k=n+w;
          for (int r=0; r<w; ++r) for (int c=0; c<w; ++c) x[r+1][c+1]=y[r][c]+k;
          makeEvenBorder(x,bSquare,n); ++*borders; ok=T;
        }
        freeRectangle(&y, w);
      }
    }
  }
  return ok;
} // makeEven

bool makeEvenU(int **x, const int n, int *bn, int *borders) { // Unequal magic sum blocks.
//   ---------
  const int f=*bn, a=n/f; bool ok=F; *borders=-1;
  
  if (a==4) ok=get4x4(bSquare);
  else if (a==6) ok=(f==3)?get6x6(bSquare,random(2)):getRectangle(bSquare,a,a,normalMagic);
  if (ok) {
    if ((n==18)||((n==12)&&random(2))) ok=makeUnequalB(x,n,bSquare);
    else ok=makeUnequalA(x,n,a,bSquare);
  } 
  return ok;
} // makeEvenU

bool makeSquare(int **x, const int n, int *bn, int *borders, int *bb, FILE *wtfp, FILE *whfp) {
//   ----------
  bool out=(wtfp!=NULL), ok=F;
  if (n&1) {
    ok=makeOdd(x,n,bn,borders);
  } else {
    if (*bn==0) {
      *bn=(n==4)||random(2)?4:6;
      if ((n>11)&&!random(3)) {
        if ((n&3)==2) {
          if ((n%6)==0) { const int a=n/6; if ((a==3)|(a==9)|((a%3)!=0)) *bn=n/6; }
        } else { const int a=n/4; if ((a==3)|(a==9)|((a&1)&((a%3)!=0))) *bn=n/4; }
      }
    }
    ok=(*bn&1)?makeEvenU(x,n,bn,borders):makeEven(x,n,bn,borders,bb);
  }
  if (ok&&out&&pushSquare(x,n,stack)) {
    ok=txtPrint(x,n,wtfp)&&htmlPrint(x,n,*bn,*borders,*bb,whfp); writeError=!ok;
  }
  return ok;
} // makeSquare
//=================================================== main ==============================================

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

int main() {
//  ----
  int n=0; bool inputSize=T; seed_rand(); bool ok=T;

  if (openDir()) { 
    do {     
      if (inputSize) { printf("\nSquare order? "); n=getSize(); }
      if (validOrder(n)) {
        if (allocateStore(n)) {
          int howMany=1; const int limit=(n==3)?8:maxSquares; ok=F;
          if (n>1) { printf("\nHow many? "); howMany=getSize(); } if (howMany<1) howMany=1;
          else if (howMany>limit) { printf("Limit is %d.\n", limit); howMany=limit; }
	        FILE *wtfp=NULL;
	        if ((wtfp=open_File(n,"txt"))!=NULL) {
            FILE *whfp=NULL;
	          if ((whfp=open_File(n,"html"))!=NULL) {
              writeError=F;
              if (writeHead(n, whfp)) {
                int loop=10000; squaresMade=0;
                while (loop--) {
                  int bn=0, borders=0, bb=0;
                  makeSquare(xSquare,n,&bn,&borders,&bb,wtfp, whfp);
                  if (writeError|(squaresMade==howMany)) break;
                }
                if (!writeError) writeError=!writeFoot(whfp);
              } else writeError=T;
              if (writeError) { perror("\a\nError writing file"); break; } else ok=T;
              fclose(whfp);
            }
            fclose(wtfp);
	        }
        } else {
	        printf("\a\nERROR: Storage allocation failed.\n"); ok= F; break;
        }
      } // if (validOrder(n))
      if (ok) {
        printf("\nContinue? y (yes) or n (no) or the square order: ");
        if (getYorOrder(&n)) inputSize=(n==0); else break;
      } else break;
    } while (T);
  } // if (openDir())

  freeStore(n); printf("\nPress a key to close the console.");
  while (!_kbhit()) Sleep(250); return ok?EXIT_SUCCESS:EXIT_FAILURE;
} // main