#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
			/* number of rows and columns in the magic square; larger values will
			yield slow progress .... */
#define N 6
			/* As the sum of the numers 1..N*N in the square
			is  N*N*(1+N*N)/2, and the number of rows and columns N,
			the magic sum must be N*(1+N*N)/2 */
#define MAGSUM ((N*(1+N*N))/2)
			/* recognize some trivial symmetry operations (rotation, mirror of a
			square) to generate magic squares out of known ones */
#define OCT
			/* place the output in a file ...*/
/*#define CREATEFILE*/
			/* generate magic squares such that the first one is chosen randomly out of the set of all possible ones */
#define RANDOMSEED

static count=0,		/* number of magic squares found */
       first=1, last=N*N ;	/* interval of values that are placed in the upper left corner
				of the magic square. 1 and N*N mean, the total amount of possible values
				is scanned. A smaller interval may be chosen to count the number in 'parallel'
				sessions and add the results later on . */
				/* Image the four vertices IND1, IND2, IND3 and IND4 that are mapped onto
				each other by mirroring at the horizontal and vertical middles of the square
				or by rotating the square around its center. Chose the quadruple of
				numbers with smallest number IND4 (to have fastest access to all of these
				vertices, if we have reached don==IN4+1 below). There are 24=4*3*2*1
				permutations of 4 different numbers placed at these 4 vertices. E.g for
				the numbers 1-2-3-4 we have
				1-2-3-4, 1-3-2-4, 4-3-2-1, ...,3-4-1-2, 2-4-1-3   (8 mirrors)
				1-3-4-2, 1-4-3-2, 2-3-4-1, 2-4-3-1, 4-1-2-3, 3-1-2-4, 3-2-1-4,4-2-1-3 (8 mirrors)
				1-2-4-3, 1-4-2-3, 4-1-3-2, 2-1-3-4, 3-2-4-1, 3-4-2-1, 2-3-1-4, 3-4-1-2 (8 mirrors).
				The code, if OCT is defined, means to reject 7/8 of them and only retain the
				non_equivalent 1-2-3-4, 1-3-4-2 and 1-2-4-3 that cannot be created by
				mutual mirroring and rotations. */
#ifdef OCT
#if N%2==0
		/* row N/2-1, column N/2-1 */
#define IND1 ((N+1)*(N/2-1))
		/* row N/2-1, column N/2 */
#define IND2 ((N-1)*(N/2))
		/* row N/2, column N/2-1 */
#define IND3 ((N+1)*N/2-1)
#define IND4 ((N+1)*N/2)
#else
		/* row N/2-1, column N/2 */
#define IND1 ((N+1)*(N/2)-N)
		/* row N/2, column N/2-1 */
