/*
 *  File:    SODLS.cpp
 *  Author:  S Harry White
 *  Created: 2015-12-08
 *  Updated: 2018-12-20
 *  Updated: 2021-01-14
 *    Add C values 12, 14, 16 for method makeFromKplusC. Order 22,623 required 14 or 16.
 *  Updated: 2021-04-12
 *    Tidy code. Add comments.
 *  Updated: 2021-06-19
 *    Increase application of x k's+1.
 *    Add option to output SODLS only, magic square only, or both SODLS and magic square.
 *  Updated: 2021-07-07
 *    Include "diagonal with center 4" SODLS examples to replace code for making them.
 *    Add examples made from Latin squares to previous ones made by backtracking
 *  Updated: 2021-07-10
 *    Replace examples made from Latin squares with code.
 *  Updated: 2021-07-15
 *    Re-name from MagicSquaresSODLS.
 *    Add code for center 4, 5, 7, 8, 9, 10, 11, 12 sub-square SODLS.
 *  Updated: 2021-07-26
 *    Dispense with making copies of center sub-square SODLS. It was not really an optimization.
 *  Updated: 2022-01-21
 *    Replace & with && where necessary.
 * Updated: 2023-01-26
 *   Change outputFile from bufSize to outSize.
 */

/*
 * Makes self-orthogonal diagonal Latin squares, (SODLS), and their magic squares for all
 * valid orders, (i.e., not order 2, 3, or 6). Makes a minimum of 1 SODLS for any order,
 * many SODLS for some orders. See:
 *
 *                    https://budshaw.ca/SODLS.html
 *                    https://budshaw.ca/addenda/SODLSmethods.html
 *                    https://budshaw.ca/addenda/SODLSnotes.html
 *
 * The magic squares M are made from the SODLS Q and its transpose as:
 *
 *                M[r][c]=N*Q[r][c]+Q[c][r]+1.
 *
 * Note: Associative, pandiagonal, and ultramagic below refer to SODLS that make these
 *       kinds of magic square.
 *
 * There are no pandiagonal Latin squares, (Knut Vik designs), for even orders or odd orders
 * that are a multiple of 3. For these orders, the SODLS that make pandiagonal and ultramagic
 * squares are only weakly pandiagonal, (SOWPLS), meaning that the sum of the n numbers in
 * each of the 2n diagonals is the same.
 *
 * Odd Order
 * ---------
 * Diagonal, centers 5, 7, 9, 11 are made like orthogonal order 18 by Natalia Makarova,
 * with reference to Colbourn and Dinitz.
 *
 * Odd Order, not multiple of 3:
 * -----------------------------
 *  Note: Pandiagonal and ultramagic here refer to Knut Vik designs. 
 *
 *  Pandiagonal, ultramagic, and associative are easily made:
 *
 *          . fill row 0: for pandiagonal use natural order
 *                        for ultramagic use (n-column) % n
 *          . start the next rows with the last 2 of the previous row
 *
 *                  pan              ultra            assoc
 *               0 1 2 3 4         0 4 3 2 1        0 2 3 4 1
 *               3 4 0 1 2         2 1 0 4 3        1 3 4 0 2
 *               1 2 3 4 0         4 3 2 1 0        4 1 2 3 0
 *               4 0 1 2 3         1 0 4 3 2        2 4 0 1 3
 *               2 3 4 0 1         3 2 1 0 4        3 0 1 2 4
 *
 * The pandiagonal method is borrowed from Inder Taneja.
 * The ultramagic method is based on the formula Q[r][c] = (2r-c)%n for prime orders ≥ 5 by
 * Cheng-Xu Xu and Zhun-Wei Lu.
 * Other associative are made by swapping rows/columns of the ultramagic.
 * 
 * Small associative squares are also made by a backtracking procedure.
 * Orders 11 and 13 associative made by backtracking are included in the program.
 * 
 * Non-prime order squares such as 25, 35, 49, 55 are sometimes made as composites
 * of smaller ones.
 * 
 * Another method is x ks+1 such as 17 (4 4s+1), 29 (4 7s+1) (7 4s+1), 49 (4 12s+1) (12 4s+1)
 * similar to an order 21 by Mitsutoshi Nakamura.
 *
 * Odd order, multiple of 3
 * -------------------------
 * Note: Pandiagonal and ultramagic here refer to weakly pandiagonal.
 * 
 * For order 9, an associative is made by backtracking, a pandiagonal is made based on one by
 * Inder Taneja, and an ultramagic made by backtracking is included in the program.
 *
 * An order 15 associative made by backtracking is included in the program.
 * An order 27 pandiagonal is made based on one by Cheng-Xu Xu and Zhun-Wei Lu.
 * 
 * Other methods, (extending concepts from squares by Mitsutoshi Nakamura):
 * 
 *     x 5s+4 such as 39 (7 5s+4), 69 (13 5s+4), 129 (25 5s+4), 159 (31 5s+4)
 *     x ks+1 such as 21 (5 4s+1), 33 (4 8s+1), 183 (13 14s+1), 633 (8 79s+1)
 *     x ks+k/2 such as 27 (5 5s+2), 93 (5 17s+8), 123 (9 13s+6), 753 (5 137s+68)
 *     x ks+c for c=4, 8, 10, 12, 14, 16, such as: 303 (13 23s+4), 543 (7 77s+4) (13 41s+10),
 *        843 (17 49s+10), 1983 (25 79s+8), 6663 (11 605s+8), 22623 (23 983s+14) (37 611s+16) 
 * 
 * Some squares, such as 45, 63, 75, 81 are made as composites of smaller ones.
 *
 * Even Order
 * ----------
 * Diagonal, centers 4, 8, 10, 12 are made like orthogonal order 18 by Natalia Makarova,
 * with reference to Colbourn and Dinitz. 
 *
 * Doubly Even Order
 * -----------------
 * Note: Pandiagonal and ultramagic here refer to weakly pandiagonal.
 * 
 * Small order associative squares are made by a backtracking procedure.
 * 
 * An order 12 associative made by backtracking is included in the program.
 * Order 24 is made as 4 5s+4 similar to a 24 magic by Mitsutoshi Nakamura.
 * An order 8 ultramagic and an order 16 associative made by backtracking are included for variety.
 * 
 * Associative squares are sometimes converted to pandiagonal squares by Planck's A-D method.
 * 
 * Some squares, such as 16, 20, 28, 32 are made as composites of smaller ones.
 * 
 * Method x ks+1 such as 36 (5 7s+1), 56 (5 11s+1), 64 (7 9s+1) (9 7s+1) is also used.
 *
 * Diagonal, center 4 squares are made for orders 16, 20, 24 and 28. 
 *
 * Singly Even Order
 * -----------------
 * Note: There are no associative or pandiagonal Latin squares of singly-even order.
 * 
 * An order 10 magic, made by backtracking, is included in the program. Two order 14 magic,
 * with center 4x4, are included. One is adapted from an order 14 by Inder Taneja based on
 * Frank E. Bennett, Beiliang Du and Hantao Zhang; the other order 14 was made by backtracking.
 * 
 * Diagonal, center 4 for orders 18, 22, 26 and 30 are included.
 * These are similar to an order 18 by Mitsutoshi Nakamura. 
 * 
 * Other methods, (extending concepts from Mitsutoshi Nakamura's squares):
 * 
 *     x 5s+4 such as 54 (10 5s+4), 94 (18 5s+4), 174 (34 5s+4), 214 (42 5s+4)
 *     x 8s+2 such as 34 (4 8s+2), 42 (5 8s+2), 58 (7 8s+2), 66 (8 8s+2)
 *     x ks+1 such as 46 (5 9s+1), 78 (7 11s+1), 86 (5 17s+1), 118 (9 13s+1)
 *     x ks+k/2 such as 22 (4 5s+2), 38 (5 7s+3), 446 (23 19s+9), 838 (19 43s+21)
 *     x ks+c for c=4, 8, 10, 12, 14, 16, such as: 158 (7 22s+4) (37 4s+10),
 *        278 (67 4s+10) (19 14s+12), 718 (7 102s+4) (71 10s+8) (59 12s+10),
 *        758 (13 58s+4) (25 30s+8) (11 68s+10) (17 44s+10) (31 24s+14) (53 14s+16)
 * 
 * x ks+k/2 are similar to an order 22 by Mitsutoshi Nakamura.
 * 
 * Some squares, such as 50, 70, 90, 98 are made as composites of smaller ones.
 * 
 */

#include "stdafx.h"
#include <conio.h>
#include <direct.h>
#include <errno.h>
#include <fcntl.h>
#include <io.h>
#include <math.h>
#include <share.h>
#include <sys\stat.h>
#include <time.h>
#include <Windows.h>

const int ptrSize=sizeof(void *); int **Msquare=NULL, **Qsquare=NULL, **Tsquare=NULL, allocatedSize=0;

const bool F=false, T=true;
bool **rUsed=NULL, **cUsed=NULL, **qrUsed=NULL, *bdUsed=NULL, *fdUsed=NULL, *numberUsed=NULL;

enum sqKind { none, magic, assoc, pan, ultra };
//char *kindName[]={ "none", "magic", "assoc", "pan", "ultra" };
typedef int kindArray[ultra+1];
//================================================== store ==============================================

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

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

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

void freeInts(int **line, const int size) { free(*line); *line=NULL; }
//   -------

void freeBools(bool **line, const int size) { free(*line); *line=NULL; }
//   ---------

void freePtrs(void ***line, const int size) { free(*line); *line=NULL; }
//   --------

void freeStore(const int size) {
//   ---------
  freeIntSquare(&Msquare,size); freeIntSquare(&Qsquare,size); freeIntSquare(&Tsquare,size);
  freeBoolSquare(&rUsed,size); freeBoolSquare(&cUsed,size);
  freeBoolSquare(&qrUsed,size); freeBools(&bdUsed,size);
  freeBools(&fdUsed,size); freeBools(&numberUsed,size*size+1); allocatedSize=0;
} // freeStore

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

bool allocateBools(bool **line, const int size) {
//   -------------
  *line=(bool*)malloc(size*sizeof(bool)); if (*line==NULL) return reportError(allocFail); return T;
} // allocateBools

bool allocatePtrs(void ***line, const int size) {
//   ------------
  *line=(void **)malloc(size*ptrSize); if (*line==NULL) return reportError(allocFail); return T;
} // allocatePtrs

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

bool allocateIntSquare(int*** square, const int size) {
//   -----------------
  bool ok; *square=(int**)malloc(size*ptrSize); 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) freeIntSquare(square,numAllocated);
  }
  if (!ok) return reportError(allocFail); return T;
} // allocateIntSquare

bool allocateStore(const int size) {
//   -------------
  bool ok=T;
  if (size>allocatedSize) {
    freeStore(allocatedSize); allocatedSize=size;
    if (ok=allocateIntSquare(&Msquare,size))
      if (ok=allocateIntSquare(&Qsquare,size))
        if (ok=allocateIntSquare(&Tsquare,size))
          if (ok=allocateBoolSquare(&rUsed,size))
            if (ok=allocateBoolSquare(&cUsed,size))
              if (ok=allocateBoolSquare(&qrUsed,size))
                if (ok=allocateBools(&bdUsed,size))
                  if (ok=allocateBools(&fdUsed,size)) ok=allocateBools(&numberUsed,size*size+1);
    if (!ok) { bell=T; printf("allocateStore %d\n",size); reportError(allocFail); freeStore(size); }
  }
  return ok;
} // allocateStore
//============================================ input ================================================

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

bool getY() {
//   ----
  int c; do { c=getchar(); } while ((c==' ')|(c=='\t')|(c=='\n'));
  clearLine(c); return (c=='Y')|(c=='y');
} // getY

bool validOrder(const int n, kindArray numSq) {
//   ----------
  if (n<=0) { bell=T; return reportError("N must be a positive integer"); }
  if ((n==2)|(n==3)|(n==6)) {
    if (numSq!=NULL) ++numSq[none];
    char s[50]; sprintf_s(s, 50, "There is no %dx%d SODLS", n, n); return reportError(s);
  };
  return T;
}// validOrder

bool checkRange(int *p, int *q) {
//   ----------
  if (*p<0) return validOrder(*p,NULL); if (*q<0) return validOrder(*q,NULL);
  if (*q<*p) { const int t=*p; *p=*q; *q=t; } if (*p==0) *p=*q; return T;
} // checkRange

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

bool getInts(int c, int *p, int *q) {
//   -------
  bool chs=F, result=F; *p=*q=0; if (c==' ') do { c=getchar(); } while ((c==' ')|(c=='\t'));
  if (c=='-') { chs=T; c=getchar(); }
  if ( ('0'<=c)&(c<='9') ) {
    int i=c-'0'; while ( ('0'<=(c=getchar()))&&(c<='9') ) i=i*10+c-'0';
    if (chs) i=-i; *p=i; chs=F;
    if ((c==' ')|(c=='\t')) {
      do { c=getchar(); } while ((c==' ')|(c=='\t')); if (c=='-') { chs=T; c=getchar(); }
      if ( ('0'<=c)&(c<='9') ) {
        int i=c-'0'; while ( ('0'<=(c=getchar()))&&(c<='9') ) i=i*10+c-'0'; if (chs) i=-i; *q=i;
      }
    }
    result=checkRange(p,q);
  }   
  clearLine(c); return result;
} // getInts