#define IND2 ((N+1)*(N/2)-1)
#define IND3 ((N+1)*(N/2)+1)
#define IND4 ((N+1)*(N/2)+N)
#endif
#endif
main(int argc, char * argv[])
{
	int magic[N*N] ,indx;
	void permut(const int *,const) ;
			/* initial permutation of numbers */
#ifdef RANDOMSEED
	void randstart(int *) ;
	randstart(magic) ;
#else
	if( argc == 3 )
	{
		first=atoi(argv[1]) ;
		last=atoi(argv[2]) ;
	}
				/* now place the numbers first, first+1 etc at the front (smaller indices)
				of magic[]. This will let them be scanned first with the kind of permutation
				chosen */
	for( indx=first; indx <= N*N ; indx++)
		magic[indx-first]=indx ;
				/* and let the rest of them fill the tailing indices of the field */
	for( indx=1; indx < first ; indx++)
		magic[N*N+indx-first]=indx ;
#endif
	permut(magic,0) ;
	printf("N=%d, count=%d\n",N,count) ;
}
/******************************
 Run thru all permutations of numbers in 'field', but
leave the first 'done' numbers unchanged.
****************/
static void permut( const int *field, const done)
{
	int newf[N*N] ,
	    ismagic(const int*);
				/* do not run into the subspace of permutations not selected
				on the command line */
	if( field[0]<first || field[0]>last)
		return ;
#ifdef OCT
	if( done == IND2+1 && field[IND2] < field [IND1] )
		return ;
	if( done == IND3+1 && field[IND3] < field [IND2] )
		return ;
	if( done == IND4+1 && field[IND4] < field [IND1] )
		return ;
#else
				/* consider not the permutations mirrored along the 2nd
				diagonal, but multiply the number of magic squares with 2 at the end .
				The decision, which of the two mirrored permutations to take, is
				done by an ordering (equality of 2 numbers in the square cannot occur),
				i.e. the one with the greater number in the upper right triangle is taken.*/
	if(done == N+1)
		if( field[1]> field[N] )
			return ;
#endif
				/* here's the case when the recursion ends and we
				have a full square to be checked with respect to 'magic' */
	if( done == N*N-1)
#ifdef OCT
		count += 8*ismagic(field) ;
#else
		count += 2*ismagic(field) ;
#endif
	else	
	{
		int exch ,tmp ;
		void permut(const int *,const) ;
				/* if we are looking for the last number in the
				current row, this is fixed by demanding it completing the magic number. */
		if( done%N == N-1  )
		{
			int needc(const int *,const),
			    fill ,
			    indexin(const int *, const, const );
			fill=needc(field,done/N) ;
						/* do not allow negative numbers or numbers greater than N*N
						in the magic square. These 2 lines are not really necessary,
						because indexin() would detect this, too, but it hopefully
						speeds up the code, if indexin() is not invoked too often. */
			if( fill < 1 || fill > N*N)
				return ;
			else if( (exch=indexin(field,fill,done)) == -1 )
				return;
			else			/* just the same code as below, but only for the one permutation */
			{
				memcpy((void *)newf,field,N*N*sizeof(int) ) ;
				newf[exch]=newf[done] ; newf[done]=fill ;
				permut(newf,done+1) ;
			}
		}
				/* if we're permuting the last row ... */
		else if ( done >= (N-1)*N)
		{
			int needr(const int *,const),
			    fill , col,
			    indexin(const int *, const, const );
			memcpy((void *)newf,field,N*N*sizeof(int) ) ;
			for( col= done%N ; col < N-1 ; col++)
			{
				fill=needr(newf,col) ;
				if( (exch=indexin(newf,fill,N*(N-1)+col)) == -1 )
					return;
				else			/* just the same code as below, but only for the one permutation */
				{
					newf[exch]= newf[N*(N-1)+col] ; newf[N*(N-1)+col]=fill ;
				}
			}
			permut(newf,N*N-1) ;
		}
					/* column 1 in first but last row must fulfill sum rules on
					the diagonal and serve to fulfill sum rule on 0th column */
		else if ( done == (N-2)*N+1 )
		{
			int needr(const int*, const), needd(const int * ),
			    fill_last, lastr, fill_now, nowr,
			    indexin(const int*, const,const) ;
			fill_now=needd(field) ;
			if( fill_now < 1 || fill_now > N*N || fill_last==field[N*(N-1)] )
				return ;
			else if( (nowr=indexin(field,fill_now,done)) == -1 )
				return;
			memcpy((void *)newf,field,N*N*sizeof(int) ) ;
					/* place these values at their correct site.
					Pay attention on not exchanging newf[done] and newf[N*(N-1)] 2 times */
			newf[nowr]=newf[done] ; newf[done]=fill_now ;
			fill_last=needr(newf,1) ;
			if( fill_last < 1 || fill_last > N*N || fill_last==newf[N*(N-1)] )
				return ;
			else if( (lastr=indexin(newf,fill_last,done+1)) == -1 )
				return;
			newf[lastr]=newf[N*(N-1)+1] ; newf[N*(N-1)+1]=fill_last ;
			permut(newf,done+1) ;
		}
		else if ( done >= (N-2)*N )	/* if we're permuting the first but last row (col 0 or 2,3, and above) */
		{
			int needr(const int *,const),
		   	    fill , lastr ,
			    indexin(const int *, const, const );
			const col=done%N ;
			for( exch = done; exch < N*N ; exch++)
			{
				if( exch%N < col ) 
					continue ;
				memcpy((void *)newf,field,N*N*sizeof(int) ) ;
				tmp= newf[done] ; newf[done]=newf[exch] ; newf[exch]=tmp ;
							/* Look which elements must be reserved for the last row of
							the current column */
	 			fill=needr(newf,col) ;
				if( fill < 1 || fill > N*N )
					continue ;
				if( (lastr=indexin(newf,fill,done+1)) == -1 )
					continue;
				if ( lastr%N < col)
					continue ;
				newf[lastr]=newf[(N-1)*N+col] ; newf[(N-1)*N+col]=fill ;
				permut(newf,done+1) ;
			}
		}
		else
			for( exch = done; exch < N*N ; exch++)
			{
				memcpy((void *)newf,field,N*N*sizeof(int) ) ;
				tmp= newf[done] ; newf[done]=newf[exch] ; newf[exch]=tmp ;
				permut(newf,done+1) ;
			}
	}
}
/**********************************
 return 1, if the numbers in field[] are a magic square,
 else 0
**********************************/
static ismagic(const int *field)
{
	int col, row, rsum ;
			/* check sum of each row */
	for(row=1 ; row < N ; row++)
	{
		for( rsum = col=0 ; col < N ; col++)
			rsum += field[col+row*N] ;
		if ( rsum != MAGSUM )
			return 0 ;
	}
			/* check diagonal 0, 1+N, 2+2N ... */
	for(col=rsum=0 ; col < N ; col++)
		rsum += field[col*(1+N)] ;	/* col+col*N the index */
	if ( rsum != MAGSUM )
		return 0 ;
			/* check diagonal N(N-1), (N-1)**2, ... */
	for(col=rsum=0 ; col < N ; col++)
		rsum += field[col+(N-1-col)*N] ;
	if ( rsum != MAGSUM )
		return 0 ;
			/* check sum of each column */
	for(col=0 ; col < N ; col++)
	{
		for( rsum = row=0 ; row < N ; row++)
			rsum += field[col+row*N] ;
		if ( rsum != MAGSUM )
			return 0 ;
	}
			/* the code below displays the magic squares found */
#ifdef CREATEFILE
	{
		char fname[9] ; FILE *outf ;
		sprintf(fname,"mag%1d.dat",N) ;
		outf=fopen(fname,"a+") ;
		fprintf(outf,"\n") ;
		for(row=0; row < N ; row++)
		{
			for(col=0; col < N ; col++)
				fprintf(outf,"%3d",field[col+row*N]) ;
			fprintf(outf,"\n") ;
		}
		fclose(outf) ;
	}
#endif
	{
		printf("\n") ;
		for(row=0; row < N ; row++)
		{
			for(col=0; col < N ; col++)
				printf("%3d",field[col+row*N]) ;
			printf("\n") ;
		}
	}
			/* return 1, if the square passed the row, column
			and the 2 diagonal tests */
	return 1 ;
}
/**********************************
 return the number that is needed as the last row
 of column 'col' to have the column sum up to the magic number
********************************************/
static needr(const int *field, const col)
{
	int sum, row ;
	for( sum=MAGSUM, row=0 ; row < N-1 ; row++)
		sum -= field[col+row*N] ;
	return sum ;
}
/**********************************
 return the number that is needed as the last column
 of row 'row' to have the row sum up to the magic number
********************************************/
static needc(const int *field, const row)
{
	int sum, col ;
	for( sum=MAGSUM, col=0 ; col < N-1 ; col++)
		sum -= field[col+row*N] ;
	return sum ;
}
/**********************************
 return the number that is needed as the column '1'
 of the last but first row to have the diagonal sum up to the magic number
********************************************/
static needd(const int *field)
{
	int dsum, col ;
	for( dsum=MAGSUM-field[(N-1)*N], col=2 ; col < N ; col++)
		dsum -= field[col+(N-1-col)*N] ;
	return dsum ;
}
/**************************************
 return the index of the occurrence of number 'wanted'
 in the array 'field' up from position 'done'. If the
 number is not found, -1 is returned.
************************************/
static indexin(const int *field, const register wanted, const done)
{
	register posit ;
	for(posit=done ; posit < N*N ; posit++)
		if( field[posit] == wanted )
			return posit ;
	return -1 ;
}
#ifdef RANDOMSEED
/******************************************************
 generate a random starting permutation
***********************************************/
static void randstart(int *magic)
{
	unsigned int now;
	int indx ;
				/* start with the simple sequential ordering */
	for( indx=1; indx <= N*N ; indx++)
		magic[indx-1]=indx ;
				/* generate random seed for rand() */
	now=time(NULL) ;
	srand(now) ;
	for( indx=0 ; indx < N*N ; indx++)
	{
		int randindx=rand() ,
		    tmp ;
		const double matchslope=(N*N-0.5-indx)/RAND_MAX ;
		randindx= matchslope*randindx+indx ;
		tmp=magic[randindx] ;
		magic[randindx]=magic[indx] ;
		magic[indx]=tmp ;
	}
}
#endif
#undef N
#undef MAGSUM