bool getNorOrders(int *p, int *q) {
//   ------------
  int c; do { c=getchar(); } while ((c==' ')|(c=='\t')|(c=='\n') );
  if ( (c=='N')|(c=='n') ) return F; return getInts(c,p,q);
} // getNorOrders
//============================================== check squares =============================================

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

bool isPandiagonal(int **x, const int n, const int magicSum) {
//   -------------
  if ((n&3)==2) return F; const int m=n-1;
  for (int i=0; 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][m-cModn]; }
    if ((sumYX!=magicSum)|(sumXY!=magicSum)) return F;
  }
  return T;
} // isPandiagonal

bool isMagic(int **M, const int n, const int nn, const int magicSum) {
//   -------
  const int m=n-1; int sumX, sumY, sumXY=0, sumYX=0; for (int i=0; i<=nn; i++) numberUsed[i]=F;
  for (int i=0; i<n; i++) {
    sumX=0; sumY=0; for (int j=0; j<n; j++) { const int t=M[i][j]; numberUsed[t]=T; sumX+=t; sumY+=M[j][i]; }
    if ((sumX!=magicSum)|(sumY!=magicSum)) { return F; } sumXY+=M[i][m-i]; sumYX+=M[i][i];
  }
  for (int i=1; i<=nn; i++) if (!numberUsed[i]) return F; return ((sumXY==magicSum)&(sumYX==magicSum));
} // isMagic

bool isUniqueQ(int **Q, const int n) {
//   ---------
  const int m=n-1;
  for (int i=0; i<n; ++i) {
    for (int j=0; j<n; ++j) { rUsed[i][j]=F; cUsed[i][j]=F; qrUsed[i][j]=F; } bdUsed[i]=F; fdUsed[i]=F;
  }
  for (int i=0; i<n; ++i) {
    for (int j=0; j<n; ++j) {
      const int q=Q[i][j], r=Q[j][i];
      if (rUsed[i][q]) { printf("rUsed[%d][%d]\n",i,q); return F; } rUsed[i][q]=T;
      if (cUsed[j][q]) { printf("cUsed[%d][%d]\n",j,q); return F; } cUsed[j][q]=T;
      if (qrUsed[q][r]) { printf("qrUsed[%d][%d]\n",q,r); return F; } qrUsed[q][r]=T; 
    }
    const int b=Q[i][i], f=Q[i][m-i];
    if (bdUsed[b]) { printf("bdUsed[%d]\n",b); return F; } bdUsed[b]=T;
    if (fdUsed[f]) { printf("fdUsed[%d]\n",f); return F; } fdUsed[f]=T;
  }
  return T;
} // isUniqueQ

bool isCorrect(int **M, int **Q, const int n, kindArray numSq) {
//   ---------
  const int nn=n*n, S2=nn+1, magicSum=n&1?n*(S2/2):(n/2)*S2; sqKind kind=none; 
  if (isUniqueQ(Q,n)&&isMagic(M,n,nn,magicSum)) {
    kind=magic;
    if (n!=1) {
      if (isAssociative(M,n,S2)) kind=isPandiagonal(M,n,magicSum)?ultra:assoc;
      else if (isPandiagonal(M,n,magicSum)) kind=pan;
    }
  }
  ++numSq[kind];
  printf("Square is %s.\n", kind==magic?"magic":kind==assoc?"associative":kind==pan?"pandiagonal":
                            kind==ultra?"ultramagic":"NOT magic");
  return kind!=none;
} // isCorrect
//=============================================== out ================================================

const int bufSize=1024; char outputDir[bufSize];
bool openDir() {
//   -------
  char buf[bufSize], msg[bufSize]; const char *baseName="SODLS"; strcpy_s(buf,bufSize,baseName);
  int sub=0; 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) {
    printf("Out files are in folder %s\n\n",buf); sprintf_s(outputDir,bufSize,"%s",buf); return T;
  }
  sprintf_s(msg,bufSize,"Can't open folder %s",buf); perror(msg); return F;
} // openDir

const int outSize=bufSize+50; char outputFile[outSize]; const int outMaxN=1000;
int openOutput(const int n) {
//  ----------
  const int baseSize=bufSize+25; char buf[outSize], baseName[baseSize];
  int wfd=-1; int sub=0; sprintf_s(baseName,baseSize,"%d",n); sprintf_s(buf,outSize,"%s.txt",baseName);
  do {
    if (_access_s(buf, 00)==ENOENT) break; sprintf_s(buf, outSize, "%s_%d.txt", baseName, ++sub);
  } while (T);
  if (_sopen_s(&wfd, buf, _O_CREAT|_O_WRONLY, _SH_DENYNO, _S_IWRITE)==0) {
    if (n<=outMaxN) printf(".. out file %s\n", buf); sprintf_s(outputFile,outSize,"%s",buf);
  } else {
    char msg[outSize+50]; sprintf_s(msg, outSize+50, "\a\nCan't open for write %s", buf); perror(msg);
  }
  return wfd;
} // openOutput

int fieldWidth(int v) {
//  ----------
  if (v<10000) if (v<100) return v<10?1:2; else return v<1000?3:4; return v<100000?5:v<1000000?6:7;
} // fieldWidth

const int writeSize=8000100; char outBuffer[writeSize];
void printSq1(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=1, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 2, "%1d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 3,"%2d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq1

void printSq2(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=2, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 3, "%2d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 4, "%3d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq2

void printSq3(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=3, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 4, "%3d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 5, "%4d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq3

void printSq4(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=4, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 5, "%4d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 6, "%5d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq4

void printSq5(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=5, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 6, "%5d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 7, "%6d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq5

void printSq6(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=6, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 7, "%6d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 8, "%7d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq6

void printSq7(int **S, const int n) {
//   --------
  char *s=outBuffer; const int FW0=7, FW=FW0+1;
  for (int i=0; i<n; i++) { sprintf_s(s, 8, "%7d", S[i][0]); s+=FW0;
    for (int j=1; j<n; j++) { sprintf_s(s, 9, "%8d", S[i][j]); s+=FW; } *s++='\n';
  } *s++='\n';
} // printSq7

typedef void (*t_printSq)(int **, const int);
static t_printSq printSq[]={ NULL, printSq1, printSq2, printSq3, printSq4, printSq5, printSq6, printSq7 };

bool printSquare(int **S, const int n, const int wfd, const bool square) {
//   -----------
  const int nn=n*n, FW0=fieldWidth(square?nn:n-1), FW=FW0+1, squareBytes=nn*FW+1;
  printSq[FW0](S,n); if (_write(wfd, outBuffer, squareBytes)==squareBytes) return T;
  printf("\aOutput error.\n"); return F;
} // printSquare

bool outputError, outMsg;
bool outputSquares(int **M, int **Q, const int n, const int wfd, kindArray numSq, const int out) {
//   -------------
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) M[r][c]=n*Q[r][c]+Q[c][r]+1;
  if (isCorrect(M,Q,n,numSq)) {
    bool ok=T;
    switch (out) {
      case 0: if (!outMsg) { printf("Output is limited to order %d.\n",outMaxN); outMsg=T; } break;
      case 1: ok=printSquare(Q,n,wfd,F); break;
      case 2: ok=printSquare(M,n,wfd,T); break;
      case 3: ok=printSquare(Q,n,wfd,F)&&printSquare(M,n,wfd,T); break;
      default: ok=F; break;
    }
    if (!ok) { printf("Output error.\n"); outputError= T; } putchar('\n'); return ok;
  }
  printf("\aProgram Error: fail n = %d\n",n); return F;
} // outputSquares
//============================================ common routines =============================================

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

enum doHow { noHow, even, odd, _5p4, _8p2, kp1, kpf, kpc, cn4, cn5, cn6, cn7, cn8, cn9, cnA, cnB, cnC };
//char *howName[]={ "noHow","even","odd","_5p4","_8p2","kp1","kpf","kpc",
//                  "cn4","cn5","cn6","cn7","cn8","cn9","cnA","cnB","cnC" };

#define not236(k) (((k)==1)|((k)==4)|((k)==5)|((k)>6))

#define not1236(k) (((k)==4)|((k)==5)|((k)>6))

#define not12356(k) (((k)==4)|((k)>6))

#define not1236d3(k) (((k)%3!=0)&(((k)==4)|((k)==5)|((k)>6)))

#define not1236odd(k) (((k)&1)&(((k)==4)|((k)==5)|((k)>6)))

doHow canDo(const int n) {
//    -----
  bool Odd=(n&1), Odd3=Odd&((n%3)==0), Seven=(n&3)==2;
  doHow how[50]; const int nd4=n/4, m=n-1; int num=0;
  if (Odd3||Seven||(random(Odd?3:5)==0)) for (int a=4; a<=nd4; ++a) { if (a==6) ++a;
    if (((m%a)==0)&not12356(m/a)) { how[num++]=kp1; break;}
  }
  if (n<7) { how[num++]=(doHow)((n==4)?even:(n==1)|(n==5)?odd:none); }
  else if (Odd) {
    if ((n>16)&(n<46)&&(random(3)==0)) {
      if (n<23) how[num++]=cn5;
      else { bool b=randomBool();
        if (n<29) how[num++]=b?cn5:cn7;
        else if (n<35) { const int b3=random(3); how[num++]=(b3==0)?cn5:(b3==1)?cn7:cn9; }
        else if (n==35) { const int b4=random(4); how[num++]=(b4==0)?cn5:(b4==1)?cn7:(b4==2)?cn9:cnB; }
        else if (n==37) { const int b3=random(3); how[num++]=(b3==0)?cn7:(b3==1)?cn9:cnB; }
        else if (n==39) how[num++]=b?cn9:cnB; else how[num++]=cnB;
      }
    }
    if (Odd3) { // multiple of 3
      if ((((n%9)==0)&not236(n/9))||(((n%15)==0)&not236(n/15))) how[num++]=odd;
      else if (n>20) {
        if (((n%27)==0)&not236(n/27)) how[num++]=odd;
        const int nd5=n/5; for (int a=5; a<=nd5; a+=2) if ((n%a)==0) {
          if (not1236(a)&not1236(n/a)) { how[num++]=odd; break; }
        }
        const int j=n-4; if ((j%5==0)&not1236d3(j/5)) how[num++]=_5p4;
        for (int a=5; a<=nd4; a+=2) { if ((a%3)==0) a+=2;
          const int ad2=a/2, nmad2=n-ad2; if (((nmad2%a)==0)&not1236odd(nmad2/a)) { how[num++]=kpf; break; }
        }
        int c=4; while (c<17) { bool found=F;
          const int z=n-c; for (int a=c+1; a<=nd4; a+=2) { if ((a%3)==0) a+=2;
            if (((z%a)==0)&not12356(z/a)) { how[num++]=kpc; found=T; break; }
          }
          if (found) break; c+=(c==4)?4:2;
        }
      }
    } else how[num++]=odd; // not multiple of 3
  } else {
    if ((n>15)&(n<49)&&(random(3)==0)) {
      if (n<28) how[num++]=cn4; else { bool b=randomBool();
        if (n<34) how[num++]=b?cn4:cn8; else if (n<38) how[num++]=b?cn8:cnA;
        else if (n==38) how[num++]=cnA; else if (n<44) how[num++]=b?cnA:cnC; else how[num++]=cnC;
      }
    }
    if (Seven) { // singly even
      if ((n>17)&(n<31)) how[num++]=cn4; if ((n==10)|(n==14)) how[num++]=even;
      else if (n>19) {
        const int nd5=n/5; for (int a=5; a<=nd5; a+=2) if ((n%a)==0) {
          if (not1236(a)&not1236(n/a)) { how[num++]=even; break; }
        }
        const int l=n-2; if (n>30) { if (((l%8)==0)&not1236(l/8)) how[num++]=_8p2; }
        if ((l%5==0)&not1236(l/5)) how[num++]=kpf; const int j=n-4; if ((j%5==0)&not1236(j/5)) how[num++]=_5p4;
        for (int a=7; a<=nd4; a+=2) { if ((a%3)==0) a+=2;
          const int ad2=a/2, nmad2=n-ad2; if (((nmad2%a)==0)&not1236odd(nmad2/a)) { how[num++]=kpf; break; }
        }
        int c=4; while (c<17) { bool found=F;
          const int z=n-c; for (int a=c+1; a<=nd4; a+=2) { if ((a%3)==0) a+=2;
            if (((z%a)==0)&not12356(z/a)) { how[num++]=kpc; found=T; break; }
          }
          if (found) break; c+=(c==4)?4:2;
        }
      }
    } else { /* doubly even */ if (n==24) { how[num++]=_5p4; how[num++]=kpc; } else how[num++]=even; }
  }
  return (num==0)?noHow:how[random(num)];
} // canDo

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

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 sqKind favor) {
//     ------------
  if (favor==pan) { for (int c=0; c<n; ++c) Q[0][c]=c; }
  else { /* assoc, ultra */ Q[0][0]=0; for (int c=1; c<n; ++c) Q[0][c]=n-c; }
  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


sqKind make(int **, const int, const doHow, const sqKind, const bool msg);

void permuteEndToK(int **Z, const int k) {
//   -------------
  const int v=Z[k][k];
  if (v!=k) for (int r=0; r<=k; ++r) for (int c=0; c<=k; ++c) 
    if (Z[r][c]==v) Z[r][c]=k; else if (Z[r][c]==k) Z[r][c]=v;
} // permuteEndToK

void permuteToNBD(int **X, const int n) { // natural order back diagonal
//   ------------
  bool nbd=T; for (int C=0; C<n; ++C) if (X[C][C]!=C) { nbd=F; break; }
  if (!nbd) {
    int *Z=Msquare[0]; for (int C=0; C<n; ++C) Z[X[C][C]]=C;
    for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) X[r][c]=Z[X[r][c]];
  }
} // permuteToNBD

sqKind Make(int **X, const int n, sqKind kind, const bool nbd, const bool msg) {
//     ----
  if ((n&1)&((n%3)!=0)&(kind>=assoc)) {
    if (msg) printf("make %d\n",n); kind=getRowShift2(X,n,kind);
    if (nbd&(kind!=ultra)&(kind!=none)) permuteToNBD(X,n);
  } else {
    kind=make(X,n,canDo(n),kind,msg); if ((kind!=none)&nbd) permuteToNBD(X,n);
  }
  return kind;
} // Make

void incrAndShift(int **Q, const int n, const int u, const int a) {
//   ------------
  const int N=n+a, M=N-1, t=u+a-1;
  // Increment numbers>=u.
  for (int r=0; r<n; ++r) for (int c=0; c<n; ++c) if (Q[r][c]>=u) Q[r][c]+=a;
  // Move right cols, bottom rows to make extra cols, rows in middle
  for (int r=0; r<n; ++r) for (int c=M; c>t; --c) Q[r][c]=Q[r][c-a];
  for (int r=M; r>t; --r) for (int c=0; c<N; ++c) Q[r][c]=Q[r-a][c];
} // incrAndShift

sqKind putSquare(int **Q, const int u, const int Tn, sqKind kind) {
//     ---------
  int **Z=NULL;
  if (allocateIntSquare(&Z,Tn)) {
    if (kind==none) kind=(sqKind)(assoc+random(3)); kind=make(Z,Tn,Tn&1?odd:even,kind,F);
    if (kind!=none) { for (int r=0, s=u; r<Tn; ++r,++s) for (int c=0,d=u; c<Tn; ++c,++d) Q[s][d]=Z[r][c]+u; }
    freeIntSquare(&Z,Tn); 
  } else kind=none;
  return kind;
} // putSquare

void moveLeftUp(int **Q, const int N, const int u, const int t, const int nT) {
//   ----------
  const int ut=u+t; int *Z=Msquare[0], i=0;
  for (int r=0; r<N; ++r) {
    i=0; for (int c=u-nT; c<u; ++c) Z[i++]=Q[r][c];
    for (int c=u; c<ut; ++c) Q[r][c-nT]=Q[r][c]; i=0; for (int c=ut-nT; c<ut; ++c) Q[r][c]=Z[i++];
  }
  void **P=NULL; allocatePtrs(&P,nT); i=0; for (int r=u-nT; r<u; ++r) P[i++]=Q[r];
  for (int r=u; r<ut; ++r) Q[r-nT]=Q[r]; i=0; i=0; for (int r=ut-nT; r<ut; ++r) Q[r]=(int*)P[i++]; freePtrs(&P,nT);
} // moveLeftUp

void compose(int **Q, int **A, int **B, const int nA, const int nB) {
//   -------
  for (int i=0; i<nB; ++i) { const int inA=i*nA; for (int r=0; r<nA; ++r) { const int inAr=inA+r;
    for (int j=0; j<nB; ++j) { const int jnA=j*nA, x=nA*B[i][j]; for (int c=0; c<nA; ++c) Q[inAr][jnA+c]=A[r][c]+x;
  }}}
} // compose

sqKind makeComposite(int **Q, const int nA, const int nB, const bool nbdA, const bool nbdB, const bool msg) {
//     -------------
  int **A=NULL, **B=NULL; sqKind kind=none; const int BisA=(nA==nB)&(nbdA==nbdB);
  if (allocateIntSquare(&A,nA)) {
    if (BisA||allocateIntSquare(&B,nB)) {
      if ((kind=Make(A,nA,ultra,nbdA,msg))!=none) {
        if (BisA) { B=A; } else { kind=Make(B,nB,nB==4?assoc:ultra,nbdB,msg); }
        if (kind!=none) { compose(Q,A, B,nA,nB); kind=magic; }
      }
      if (B!=A) freeIntSquare(&B,nB);
    }
    freeIntSquare(&A,nA);
  }
  return kind;
} // makeComposite
//============================================ special SODLS ================================================

// Squares 8, 9 included for variety.
char ultra8[][8]={ // backtrack
  {0,2,7,5,4,6,3,1},{3,1,4,6,7,5,0,2},{4,6,3,1,0,2,7,5},{7,5,0,2,3,1,4,6},
  {1,3,6,4,5,7,2,0},{2,0,5,7,6,4,1,3},{5,7,2,0,1,3,6,4},{6,4,1,3,2,0,5,7}
};

char ultra9[][9]={ // backtrack
  {0,6,7,4,1,2,8,5,3},{7,1,4,2,5,8,3,6,0},{4,5,2,8,6,3,0,1,7},
  {2,8,6,3,0,1,7,4,5},{6,0,3,1,4,7,5,8,2},{3,4,1,7,8,5,2,0,6},
  {1,7,8,5,2,0,6,3,4},{8,2,5,0,3,6,4,7,1},{5,3,0,6,7,4,1,2,8}
};

// Square 10 not otherwise supported.
char magic10[][10]={ // backtrack
  {0,5,9,1,6,8,4,3,2,7},{2,1,4,7,0,9,5,6,3,8},{1,8,2,9,3,4,7,5,6,0},{4,6,8,3,7,2,0,9,1,5},
  {3,7,0,5,4,6,1,8,9,2},{6,0,7,8,9,5,3,2,4,1},{5,9,3,2,8,1,6,0,7,4},{9,2,1,4,5,3,8,7,0,6},
  {7,4,5,6,2,0,9,1,8,3},{8,3,6,0,1,7,2,4,5,9}
};

// Squares 11, 12 and 13 made by backtracking, included for speed.
char assoc11[][11]={
  {0,2,3,4,10,7,9,8,6,5,1},{3,1,6,7,9,8,4,5,10,0,2},{1,8,2,9,7,4,5,10,3,6,0},
  {5,0,8,3,2,10,7,4,9,1,6},{7,6,5,10,4,1,2,9,0,3,8},{6,3,9,8,0,5,10,2,1,7,4},
  {2,7,10,1,8,9,6,0,5,4,3},{4,9,1,6,3,0,8,7,2,10,5},{10,4,7,0,5,6,3,1,8,2,9},
  {8,10,0,5,6,2,1,3,4,9,7},{9,5,4,2,1,3,0,6,7,8,10}
};

char assoc12[][12]={
  {0,2,3,4,9,8,7,6,11,10,5,1},{3,1,7,6,8,9,5,10,4,11,0,2},{1,9,2,5,7,10,8,11,0,3,6,4},
  {5,4,10,3,11,0,9,8,2,6,1,7},{6,7,11,10,4,1,0,5,9,2,8,3},{11,0,6,8,2,5,4,1,10,7,3,9},
  {2,8,4,1,10,7,6,9,3,5,11,0},{8,3,9,2,6,11,10,7,1,0,4,5},{4,10,5,9,3,2,11,0,8,1,7,6},
  {7,5,8,11,0,3,1,4,6,9,2,10},{9,11,0,7,1,6,2,3,5,4,10,8},{10,6,1,0,5,4,3,2,7,8,9,11}
};

char assoc13[][13]={
  {0,2,3,4,5,6,7,11,10,12,8,9,1},{3,1,8,7,6,12,10,5,9,11,2,0,4},
  {1,10,2,12,9,11,4,8,0,7,3,5,6},{5,0,6,3,11,10,12,9,4,2,1,7,8},
  {2,11,12,8,4,9,1,10,5,3,7,6,0},{7,8,0,6,1,5,9,4,12,10,11,3,2},
  {9,4,7,11,10,0,6,12,2,1,5,8,3},{10,9,1,2,0,8,3,7,11,6,12,4,5},
  {12,6,5,9,7,2,11,3,8,4,0,1,10},{4,5,11,10,8,3,0,2,1,9,6,12,7},
  {6,7,9,5,12,4,8,1,3,0,10,2,11},{8,12,10,1,3,7,2,0,6,5,4,11,9},
  {11,3,4,0,2,1,5,6,7,8,9,10,12}
};

// Two magic order 14, with center 4x4; not otherwise supported.
char magic14a[][14]={ // Adapted from Inder Taneja based on Bennett, Du and Zhang.
  {0,12,6,10,9,13,11,2,1,5,8,7,4,3}, {8,1,10,7,6,9,12,4,0,3,11,13,2,5},
  {4,5,2,9,7,11,13,0,3,1,12,10,8,6}, {6,10,8,3,5,4,0,13,2,7,1,12,9,11},
  {5,13,3,0,4,12,2,10,9,11,7,1,6,8}, {10,0,12,2,1,5,7,8,6,13,3,9,11,4},
  {1,2,9,11,0,8,6,5,7,12,4,3,13,10}, {3,11,13,12,2,6,8,7,5,4,9,0,10,1},
  {2,3,11,13,12,7,5,6,8,10,0,4,1,9}, {11,8,5,1,13,3,10,12,4,9,6,2,0,7},
  {13,6,7,4,11,1,9,3,12,0,10,8,5,2}, {9,7,0,5,8,10,4,1,13,6,2,11,3,12},
  {7,9,4,6,3,2,1,11,10,8,13,5,12,0}, {12,4,1,8,10,0,3,9,11,2,5,6,7,13}
};

char magic14b[][14]={ // Made by backtracking.
  {0,12,9,6,5,1,13,2,3,8,7,4,11,10}, {5,1,0,2,6,11,12,4,13,7,3,10,9,8},
  {6,5,2,11,12,3,9,10,4,1,0,13,8,7}, {11,6,5,3,13,10,2,9,1,4,12,8,7,0},
  {9,13,6,5,4,12,0,1,10,11,8,7,2,3}, {3,2,11,1,10,5,8,6,7,12,13,0,4,9},
  {1,4,10,9,3,7,6,8,5,0,11,12,13,2}, {4,10,13,12,0,8,5,7,6,2,9,1,3,11},
  {10,0,3,4,11,6,7,5,8,13,2,9,1,12}, {12,11,7,8,1,13,10,3,2,9,5,6,0,4},
  {13,7,8,0,2,9,4,11,12,3,10,5,6,1}, {7,8,4,13,9,2,3,12,0,10,1,11,5,6},
  {8,3,1,10,7,0,11,13,9,6,4,2,12,5}, {2,9,12,7,8,4,1,0,11,5,6,3,10,13}
};

// Square 15 included for speed.
char assoc15[][15]={ // backtrack
  {0,2,3,4,5,6,7,8,9,10,11,12,13,14,1},{3,1,4,5,8,7,9,12,10,11,14,13,6,0,2},
  {1,7,2,10,6,9,8,13,11,4,12,14,3,5,0},{5,10,6,3,11,12,1,9,0,14,13,2,7,8,4},
  {2,11,12,7,4,14,0,3,1,13,5,6,9,10,8},{4,12,0,13,10,5,11,14,2,6,7,8,1,9,3},
  {8,13,14,9,12,2,6,10,4,1,0,3,5,11,7},{9,8,10,14,13,11,2,7,12,3,1,0,4,6,5},
  {7,3,9,11,14,13,10,4,8,12,2,5,0,1,6},{11,5,13,6,7,8,12,0,3,9,4,1,14,2,10},
  {6,4,5,8,9,1,13,11,14,0,10,7,2,3,12},{10,6,7,12,1,0,14,5,13,2,3,11,8,4,9},
  {14,9,11,0,2,10,3,1,6,5,8,4,12,7,13},{12,14,8,1,0,3,4,2,5,7,6,9,10,13,11},
  {13,0,1,2,3,4,5,6,7,8,9,10,11,12,14}
};

// Associative 16 included for variety.
char assoc16[][16]={ // backtrack
  {0,2,3,4,5,6,7,8,9,10,11,12,13,14,15,1},{3,1,4,5,6,7,11,14,10,12,13,15,8,9,0,2},
  {1,8,2,10,14,12,4,13,15,5,7,6,9,3,11,0},{5,11,14,3,1,9,13,12,4,0,15,7,2,8,10,6},
  {2,13,8,11,4,15,12,3,14,6,1,5,0,10,9,7},{7,12,9,1,11,5,2,10,0,14,4,13,6,15,8,3},
  {4,14,13,8,15,10,6,2,11,7,12,3,1,0,5,9},{10,9,11,0,13,1,14,7,6,15,2,8,3,5,12,4},
  {11,3,10,12,7,13,0,9,8,1,14,2,15,4,6,5},{6,10,15,14,12,3,8,4,13,9,5,0,7,2,1,11},
  {12,7,0,9,2,11,1,15,5,13,10,4,14,6,3,8},{8,6,5,15,10,14,9,1,12,3,0,11,4,7,2,13},
  {9,5,7,13,8,0,15,11,3,2,6,14,12,1,4,10},{15,4,12,6,9,8,10,0,2,11,3,1,5,13,7,14},
  {13,15,6,7,0,2,3,5,1,4,8,9,10,11,14,12},{14,0,1,2,3,4,5,6,7,8,9,10,11,12,13,15}
};

// Diagonal with center 4.
char dc4s18[][18]={ // backtrack
  {0,11,16,12,5,8,7,1,17,2,3,10,9,4,15,13,6,14},{7,1,17,15,11,6,8,2,16,3,4,9,5,14,12,0,13,10},
  {8,7,2,16,14,17,0,3,15,4,5,6,13,11,1,12,10,9},{1,8,7,3,15,13,16,4,14,5,6,12,17,2,11,10,9,0},
  {15,2,8,7,4,14,12,5,13,6,0,16,3,17,10,9,1,11},{11,14,3,8,7,5,13,6,12,0,1,4,16,10,9,2,17,15},
  {12,17,13,4,8,7,6,0,11,1,2,15,10,9,3,16,14,5},{3,4,5,6,0,1,2,7,10,8,9,17,11,12,13,14,15,16},
  {16,15,14,13,12,11,17,9,8,10,7,2,1,0,6,5,4,3},{17,16,15,14,13,12,11,10,7,9,8,1,0,6,5,4,3,2},
  {4,5,6,0,1,2,3,8,9,7,10,13,14,15,16,17,11,12},{13,6,12,2,9,10,1,15,5,14,16,11,7,8,17,3,0,4},
  {5,13,1,9,10,0,14,16,4,15,17,3,12,7,8,11,2,6},{14,0,9,10,6,15,4,17,3,16,11,5,2,13,7,8,12,1},
  {6,9,10,5,16,3,15,11,2,17,12,0,4,1,14,7,8,13},{9,10,4,17,2,16,5,12,1,11,13,14,6,3,0,15,7,8},
  {10,3,11,1,17,4,9,13,0,12,14,8,15,5,2,6,16,7},{2,12,0,11,3,9,10,14,6,13,15,7,8,16,4,1,5,17}
},

     dc4s22[][22]={ // backtrack
  {0,16,14,8,15,6,17,10,9,1,21,2,3,12,11,5,20,19,4,13,7,18},{9,1,15,13,0,14,7,16,10,2,20,3,4,11,6,19,18,5,21,8,17,12},
  {10,9,2,14,21,1,13,8,15,3,19,4,5,7,18,17,6,20,0,16,12,11},{14,10,9,3,13,20,2,21,0,4,18,5,6,17,16,7,19,1,15,12,11,8},
  {1,13,10,9,4,21,19,3,20,5,17,6,7,15,8,18,2,14,12,11,0,16},{19,2,21,10,9,5,20,18,4,6,16,7,8,0,17,3,13,12,11,1,15,14},
  {5,18,3,20,10,9,6,19,17,7,15,8,0,16,4,21,12,11,2,14,13,1},{16,6,17,4,19,10,9,7,18,8,14,0,1,5,20,12,11,3,13,21,2,15},
  {17,15,7,16,5,18,10,9,8,0,13,1,2,19,12,11,4,21,20,3,14,6},{3,4,5,6,7,8,0,1,2,9,11,12,10,21,13,14,15,16,17,18,19,20},
  {20,19,18,17,16,15,14,13,21,12,10,9,11,2,1,0,8,7,6,5,4,3},{21,20,19,18,17,16,15,14,13,10,12,11,9,1,0,8,7,6,5,4,3,2},
  {4,5,6,7,8,0,1,2,3,11,9,10,12,14,15,16,17,18,19,20,21,13},{18,7,20,21,2,4,11,12,1,19,5,14,17,13,9,10,3,8,16,15,6,0},
  {6,21,13,1,3,11,12,0,19,20,4,15,18,8,14,9,10,2,7,17,16,5},{13,14,0,2,11,12,8,20,5,21,3,16,19,4,7,15,9,10,1,6,18,17},
  {15,8,1,11,12,7,21,4,14,13,2,17,20,18,3,6,16,9,10,0,5,19},{7,0,11,12,6,13,3,15,16,14,1,18,21,20,19,2,5,17,9,10,8,4},
  {8,11,12,5,14,2,16,17,6,15,0,19,13,3,21,20,1,4,18,9,10,7},{11,12,4,15,1,17,18,5,7,16,8,20,14,6,2,13,21,0,3,19,9,10},
  {12,3,16,0,18,19,4,6,11,17,7,21,15,10,5,1,14,13,8,2,20,9},{2,17,8,19,20,3,5,11,12,18,6,13,16,9,10,4,0,15,14,7,1,21}
},

     dc4s26[][26]={ // backtrack
  {0,17,20,16,19,10,21,5,4,12,11,1,25,2,3,14,13,24,6,7,9,23,15,18,8,22},
  {11,1,16,19,15,18,0,20,6,5,12,2,24,3,4,13,23,7,8,10,22,25,17,9,21,14},
  {12,11,2,15,18,25,17,1,19,7,6,3,23,4,5,22,8,9,0,21,24,16,10,20,14,13},
  {7,12,11,3,25,17,24,16,2,18,8,4,22,5,6,9,10,1,20,23,15,0,19,14,13,21},
  {9,8,12,11,4,24,16,23,15,3,17,5,21,6,7,0,2,19,22,25,1,18,14,13,20,10},
  {16,10,9,12,11,5,23,15,22,25,4,6,20,7,8,3,18,21,24,2,17,14,13,19,0,1},
  {5,15,0,10,12,11,6,22,25,21,24,7,19,8,9,17,20,23,3,16,14,13,18,1,2,4},
  {23,6,25,1,0,12,11,7,21,24,20,8,18,9,10,19,22,4,15,14,13,17,2,3,5,16},
  {19,22,7,24,2,1,12,11,8,20,23,9,17,10,0,21,5,25,14,13,16,3,4,6,15,18},
  {22,18,21,8,23,3,2,12,11,9,19,10,16,0,1,6,24,14,13,15,4,5,7,25,17,20},
  {18,21,17,20,9,22,4,3,12,11,10,0,15,1,2,23,14,13,25,5,6,8,24,16,19,7},
  {3,4,5,6,7,8,9,10,0,1,2,11,13,12,14,25,15,16,17,18,19,20,21,22,23,24},
  {24,23,22,21,20,19,18,17,16,15,25,14,12,13,11,2,1,0,10,9,8,7,6,5,4,3},
  {25,24,23,22,21,20,19,18,17,16,15,13,11,14,12,1,0,10,9,8,7,6,5,4,3,2},
  {4,5,6,7,8,9,10,0,1,2,3,12,14,11,13,24,25,15,16,17,18,19,20,21,22,23},
  {17,7,18,4,22,16,3,6,13,14,1,23,5,20,21,15,11,12,19,24,0,9,25,2,10,8},
  {6,19,3,23,17,2,5,13,14,0,18,24,4,21,22,7,16,11,12,20,25,10,8,15,1,9},
  {20,2,24,18,1,4,13,14,10,19,5,25,3,22,23,8,6,17,11,12,21,15,9,7,16,0},
  {1,25,19,0,3,13,14,9,20,4,21,15,2,23,24,10,7,5,18,11,12,22,16,8,6,17},
  {15,20,10,2,13,14,8,21,3,22,0,16,1,24,25,18,9,6,4,19,11,12,23,17,7,5},
  {21,9,1,13,14,7,22,2,23,10,16,17,0,25,15,4,19,8,5,3,20,11,12,24,18,6},
  {8,0,13,14,6,23,1,24,9,17,22,18,10,15,16,5,3,20,7,4,2,21,11,12,25,19},
  {10,13,14,5,24,0,25,8,18,23,7,19,9,16,17,20,4,2,21,6,3,1,22,11,12,15},
  {13,14,4,25,10,15,7,19,24,6,9,20,8,17,18,16,21,3,1,22,5,2,0,23,11,12},
  {14,3,15,9,16,6,20,25,5,8,13,21,7,18,19,12,17,22,2,0,23,4,1,10,24,11},
  {2,16,8,17,5,21,15,4,7,13,14,22,6,19,20,11,12,18,23,1,10,24,3,0,9,25}
},

     dc4s30[][30]={ // backtrack
  {0,22,24,11,23,27,28,8,6,5,4,14,13,1,29,2,3,16,15,20,9,12,25,19,18,17,21,10,7,26},
  {13,1,21,23,12,22,26,27,9,7,6,5,14,2,28,3,4,15,19,10,0,24,18,17,29,20,11,8,25,16},
  {14,13,2,20,22,0,21,25,26,10,8,7,6,3,27,4,5,18,11,1,23,17,29,28,19,12,9,24,16,15},
  {7,14,13,3,19,21,1,20,24,25,11,9,8,4,26,5,6,12,2,22,29,28,27,18,0,10,23,16,15,17},
  {9,8,14,13,4,18,20,2,19,23,24,12,10,5,25,6,7,3,21,28,27,26,17,1,11,22,16,15,29,0},
  {11,10,9,14,13,5,17,19,3,18,22,23,0,6,24,7,8,20,27,26,25,29,2,12,21,16,15,28,1,4},
  {1,12,11,10,14,13,6,29,18,4,17,21,22,7,23,8,9,26,25,24,28,3,0,20,16,15,27,2,5,19},
  {21,2,0,12,11,14,13,7,28,17,5,29,20,8,22,9,10,24,23,27,4,1,19,16,15,26,3,6,18,25},
  {19,20,3,1,0,12,14,13,8,27,29,6,28,9,21,10,11,22,26,5,2,18,16,15,25,4,7,17,24,23},
  {27,18,19,4,2,1,0,14,13,9,26,28,7,10,20,11,12,25,6,3,17,16,15,24,5,8,29,23,22,21},
  {8,26,17,18,5,3,2,1,14,13,10,25,27,11,19,12,0,7,4,29,16,15,23,6,9,28,22,21,20,24},
  {26,9,25,29,17,6,4,3,2,14,13,11,24,12,18,0,1,5,28,16,15,22,7,10,27,21,20,19,23,8},
  {23,25,10,24,28,29,7,5,4,3,14,13,12,0,17,1,2,27,16,15,21,8,11,26,20,19,18,22,9,6},
  {3,4,5,6,7,8,9,10,11,12,0,1,2,13,15,14,16,29,17,18,19,20,21,22,23,24,25,26,27,28},
  {28,27,26,25,24,23,22,21,20,19,18,17,29,16,14,15,13,2,1,0,12,11,10,9,8,7,6,5,4,3},
  {29,28,27,26,25,24,23,22,21,20,19,18,17,15,13,16,14,1,0,12,11,10,9,8,7,6,5,4,3,2},
  {4,5,6,7,8,9,10,11,12,0,1,2,3,14,16,13,15,28,29,17,18,19,20,21,22,23,24,25,26,27},
  {22,19,23,28,10,2,11,6,29,21,15,16,1,27,5,20,25,17,13,14,24,9,3,4,26,18,8,12,0,7},
  {20,24,29,9,1,10,5,17,22,15,16,0,23,28,4,21,26,6,18,13,14,25,8,2,3,27,19,7,11,12},
  {25,17,8,0,9,4,18,23,15,16,12,24,21,29,3,22,27,11,5,19,13,14,26,7,1,2,28,20,6,10},
  {18,7,12,8,3,19,24,15,16,11,25,22,26,17,2,23,28,9,10,4,20,13,14,27,6,0,1,29,21,5},
  {6,11,7,2,20,25,15,16,10,26,23,27,19,18,1,24,29,4,8,9,3,21,13,14,28,5,12,0,17,22},
  {10,6,1,21,26,15,16,9,27,24,28,20,5,19,0,25,17,23,3,7,8,2,22,13,14,29,4,11,12,18},
  {5,0,22,27,15,16,8,28,25,29,21,4,9,20,12,26,18,19,24,2,6,7,1,23,13,14,17,3,10,11},
  {12,23,28,15,16,7,29,26,17,22,3,8,4,21,11,27,19,10,20,25,1,5,6,0,24,13,14,18,2,9},
  {24,29,15,16,6,17,27,18,23,2,7,3,11,22,10,28,20,8,9,21,26,0,4,5,12,25,13,14,19,1},
  {17,15,16,5,18,28,19,24,1,6,2,10,25,23,9,29,21,0,7,8,22,27,12,3,4,11,26,13,14,20},
  {15,16,4,19,29,20,25,0,5,1,9,26,18,24,8,17,22,21,12,6,7,23,28,11,2,3,10,27,13,14},
  {16,3,20,17,21,26,12,4,0,8,27,19,15,25,7,18,23,14,22,11,5,6,24,29,10,1,2,9,28,13},
  {2,21,18,22,27,11,3,12,7,28,20,15,16,26,6,19,24,13,14,23,10,4,5,25,17,9,0,1,8,29}
};

sqKind getUltra8(int **Q) { for (int r=0; r<8; ++r) for (int c=0; c<8; ++c) Q[r][c]=ultra8[r][c]; return ultra; }
//     ---------

void permute3L(int *R, const int n) {
//   ---------
   for (int c=0; c<n; c+=3) { const int t=R[c]; R[c]=R[c+1]; R[c+1]=R[c+2]; R[c+2]=t; }
} // permute3L

void permute3R(int *R, const int n) {
//   ---------
   for (int c=0; c<n; c+=3) { const int t=R[c+2]; R[c+2]=R[c+1]; R[c+1]=R[c]; R[c]=t; }
} // permute3R

void shiftLeftCopy(int *fromR, int *toR, const int n, const int x) {
//   -------------
  const int nmx=n-x; for (int c=0; c<n; ++c) { const int k=c<nmx?c+x:c+x-n; toR[c]=fromR[k]; }
} // shiftLeftCopy

void shiftLeft(int *R, const int n, const int x) {
//   ---------
  int *fromR=Msquare[0]; for (int c=0; c<n; ++c) fromR[c]=R[c]; shiftLeftCopy(fromR,R,n,x);
} // shiftLeft

sqKind getPan9(int **Q) { // based on SODLS from Inder Taneja
//     -------
  const int n=9, nm3=n-3; for (int c=0; c<n; ++c) Q[0][c]=c;
  for (int r=1; r<n; ++r) { shiftLeftCopy(Q[r-1],Q[r],n,3); if ((r==3)|(r==6)) permute3R(Q[r],n); }
  return pan;
} // getPan9

sqKind getUltra9(int **Q) { for (int r=0; r<9; ++r) for (int c=0; c<9; ++c) Q[r][c]=ultra9[r][c]; return ultra; }
//     ---------

sqKind getMagic10(int **Q) { for (int r=0; r<10; ++r) for (int c=0; c<10; ++c) Q[r][c]=magic10[r][c]; return magic; }
//     ----------

sqKind getAssoc11(int **Q) { for (int r=0; r<11; ++r) for (int c=0; c<11; ++c) Q[r][c]=assoc11[r][c]; return assoc; }
//     ----------

sqKind getAssoc12(int **Q) { for (int r=0; r<12; ++r) for (int c=0; c<12; ++c) Q[r][c]=assoc12[r][c]; return assoc; }
//     ----------

sqKind getAssoc13(int **Q) {  for (int r=0; r<13; ++r) for (int c=0; c<13; ++c) Q[r][c]=assoc13[r][c]; return assoc; }
//     ----------

sqKind getMagic14(int **Q, char m14[14][14]) {
//     ---------- 
  for (int r=0; r<14; ++r) for (int c=0; c<14; ++c) Q[r][c]=m14[r][c]; return magic;
} // getMagic14

sqKind getAssoc15(int **Q) { for (int r=0; r<15; ++r) for (int c=0; c<15; ++c) Q[r][c]=assoc15[r][c]; return assoc; }
//     ----------

sqKind getAssoc16(int **Q) { for (int r=0; r<16; ++r) for (int c=0; c<16; ++c) Q[r][c]=assoc16[r][c]; return assoc; }
//     ----------

sqKind getPan27(int **Q) { // based on SODLS from Cheng-Xu Xu and Zhun-Wei Lu
//     -------
  const int n=27; for (int c=0; c<n; ++c) Q[0][c]=c;
  for (int r=1; r<n; ++r) {
    int x=((r%3)==0)?18:9; shiftLeftCopy(Q[r-1],Q[r],n,x);
    if ((r%3)==0) { // 3 6 9 12 15 18 21 24
      x=((r%9)==0)?6:3; shiftLeft(Q[r],9,x); shiftLeft(Q[r]+9,9,x); shiftLeft(Q[r]+18,9,x);
      if ((r%9)==0) permute3L(Q[r],n); // 9 18
    }
  }
  return pan;
} // getPan27

sqKind getSpecial(int **Q, const int n, const sqKind favor)  {
//     ----------
  if ((n>=17)&&isPrime(n)) return getRowShift2(Q,n,favor);
  if (n==1) { Q[0][0]=0; return magic; } // do here before n&1 check below
  switch (favor) {
    case assoc:
      if (randomBool()) {
        if (n==16) return getAssoc16(Q); // not for pan, XformA_D does not make a SODLS
        if ((n&1)&((n%3)!=0)) return getRowShift2(Q,n,assoc);
      }
      break;
    case pan: if (n==9) return getPan9(Q); break;
    case ultra:
      switch (n) {
        case 5: case 7: return getRowShift2(Q,n,ultra);
        case 8: return getUltra8(Q); case 9: return getUltra9(Q); default: break;
      }
    default: break;
  }
  switch (n) {
    case 10: return getMagic10(Q); case 11: return getAssoc11(Q); case 12: return getAssoc12(Q);
    case 13: return getAssoc13(Q); case 14: return getMagic14(Q,randomBool()?magic14a:magic14b);
    case 15: return getAssoc15(Q); case 27: return getPan27(Q); default: break;
  }
  return none;
} // getSpecial
//============================================= associative ================================================

const int B[]={ 1, 2, 4, 8, 16, 32, 64, 128, 256 };
void initSquares(int **Q, const int n, int rU[9], int cU[9], int *fdU, int qrU[9]) {
//   -----------
  for (int i=0; i<n; ++i) { for (int j=0; j<n; ++j) Q[i][j]=-1; rU[i]=0; cU[i]=0; qrU[i]=0; }
  for (int r=0; r<n; ++r) { Q[r][r]=r; rU[r]|=B[r]; cU[r]|=B[r]; qrU[r]|=B[r]; } *fdU=(n&1)?B[n/2]:0;
} // initSquares

sqKind ASODLSf(int **Q, const int n, const sqKind favor)  {
//     -------
  sqKind kind=getSpecial(Q,n,favor); if (kind!=none) return kind;
  // Backtrack to make associative SODLS.
  if (n>9) return none; // Slow or no result
  const int odd=n&1, g=n/2, M=n-1, L=M-1, U=g-1, Ud=odd?g:U, Ur=odd?U:(U-1); int rU[9], cU[9], fdU, qrU[9];
  int r=0, c=M, v=1; enum { row, col, fd } next=fd; initSquares(Q,n,rU,cU,&fdU,qrU);
  do {
    const int ov=Q[r][c], oV=M-ov; Q[r][c]=-1;
    if (next==row) {
      const int R=M-r, C=M-c; if (ov>=0) { rU[r]^=B[ov]; cU[c]^=B[ov]; rU[R]^=B[oV]; cU[C]^=B[oV]; }
      if (v>M) {
        if (c==(r+1)) { if (r==0) { r=U, c=M-U; next=fd; } else { next=col; c=r-1; r=M-c-1; } }
        else { next=col; const int t=r; r=c; c=t; --r; } v=Q[r][c]+1;
      } else {
        int V=M-v;
        while (v<=M) { if (!((rU[r]&B[v])||(cU[c]&B[v])||(rU[R]&B[V])||(cU[C]&B[V]))) break; ++v; --V; }
        if (v<=M) {
          Q[r][c]=v; rU[r]|=B[v]; cU[c]|=B[v]; rU[R]|=B[V]; cU[C]|=B[V];
          next=col; const int t=r; r=c; c=t; v=0;
        }
      }
    } else if (next==col) {
      const int R=M-r, C=M-c; int vs=Q[c][r], Vs=M-vs;
      if (ov>=0) {
        rU[r]^=B[ov]; cU[c]^=B[ov]; qrU[vs]^=B[ov]; qrU[ov]^=B[vs]; 
        rU[R]^=B[oV]; cU[C]^=B[oV]; qrU[Vs]^=B[oV]; qrU[oV]^=B[Vs];
      }
      if (v>M) { next=row; const int t=r; r=c; c=t; v=Q[r][c]+1; }
      else {
        int V=M-v;
        while (v<=M) { if (!((rU[r]&B[v])||(cU[c]&B[v])||(rU[R]&B[V])||(cU[C]&B[V])||(qrU[v]&B[vs]))) break; ++v; --V; }
        if (v<=M) {
          Q[r][c]=v; vs=Q[c][r], Vs=M-vs; rU[r]|=B[v]; cU[c]|=B[v]; qrU[vs]|=B[v];
          rU[R]|=B[V]; cU[C]|=B[V]; qrU[Vs]|=B[V]; qrU[v]|=B[vs]; qrU[V]|=B[Vs];
          next=row; const int t=r; r=c; c=t;
          if (c==(L-r)) {
            if (r==Ur) { for (int i=1; i<n; ++i) for (int j=n-i; j<n; ++j) Q[i][j]=M-Q[M-i][M-j]; kind=assoc; break; }
            else { ++r; c=r+1; v=0; } 
          } else { ++c; v=0; }
        }
      }
    } else if (next==fd) {
      const int R=c, C=r;
      if (ov>=0) {
        rU[r]^=B[ov]; cU[c]^=B[ov]; fdU^=B[ov]; qrU[oV]^=B[ov];
        rU[R]^=B[oV]; cU[C]^=B[oV]; fdU^=B[oV]; qrU[ov]^=B[oV];
      }
      if (v>M) { if (r==0) { break; } else { --r; ++c; v=Q[r][c]+1; }
      } else {
        int V=M-v;
        while (v<=M) { if (!((rU[r]&B[v])||(cU[c]&B[v])||(fdU&B[v])||(rU[R]&B[V])||(cU[C]&B[V]))) break; ++v; --V; }
        if ((r==0)&(v>U)) break; // rotation
        if ((v<=M)) {
          Q[r][c]=v; Q[R][C]=V; rU[r]|=B[v]; cU[c]|=B[v]; fdU|=B[v]; qrU[V]|=B[v];
          rU[R]|=B[V]; cU[C]|=B[V]; fdU|=B[V]; qrU[v]|=B[V]; v=0;
          if (r==U) { next=row; r=0; c=1; } else { ++r, --c; }
        }
      }
    }
  } while (T);
  return kind;
} // ASODLSf
//======================================================= pandiagonal ===============================================

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

sqKind XformA_D(int **Q, const int n) { reverseColumnsBottomHalf(Q,n); reverseRowsRightHalf(Q,n); return pan; }
//     --------

sqKind PSODLSf(int**Q, const int n, const sqKind favor) {
//     -------
  if (isPrime(n)) return getRowShift2(Q,n,favor);
  if ((n&1)&((n%3)!=0)&&randomBool()) return getRowShift2(Q,n,favor);
  const sqKind kind=ASODLSf(Q,n,favor);
  if ((favor==pan)&(kind==assoc)&!(n&3)) return XformA_D(Q,n); return kind;
} // PSODLSf
//============================================ method: x 5's+4 ===============================================

void pat5a(int **Q, const int R, const int u, const int v, int c) { c+=3;
//   -----
  const int R1=R+1, R5=R+5;
  for (int r=R; r<R5; ++r) {
    Q[v][c]=Q[r][c]; Q[r][v]=Q[r][c]; Q[r][c]=v; c+=(r==R1)?-4:1; Q[u][c]=Q[r][c]; Q[r][u]=Q[r][c];Q[r][c]=u; 
  }
} // pat5a

void pat5b(int **Q, const int R, const int x, const int w, int c) { c+=3; 
//   -----
  const int R3=R+3, R5=R+5;
  for (int r=R; r<R5; ++r) {
    Q[x][c]=Q[r][c]; Q[r][w]=Q[r][c]; Q[r][c]=w; c+=(r==R3)?4:-1; Q[w][c]=Q[r][c]; Q[r][x]=Q[r][c]; Q[r][c]=x; 
  }
} // pat5b

void pat5c(int **Q, const int R, const int x, const int w, int c) { ++c; 
//   -----
  const int R1=R+1, R5=R+5;
  for (int r=R; r<R5; ++r) {
    Q[x][c]=Q[r][c]; Q[r][w]=Q[r][c]; Q[r][c]=x; c+=(r==R1)?4:-1; Q[w][c]=Q[r][c]; Q[r][x]=Q[r][c]; Q[r][c]=w; 
  }
} // pat5c

sqKind expand5_4(int **Q, const int n) {
//     ---------
  const int N=n+4, nB=n/5, u=((nB&1)?n+5:n)/2, v=u+1, w=v+1, x=w+1, y=x+1; incrAndShift(Q,n,u,4);

  // Change 5x5 squares to add numbers u, v, w, x and the middle rows, cols numbers.
  if (N&1) { // Change squares on diagonals parallel to the forward diagonal.
    for (int R=0, S=N-10, C=5, D=N-5; R<N; R+=5, S+=5) {
      if (R==u) R+=4; if (S==N) S=0; else if (S==u) S+=4;
      pat5a(Q,R,u,v,C); pat5a(Q,S,w,x,D);
      C-=(C==0)?5-N:(C==y)? 9:5; D-=(D==y)?9:5;
    }
  } else { // Change squares on the back and forward diagonals.
    for (int R=0; R<N; R+=5) { if (R==u) R+=4; const int C=N-R-5;
      pat5a(Q,R,u,v,R); if (R<u) pat5b(Q,R,x,w,C); else pat5c(Q,R,x,w,C);
    }
  }
  sqKind kind=putSquare(Q,u,4,none);
  if (kind!=none) { if (N&1) moveLeftUp(Q,N,u,2,5); kind=magic; } return kind;
} // expand5_4

sqKind makeFrom5Plus4(int **Q, const int N, const bool msg) {
//     --------------
  const int n=N-4, nB=n/5; if (msg) printf("make %d, %d 5s+4\n",N,nB);
  return makeComposite(Q,5,nB,F,F,msg)==magic?expand5_4(Q,n):none;
} // makeFrom5Plus4
//======================================= method: x 8's+2 ==========================================

void pat8b(int **Q, const int R, const int v, const int u, const int z, const int w, int c) {
//   -----                
  ++c; const int Rp1=R+1, Rp4=R+4, Rp5=R+5, Rp8=R+8;
  for (int r=R; r<Rp4; ++r) {
    Q[u][c]=Q[r][c]; Q[r][u]=Q[r][c]; Q[r][c]=w; c+=4;
    Q[v][c]=Q[r][c]; Q[r][v]=Q[r][c]; Q[r][c]=z; c-=(r==Rp1)?7:2;
  } c-=3;
  for (int r=Rp4; r<Rp8; ++r) {
    Q[v][c]=Q[r][c]; Q[r][v]=Q[r][c]; Q[r][c]=z; c+=4;
    Q[u][c]=Q[r][c]; Q[r][u]=Q[r][c]; Q[r][c]=w; c-=(r==Rp5)?7:2;
  }
} // pat8b

void pat8c(int **Q, const int N, const int R, const int v, const int u, const int w) {
//   -----
  const int Rp4=R+4, Rp8=R+8;
  for (int C=0; C<N; C+=8) { if (C==u) C+=2;
    if (Q[R][C]==w) { int c=C;
      for (int r=R; r<Rp4; ++r) { Q[r][c++]=u; Q[r][c]=v; c+=(r==R)?5:-3; } c+=7;
      for (int r=Rp4; r<Rp8; ++r) { Q[r][c--]=u; Q[r][c]=v; c+=(r==Rp4)?-5:3; }
    }
  }
} // pat8c

sqKind expand8_2(int **Q, const int n) {
//     ---------
  const int N=n+2, nB=n/8, u=((nB&1)?n+8:n)/2, v=u+1, w=v+1, z=w+3; incrAndShift(Q,n,u,2);
  // Change 8 x 8 squares on the back diagonal and corresponding order 8 squares to
  // add u and v numbers, and fill the 2 new rows, cols.
  sqKind kind=none;
  for (int R=0; R<N; R+=8) { if (R==u) R+=2;
    if (Q[R][R]==w) { kind=putSquare(Q,u<R?u:R,10,ultra); if (kind==none) break; }
    else { pat8b(Q,R,v,u,z,w,R); pat8c(Q,N,R,v,u,w); }
  }
  if (kind!=none) { if (nB&1) moveLeftUp(Q,N,u,1,8); kind=magic; } return kind;
} // expand8_2

sqKind makeFrom8Plus2(int **Q, const int N, const bool msg) {
//     --------------
  const int n=N-2, nB=n/8; if (msg) printf("make %d, %d 8s+2\n",N,nB);
  return makeComposite(Q,8,nB,T,T,msg)==magic?expand8_2(Q,n):none;
} // makeFrom8Plus2
//====================================== method: x k's+k/2 =========================================

int getKkpf(const int n) {
//  -------
  const int l=n-2; if ((l%5==0)&not1236(l/5)) return 5;
  const int nd4=n/4; for (int a=5; a<=nd4; a+=2) { if ((a%3)==0) a+=2;
    const int ad2=a/2, m=n-ad2; if (((m%a)==0)&not1236odd(m/a)) return a;
  }
  printf("\aProgram Error: getKkpf fail n = %d\n",n); return 5;
} // getKkpf

bool patka(int **Q, const int R, const int nA, const int f, const int t) {
//   -----
  const int Tn=nA+f, RTn=R+Tn; int **Z=NULL; sqKind kind=none;
  if (allocateIntSquare(&Z,Tn)) {
    kind=Make(Z,Tn,ultra,F,F);
    if (kind!=none) for (int r=R; r<RTn; ++r) for (int c=R; c<RTn; ++c) Q[r][c]=Z[r-R][c-R]+t;
    freeIntSquare(&Z,Tn);
  }
  return kind!=none;
} // patka

void patkb(int **Q, const int R, const int nA, const int f, int u, const int t) {
//   -----
  const int RnA=R+nA, RnAm1=RnA-1, RnAm2=RnAm1-1, nAmf=nA-f, nAm1=nA-1; int w=t+nAmf+1;
  for (int C=RnA-f; C<RnA; ++C) { int c=C;
    for (int r=R; r<RnA; ++r) { Q[u][c]=Q[r][c]; Q[r][u]=Q[r][c]; Q[r][c]=w; c+=(c==RnAm1)?-nAm1:1; }
    w=(C==RnAm2)?t:w+1; ++u;
  }
} // patkb

void patkc(int **Q, const int N, const int R, const int nA, const int f, int u, const int t) {
//   -----
  const int RnA=R+nA;
  for (int D=0; D<N; D+=nA) { if (D==u) D+=f;
    if (Q[R][D]==t) { const int DnA=D+nA, Df=D+f;
      for (int C=D; C<Df; ++C) { int c=C;
        for (int r=R; r<RnA; ++r) { Q[r][c]=u; c+=2; if (c>=DnA) c-=nA; } ++u;
      }
    }
  }
} // patkc

sqKind expandk_f(int **Q, const int n, const int nA, const int f) {
//     ---------
  const int N=n+f, nB=n/nA, u=((nB&1)?n+nA:n)/2, t=u-nA; incrAndShift(Q,n,u,f);
  // Change nA x nA squares on the back diagonal to add f new numbers,
  // and fill the f new rows, cols.
  for (int R=0; R<N; R+=nA) { if (R==u) R+=f;
    if (Q[R][R]==t) { if (!patka(Q,R,nA,f,t)) return none; }
    else { patkb(Q,R,nA,f,u,t); patkc(Q,N,R,nA,f,u,t); }
  }
  return magic;
} // expandk_f

sqKind makeFromKplusF(int **Q, const int N, const bool msg) {
//     --------------
  const int nA=getKkpf(N), f=nA/2, n=N-f, nB=n/nA;
  if (msg) printf("make %d, %d %ds+%d\n",N,nB,nA,f);
  return makeComposite(Q,nA,nB,T,T,msg)==magic?expandk_f(Q,n,nA,f):none;
} // makeFromKplusF
//================================================ method: x k's+1 =============================================

int getKkp1(const int n) {
//  -------
  const int m=n-1, nd4=n/4, maxnum=100; int k[maxnum], num=0;
  for (int a=4; a<=nd4; ++a) if (num<maxnum) { if (a==6) ++a; if (((m%a)==0)&not12356(m/a)) k[num++]=a; }
  if (num>0) return k[random(num)]; printf("\aProgram Error: getKkp1 fail n = %d\n",n); return 4;
} // getKkp1

void patk(int**Q, int **Z, const int R, const int C, const int nA, const int t) {
//   ----
  const int RnA=R+nA, CnA=C+nA, s=Q[R][C];
  for (int r=R; r<RnA; ++r) for (int c=C; c<CnA; ++c) if (Z[r-R][c-C]==nA) Q[r][c]=t; else Q[r][c]=Z[r-R][c-C]+s;
  for (int r=R; r<RnA; ++r) Q[r][t]=Z[r-R][nA]+s; for (int c=C; c<CnA; ++c) Q[t][c]=Z[nA][c-C]+s;
} // patk

sqKind expandk_1(int **Q, const int n, const int nA, const int nB) {
//     ---------
  const int N=n+1, u=((nB&1)?nB+1:nB)/2*nA; incrAndShift(Q,n,u,1);
  // Change nA x nA squares on the back diagonal to add new number, and fill the new row, col.
  const int Tn=nA+1; int **Z=NULL; sqKind kind=none;
  if (allocateIntSquare(&Z,Tn)) { kind=Make(Z,Tn,ultra,F,F);
    if (kind!=none) {
      if (Z[nA][nA]!=nA) permuteEndToK(Z,nA);
      for (int R=0; R<N; R+=nA) { if (R==u) ++R; patk(Q,Z,R,R,nA,u); } Q[u][u]=u; kind=magic;
    }
    freeIntSquare(&Z,Tn);
  }
  return kind;
} // expandk_1

sqKind makeFromKplus1(int **Q, const int N, const bool msg) {
//     --------------
  const int nB=getKkp1(N), n=N-1, nA=n/nB; if (msg) printf("make %d, %d %ds+1\n",N,nB,nA);
  return makeComposite(Q,nA,nB,T,T,msg)==magic?expandk_1(Q,n,nA,nB):none;
} // makeFromKplus1
//======================================== method: x k's+c =========================================

int getKkpc(const int n, int *C) {
//  -------
  int c=4, cs[10], ks[10], num=0; while (c<=16) {
    const int z=n-c, nd4=n/4; for (int a=c+1; a<=nd4; a+=2) { if ((a%3)==0) a+=2;
      if (((z%a)==0)&not12356(z/a)) { ks[num]=a; cs[num++]=c; break; }
    } c+=(c==4)?4:2;
  }
  if (num==0) { printf("\aProgram Error: getKkpc fail n = %d\n",n); return 4; }
  const int i=random(num); *C=cs[i]; return ks[i];
} // getKkpc

sqKind expandk_c(int **Q, const int n, const int nA, const int nB, const int C) {
//     ---------
  const int N=n+C, u=((nB&1)?nB+1:nB)/2*nA, z=u+C; incrAndShift(Q,n,u,C);
  // Change nA x nA squares on C diagonals parallel to fd to add new numbers,
  // and fill the new rows, cols.
  const int Tn=nA+1; int **Z=NULL; sqKind kind=none;
  if (allocateIntSquare(&Z,Tn)) {
    kind=Make(Z,Tn,ultra,F,F);
    if (kind==none) freeIntSquare(&Z,Tn);
    else {
      if (Z[nA][nA]!=nA) permuteEndToK(Z,nA); int D0=N-nA-nA, t=z-1, loop=C;
      while (loop--) {
        for (int R=0, D=D0; R<N; R+=nA) { if (R==u) R+=C; patk(Q,Z,R,D,nA,t); D-=(D==0)?nA-N:(D==z)?nA+C:nA; }
        D0-=(D0==z)?nA+C:nA; --t;
      }
      freeIntSquare(&Z,Tn); kind=putSquare(Q,u,C,none); if (kind!=none) { moveLeftUp(Q,N,u,C/2,nA); kind=magic; }
    }
  }
  return kind;
} // expandk_c

sqKind makeFromKplusC(int **Q, const int N, const bool msg) {
//     --------------
  int C=0; const int nB=getKkpc(N,&C), n=N-C, nA=n/nB;
  if (msg) printf("make %d, %d %ds+%d\n",N,nB,nA,C);
  return (makeComposite(Q,nA,nB,T,T,msg)==magic)?expandk_c(Q,n,nA,nB,C):none;
} // makeFromKplusC
//=============================== method: with center sub-square for orders 16 .. 46 ===============================

//void printSODLS(int **x, const int n, FILE *wfp) {
////   ----------
//  for (int r=0; r<n; ++r) { for (int c=0; c<n; ++c) fprintf(wfp, "%3d", x[r][c]); putchar('\n'); }  putchar('\n'); 
//} // printSODLS

sqKind xFormLS(const int n, const int nc, int **from, int **to) {
//     -------
  const int m=n-1, e=n-nc, s=e/2;
  for (int r=0; r<s; ++r) { 
    for (int c=0; c<s; ++c) to[r][c]=from[r][c]; for (int cf=s, ct=m; cf<e; ++cf,--ct) to[r][ct]=from[r][cf];
    for (int cf=n-nc, ct=s; cf<n; ++cf,++ct) to[r][ct]=from[r][cf];
  }
  for (int rf=s, rt=m; rf<e; ++rf,--rt) {
    for (int c=0; c<s; ++c) to[rt][c]=from[rf][c]; for (int cf=s, ct=m; cf<e; ++cf,--ct) to[rt][ct]=from[rf][cf];
    for (int cf=n-nc, ct=s; cf<n; ++cf,++ct) to[rt][ct]=from[rf][cf];
  }
  for (int rf=n-nc,rt=s; rf<n; ++rf, ++rt) {
    for (int c=0; c<s; ++c) to[rt][c]=from[rf][c]; for (int cf=s, ct=m; cf<e; ++cf,--ct) to[rt][ct]=from[rf][cf];
    for (int cf=n-nc, ct=s; cf<n; ++cf,++ct) to[rt][ct]=from[rf][cf];
  }
  permuteToNBD(to,n); /*printSODLS(to, n, stdout);*/ return magic;
} // xFormLS

sqKind makeFromLS(const int n, const int nc, char *s, int **Q, int **LS) {
//     ----------
  const int j=n-nc;
  for (int i=0; i<j; ++i) 
    for (int r=0, c=i, v=s[i]; r<j; ++r) { bool incr=s[i]<j;
      LS[r][c++]=v; if (c==j) c=0; if (incr) { if (++v==j) v=0; }
    }
  for (int c=j; c<n; ++c) for (int r=0, v=s[c]; r<j; ++r) { LS[r][c]=v++; if (v==j) v=0; }
  for (int r=j, i=n; r<n; ++r, ++i) for (int c=0, v=s[i]; c<j; ++c) { LS[r][c]=v++; if (v==j) v=0; }
  putSquare(LS,j,nc,none); return xFormLS(n,nc,LS,Q);
} // makeFromLS

char
  s16c4[]={0,2,1,5,7,9,3,13,14,15,12,4,6,8,11,10,10,7,6,8},
  s18c4[]={0,14,15,16,17,3,11,13,1,6,5,12,8,7,4,10,2,9,2,13,3,4},
  s20c4[]={0,14,11,6,16,13,12,17,9,3,18,7,10,19,2,1,15,5,8,4,5,7,11,15},
  s22c4[]={0,10,17,16,12,4,11,18,19,15,8,5,13,2,6,1,20,21,3,7,9,14,2,11,14,3},
  s24c4[]={0,8,5,14,18,7,1,12,20,6,16,10,4,21,22,3,17,13,11,23,2,9,19,15,18,4,10,9},
  s26c4[]={0,14,3,22,11,4,2,13,16,19,15,20,7,17,23,9,5,24,21,12,25,18,10,8,6,1,14,20,12,2},
  s28c4[]={0,2,10,23,8,21,24,6,17,7,16,25,22,18,26,5,19,4,20,14,3,12,11,27,9,13,1,15,18,17,12,21},
  s30c4[]={0,23,20,2,6,8,17,14,9,26,27,16,7,1,22,25,13,10,12,5,29,28,11,21,15,3,4,24,19,18,16,6,9,13},
  s32c4[]={0,19,4,23,2,8,14,1,21,18,28,22,17,29,15,27,11,10,25,30,26,31,13,20,6,24,12,16,3,5,9,7,15,24,16,4};

sqKind makeWithCenter4(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 4\n",N);
  switch (N) {
    case 16: return makeFromLS(16,4,s16c4,Q,Tsquare);
    case 18: if (randomBool()) return makeFromLS(18,4,s18c4,Q,Tsquare);
      for (int r=0; r<N; ++r) for (int c=0; c<N; ++c) Q[r][c]=dc4s18[r][c]; break;
    case 20: return makeFromLS(20,4,s20c4,Q,Tsquare);
    case 22: if (randomBool()) return makeFromLS(22,4,s22c4,Q,Tsquare);
      for (int r=0; r<N; ++r) for (int c=0; c<N; ++c) Q[r][c]=dc4s22[r][c]; break;
    case 24: return makeFromLS(24,4,s24c4,Q,Tsquare);
    case 26: if (randomBool()) return makeFromLS(26,4,s26c4,Q,Tsquare);
      for (int r=0; r<N; ++r) for (int c=0; c<N; ++c) Q[r][c]=dc4s26[r][c]; break;break;
    case 28: return makeFromLS(28,4,s28c4,Q,Tsquare);
    case 30: if (randomBool()) return makeFromLS(30,4,s30c4,Q,Tsquare);
      for (int r=0; r<N; ++r) for (int c=0; c<N; ++c) Q[r][c]=dc4s30[r][c]; break;break;
    case 32: return makeFromLS(32,4,s32c4,Q,Tsquare);
    default: return none;
  }
  return magic;
} // makeWithCenter4

char
  s17c5[]={0,12,13,14,15,16,11,10,5,8,6,1,9,3,7,4,2,1,10,6,7,4},
  s19c5[]={0,2,1,5,8,3,13,12,14,15,16,17,18,9,6,7,11,4,10,11,9,8,3,6},
  s21c5[]={0,16,17,18,19,14,1,20,15,5,7,12,2,11,8,3,9,13,6,10,4,8,3,15,5,2},
  s23c5[]={0,18,19,20,21,22,7,1,17,14,8,10,4,15,3,12,6,5,9,13,2,11,16,13,3,14,4,11},
  s25c5[]={0,20,21,22,23,24,17,4,12,19,18,5,8,14,13,1,3,10,7,11,2,9,16,15,6,15,18,2,3,5},
  s27c5[]={0,22,23,24,12,11,1,3,21,7,9,20,6,15,19,25,5,26,10,4,17,14,8,18,16,2,13,12,4,3,1,10},
  s29c5[]={0,24,25,26,27,19,21,1,4,13,15,17,9,28,7,18,23,2,10,3,6,8,20,11,5,16,22,12,14,23,1,2,19,13},
  s31c5[]={0,26,27,5,12,28,3,20,29,7,13,1,19,14,23,4,8,16,22,30,25,6,2,17,15,11,9,24,10,18,21,10,21,19,22,14},
  s33c5[]={0,28,29,30,31,26,23,27,13,7,32,20,24,19,16,25,10,14,22,4,6,8,2,18,15,21,5,17,12,3,1,11,9,1,27,11,3,16},
  s35c5[]={0,30,31,17,32,24,12,20,29,2,9,15,33,18,26,11,1,4,28,21,6,19,16,34,27,3,14,22,5,8,25,10,23,7,13,1,22,20,27,11};

sqKind makeWithCenter5(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 5\n",N);
  switch (N) {
    case 17: return makeFromLS(17,5,s17c5,Q,Tsquare); case 19: return makeFromLS(19,5,s19c5,Q,Tsquare);
    case 21: return makeFromLS(21,5,s21c5,Q,Tsquare); case 23: return makeFromLS(23,5,s23c5,Q,Tsquare);
    case 25: return makeFromLS(25,5,s25c5,Q,Tsquare); case 27: return makeFromLS(27,5,s27c5,Q,Tsquare);
    case 29: return makeFromLS(29,5,s29c5,Q,Tsquare); case 31: return makeFromLS(31,5,s31c5,Q,Tsquare);
    case 33: return makeFromLS(33,5,s33c5,Q,Tsquare); case 35: return makeFromLS(35,5,s35c5,Q,Tsquare);
    default: return none;
  }
} // makeWithCenter5

char
  s23c7[]={0,16,17,18,19,20,21,22,12,8,13,4,14,10,15,9,3,7,5,1,6,11,2,14,5,12,11,7,8,6},
  s25c7[]={0,18,19,20,21,22,23,24,6,4,15,10,9,14,17,3,12,7,8,2,11,13,1,5,16,10,12,4,9,7,2,11},
  s27c7[]={0,20,21,22,23,24,25,18,10,12,4,19,5,26,13,7,14,1,8,15,9,2,17,3,16,11,6,1,9,15,17,5,6,7},
  s29c7[]={0,22,23,24,25,26,27,20,28,16,4,13,11,5,7,3,17,14,1,9,6,19,2,12,21,15,8,18,10,17,18,3,6,11,4,9},
  s31c7[]={0,24,25,26,27,28,17,6,14,29,15,5,2,30,7,13,1,12,10,22,3,9,18,20,23,11,8,19,4,21,16,13,15,2,10,1,4,8},
  s33c7[]={0,26,27,28,29,30,31,32,15,19,24,4,3,2,22,9,21,12,10,25,16,6,8,13,1,23,11,14,7,17,20,18,5,2,13,25,23,9,4,1},
  s35c7[]={0,28,29,30,31,32,23,19,3,33,2,18,8,17,16,5,1,10,27,22,
           6,20,34,24,21,12,14,25,9,4,26,15,13,7,11,19,8,10,6,11,22,5},
  s37c7[]={0,23,30,17,2,31,32,16,11,6,29,18,33,21,34,10,12,28,8,13,3,
           35,36,27,4,26,25,9,15,20,22,19,1,7,24,5,14,15,2,23,16,5,6,18};

sqKind makeWithCenter7(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 7\n",N);
  switch (N) {
    case 23: return makeFromLS(23,7,s23c7,Q,Tsquare); case 25: return makeFromLS(25,7,s25c7,Q,Tsquare);
    case 27: return makeFromLS(27,7,s27c7,Q,Tsquare); case 29: return makeFromLS(29,7,s29c7,Q,Tsquare);
    case 31: return makeFromLS(31,7,s31c7,Q,Tsquare); case 33: return makeFromLS(33,7,s33c7,Q,Tsquare);
    case 35: return makeFromLS(35,7,s35c7,Q,Tsquare); case 37: return makeFromLS(37,7,s37c7,Q,Tsquare);
    default: return none;
  }
} // makeWithCenter7

char
  s28c8[]={0,2,13,10,7,17,20,21,22,23,4,8,5,9,18,24,25,26,27,14,1,19,6,3,15,12,16,11,9,2,8,18,19,6,5,10},
  s30c8[]={0,22,23,24,25,26,27,28,29,2,17,13,1,7,20,19,11,15,5,18,10,4,9,16,12,6,8,3,21,14,8,14,19,10,3,13,18,1},
  s32c8[]={0,10,3,22,20,15,17,9,24,25,26,27,1,19,4,13,21,28,
           29,30,31,12,6,11,14,23,18,7,16,5,8,2,4,20,17,3,7,23,21,18},
  s34c8[]={0,7,5,8,25,23,10,4,26,27,28,20,1,24,22,29,18,16,6,30,
           31,19,13,32,17,33,2,9,14,11,12,3,15,21,22,7,10,20,1,13,12,16},
  s36c8[]={0,28,29,30,25,31,32,5,22,10,16,33,24,23,1,12,20,
           13,34,26,9,35,21,15,27,8,3,17,18,4,2,14,19,7,11,6,2,8,23,13,9,16,22,19};

sqKind makeWithCenter8(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 8\n",N);
  switch (N) {
    case 28: return makeFromLS(28,8,s28c8,Q,Tsquare); case 30: return makeFromLS(30,8,s30c8,Q,Tsquare);
    case 32: return makeFromLS(32,8,s32c8,Q,Tsquare); case 34: return makeFromLS(34,8,s34c8,Q,Tsquare);
    case 36: return makeFromLS(36,8,s36c8,Q,Tsquare); default: return none;
  }
} // makeWithCenter8

char
  s29c9[]={0,20,21,22,23,24,25,26,27,28,13,5,1,19,15,12,3,16,10,4,7,11,17,6,9,8,14,18,2,10,18,2,8,13,16,15,4,11},
  s31c9[]={0,22,23,24,25,26,27,28,18,29,30,9,3,16,1,8,15,19,2,20,14,6,11,17,13,10,7,5,4,21,12,12,14,19,18,5,17,11,4,8},
  s33c9[]={0,24,25,26,27,28,29,17,19,11,13,20,9,30,31,7,5,32,
           8,23,2,16,15,14,12,21,10,3,22,18,4,1,6,23,7,18,12,20,22,1,8,5},
  s35c9[]={0,26,27,28,29,30,18,8,7,14,16,5,9,17,31,32,33,25,
           13,2,34,19,6,11,15,21,10,4,23,22,1,3,20,12,24,15,11,13,19,3,18,2,16,7},
  s37c9[]={0,28,29,30,31,6,32,13,33,3,2,34,22,21,9,35,36,15,
           1,23,7,26,18,4,10,27,17,24,5,19,16,14,12,11,20,25,8,17,13,21,16,27,18,3,7,12},
  s39c9[]={0,30,31,32,33,1,3,24,13,34,22,21,4,12,16,10,6,15,
           35,36,29,14,26,37,38,19,7,5,11,18,8,20,27,9,17,28,23,25,2,6,21,16,15,7,1,18,3,14};

sqKind makeWithCenter9(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 9\n",N);
  switch (N) {
    case 29: return makeFromLS(29,9,s29c9,Q,Tsquare); case 31: return makeFromLS(31,9,s31c9,Q,Tsquare);
    case 33: return makeFromLS(33,9,s33c9,Q,Tsquare); case 35: return makeFromLS(35,9,s35c9,Q,Tsquare);
    case 37: return makeFromLS(37,9,s37c9,Q,Tsquare); case 39: return makeFromLS(39,9,s39c9,Q,Tsquare);
    default: return none;
  }
} // makeWithCenter9

char
  s34c10[]={0,23,24,25,11,16,8,17,4,12,26,27,9,5,15,28,29,30,31,
            32,1,3,13,33,22,6,7,18,19,10,2,21,20,14,17,9,18,8,23,19,4,14,12,13},
  s36c10[]={0,26,27,22,2,28,8,29,30,31,32,33,34,7,24,6,4,18,25,1,
            35,20,19,10,14,21,3,23,15,11,5,13,16,12,17,9,9,12,18,15,4,21,6,5,3,11},
  s38c10[]={0,28,19,2,5,10,21,27,29,30,31,17,32,9,16,33,20,11,26,
            1,18,12,34,7,35,36,37,24,8,14,4,15,22,6,25,13,23,3,21,13,14,11,3,23,18,16,7,9},
  s40c10[]={0,23,30,7,17,31,27,32,33,29,34,5,35,20,25,3,2,36,19,
            6,16,10,21,26,37,9,38,39,4,11,13,1,14,28,24,22,8,12,18,15,5,15,8,2,23,27,10,9,25,28},
  s42c10[]={0,14,4,32,33,34,35,15,36,21,28,12,11,37,7,31,25,5,24,10,38,
            39,40,8,2,41,13,22,20,19,27,3,29,18,6,26,16,17,1,23,30,9,26,5,30,28,11,3,7,14,15,21};

sqKind makeWithCenter10(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 10\n",N);
  switch (N) {
    case 34: return makeFromLS(34,10,s34c10,Q,Tsquare); case 36: return makeFromLS(36,10,s36c10,Q,Tsquare);
    case 38: return makeFromLS(38,10,s38c10,Q,Tsquare); case 40: return makeFromLS(40,10,s40c10,Q,Tsquare);
    case 42: return makeFromLS(42,10,s42c10,Q,Tsquare); default: return none;
  }
} // makeWithCenter10

char
  s35c11[]={0,24,25,26,27,28,29,30,15,20,5,21,4,31,32,33,34,19,23,
            9,18,14,6,2,7,8,1,16,17,13,12,3,22,10,11,1,12,6,23,4,15,9,13,21,18,20},
  s37c11[]={0,26,27,28,29,30,31,32,33,16,34,17,21,1,35,10,12,36,5,
            23,4,6,20,15,14,11,24,2,9,22,19,18,25,3,7,13,8,1,23,17,3,8,19,15,20,5,25,2},
  s39c11[]={0,28,29,7,30,3,17,31,32,23,2,33,34,26,22,35,15,27,36,
            37,8,10,38,1,20,19,13,18,24,11,4,9,14,5,16,21,12,25,6,1,21,12,25,18,7,9,2,23,3,5},
  s41c11[]={0,30,31,32,29,33,23,34,11,35,36,19,10,14,27,20,1,21,
            37,15,8,12,38,25,39,17,40,13,4,6,24,22,18,2,9,16,5,3,28,26,7,23,24,10,11,14,12,19,27,9,29,20},
  s43c11[]={0,32,33,34,24,13,35,19,18,36,21,26,37,14,5,4,12,7,20,
            38,15,39,9,28,6,40,25,41,42,10,22,3,27,1,11,17,31,16,30,29,2,8,23,18,7,16,6,9,3,29,17,26,25,30},
  s45c11[]={0,34,35,36,37,38,39,40,41,30,12,42,27,21,43,32,22,15,
            44,33,5,28,26,7,3,11,17,23,20,18,1,13,31,10,25,19,29,16,14,2,8,9,6,24,4,29,28,31,9,22,24,27,12,1,3,10};

sqKind makeWithCenter11(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 11\n",N);
  switch (N) {
    case 35: return makeFromLS(35,11,s35c11,Q,Tsquare); case 37: return makeFromLS(37,11,s37c11,Q,Tsquare);
    case 39: return makeFromLS(39,11,s39c11,Q,Tsquare); case 41: return makeFromLS(41,11,s41c11,Q,Tsquare);
    case 43: return makeFromLS(43,11,s43c11,Q,Tsquare); case 45: return makeFromLS(45,11,s45c11,Q,Tsquare);
    default: return none;
  }
} // makeWithCenter11

char
  s40c12[]={0,28,29,22,30,8,31,32,33,1,19,16,34,20,18,2,11,35,
            36,37,4,9,23,38,17,39,25,12,7,24,15,5,3,14,10,13,27,26,6,21,10,26,6,25,18,8,11,24,17,2,22,14},
  s42c12[]={0,5,20,16,26,1,30,31,32,28,33,18,29,24,34,25,7,10,
            35,36,2,27,23,9,3,37,38,39,40,41,19,8,4,17,13,15,11,22,12,6,21,14,3,5,14,28,15,8,2,27,29,24,25,20},
  s44c12[]={0,9,8,16,23,32,33,34,17,29,3,35,36,20,37,38,10,2,
            14,39,4,18,40,6,22,15,30,5,27,41,42,43,26,11,7,13,24,31,12,28,21,19,25,1,3,14,24,18,2,5,11,21,1,23,27,12},
  s46c12[]={0,32,16,19,10,17,4,29,34,35,36,37,38,15,39,40,41,21,26,
           5,27,42,1,6,20,12,18,11,43,44,24,45,23,2,14,22,8,7,3,25,31,13,9,33,28,30,15,29,27,19,9,5,10,11,1,23,24,33},
  s48c12[]={0,36,37,38,39,40,41,8,42,34,43,19,44,30,45,10,28,46,11,18,2,23,32,13,
            21,29,12,47,24,20,35,14,9,31,5,22,25,17,6,33,27,1,3,15,16,7,4,26,11,28,3,6,14,30,15,16,21,9,24,20};

sqKind makeWithCenter12(int **Q, const int N, const bool msg) {
//     ---------------
  if (msg) printf("make %d, center 12\n",N);
  switch (N) {
    case 40: return makeFromLS(40,12,s40c12,Q,Tsquare); case 42: return makeFromLS(42,12,s42c12,Q,Tsquare);
    case 44: return makeFromLS(44,12,s44c12,Q,Tsquare); case 46: return makeFromLS(46,12,s46c12,Q,Tsquare);
    case 48: return makeFromLS(48,12,s48c12,Q,Tsquare); default: return none;
  }
} // makeWithCenter12
//======================================================== make =========================================================

bool getFactors(const int n, int *nA, int *nB) {
//   ----------
  const int r=(int)sqrt((float)n); int num=0, *fac=NULL; bool ok=F;
  if (allocateInts(&fac,r)) {
    if (n&1) { for (int a=5; a<=r; a+=2) if ((n%a)==0) {
                 const int b=n/a; if (not1236(a)&not1236(b)) { fac[num++]=a; fac[num++]=b; }}
    } else for (int a=4; a<=r; ++a) if ((n%a)==0) {
             const int b=n/a; if (not1236(a)&not1236(b)) { fac[num++]=a; fac[num++]=b; }}
    if (num>0) { const int a=fac[random(num)]; *nA=a; *nB=n/a; ok=T; } freeInts(&fac,n);
  }
  return ok;
} // getFactors

sqKind compositeKind(const sqKind a, const sqKind b) {
//     -------------
  sqKind kind=magic;
  if (a==b) kind=a; else {
    if ((a==ultra)&((b==assoc)|(b==pan))) kind=b; else if ((b==ultra)&((a==assoc)|(a==pan))) kind=a;
  }
  return kind;
} // compositeKind   

sqKind make(int **Q, const int n, const doHow how, const sqKind favor, const bool msg) {
//     ----
  sqKind kind=none;
  switch (how) {   
    case _5p4: kind=makeFrom5Plus4(Q,n,msg); break; // 5s+4
    case _8p2: kind=makeFrom8Plus2(Q,n,msg); break; // 8s+2
    case  kpf: kind=makeFromKplusF(Q,n,msg); break; // ks+k/2
    case  kp1: kind=makeFromKplus1(Q,n,msg); break; // ks+1
    case  kpc: kind=makeFromKplusC(Q,n,msg); break; // ks+C = 4, 8, 10, 12, 14, 16
    case cn4: kind=makeWithCenter4(Q,n,msg); break; // orders 16 .. 32
    case cn5: kind=makeWithCenter5(Q,n,msg); break; // orders 17 .. 33
    case cn7: kind=makeWithCenter7(Q,n,msg); break; // orders 23 .. 33
    case cn8: kind=makeWithCenter8(Q,n,msg); break; // orders 28 .. 24
    case cn9: kind=makeWithCenter9(Q,n,msg); break; // orders 29 .. 39
    case cnA: kind=makeWithCenter10(Q,n,msg); break; // orders 34 .. 40
    case cnB: kind=makeWithCenter11(Q,n,msg); break; // orders 35 .. 45
    case cnC: kind=makeWithCenter12(Q,n,msg); break; // orders 40 .. 46
    default: kind=(favor==assoc)?ASODLSf(Q,n,assoc):PSODLSf(Q,n,favor);
      if (msg&(kind!=none)) printf("make %d\n",n); break;
  }
  if (kind!=none) return kind; int nA=0, nB=0;
  if (getFactors(n,&nA,&nB)) {
    if (msg) printf("make %d, %d %ds\n",n,nB,nA); int **A=NULL, **B=NULL;
    if (allocateIntSquare(&A, nA)) {
      if ((nA==nB)||allocateIntSquare(&B, nB)) {
        const sqKind kindA=make(A,nA,canDo(nA),favor,msg); sqKind kindB;
        if (nA==nB) { kindB=kindA; B=A; } else kindB=make(B,nB,canDo(nB),kindA,msg);
        if (kindB!=none) { compose(Q,A,B,nA,nB); kind=compositeKind(kindA,kindB); }
        if (B!=A) freeIntSquare(&B,nB);
      }
      freeIntSquare(&A,nA);
    }
  }
  return kind;
} // make

bool makeSquare(int **Q, const int n, kindArray numSq) { 
//   ----------
  doHow how=canDo(n);
  if (how) {
    const int f=randomBool(), a=random(3), p=random(5);
    bool ok=(make(Q,n,how,f?a?assoc:ultra:p?pan:ultra,T)!=none);
    if (!ok) { printf("\aProgram Error: fail n = %d\n",n); ++numSq[none]; } return ok;
  }
  printf("Not supported.\n\n"); ++numSq[none]; return F;
} // makeSquare
//============================================== main =============================================

void outputLocalTime() {
//   ---------------
  time_t startTime=time(NULL);
  struct tm local; localtime_s(&local,&startTime); char dateTime[100];
  size_t slen=strftime(dateTime,100,"%A %Y-%m-%d %X %Z\n\n\0",&local);  printf(dateTime);
} // outputLocalTime

void printElapsedTime(time_t startTime) {
//   ---------------- 
  const int et=(int)difftime(time(NULL),startTime);
  if (et>0) printf("\nelapsed time %d:%02d:%02d\n",et/3600,et%3600/60,et%60);
} // printElapsedTime
 
int main() {
//  ----
  bool goodInput=F, ok=F; int n1=0, n2=0, retry=3; outputLocalTime(); seed_rand(); allocatedSize=0;
  while (retry--) {
    printf("Square order or order range? "); if (getInts(' ',&n1,&n2)) { goodInput=T; break; }
  }
  if (goodInput&&openDir()) { bool another=T, fileOut=F; outMsg=F;
    do {
      if (n1<=outMaxN) fileOut=T; kindArray numSq;
      for (int i=none; i<=ultra; ++i) numSq[i]=0;
      if (allocateStore(n2)) {
        printf("Output, (1 SODLS only, 2 magic square only, or 3 both)? ");
        int out=getInt(); if (out<=0) out=1; else if (out>3) out=3; outputError=F;
        time_t startTime=time(NULL);
        for (int i=n1; i<=n2; ++i) if ((i>6)||validOrder(i,numSq)) {
          int wfd=(i>outMaxN) ? -1 : openOutput(i); int Out=(i>outMaxN) ? 0 : out; ok=(i>outMaxN)|(wfd!=-1);
	        if (ok) { ok=makeSquare(Qsquare,i,numSq);
            if (ok) ok=outputSquares(Msquare,Qsquare,i,wfd,numSq,Out);
            if (wfd!=-1) { _close(wfd); if (!ok) { remove(outputFile); if (outputError) break; } }
	        }
        }
        if (n1!=n2) printf("none %d magic %d assoc %d pan %d ultra %d\n",
          numSq[none],numSq[magic],numSq[assoc],numSq[pan],numSq[ultra]); printElapsedTime(startTime); 
      }
      printf("Continue? n (no) or square order or order range: "); another=getNorOrders(&n1,&n2);
    } while (another);
    if (!fileOut) { if (_chdir("..")==0) _rmdir(outputDir); }
  }
  freeStore(allocatedSize); printf("\nPress a key to close the console");
  while (!_kbhit()) Sleep(250); return ok ? EXIT_SUCCESS : EXIT_FAILURE;
} // main