// ====== Program SquareGrid =====================
//  Written: 31.January 2001, Matthias Engelhardt
//     Matthias.Engelhardt(AT)nbg.siemens.de

import java.math.BigInteger;

/**  Compute number of different patterns that arise from permutations on a torus.
*
*  If you have a permutation pi, i.e.an element of the symmetric group Sn,
* you may visualize it by drawing a dot at all (k, pi(k)) in a coordinate system
* or placing a chess pawn at field (k, pi(k)) on a n by n chess board, for all 1<=k<=n.
* Different permutations will lead to different patterns.

*The patterns may similar in the following sense: they may be rotations,
* reflections of torus shifts of other patterns. Torus shift means that you shift
* end around, i.e. if you arrive at n you start with 1 again.

*  Program SquareGrid computes the number of patterns which are different also
* under these operations.
*
*  Counting that number is equivalent to counting classes of sequential binary
* arrays for the sequence 000001 ((n-1) copies 0, one copy 1). That problem is
* discussed in the paper

Anne Penfold Street and Robert Day
Sequential binary arrays II: Further results on the
square grid,
pp. 392-418 of Combinatorial Mathematics IX. Proc. Ninth
Australian Conference
(Brisbane, August 1981).
Ed. E.J.Billington, S.Oates-Williams and A.P.Street.
Lecture Notes Math., 952. Springer-Verlag, 1982.

*  The program uses the formulae that are given in that
paper. The program showed,
* however, that we need two correction terms for n = 1 and n = 2.
*/

public class SquareGrid {

// Two arrays of prepared values
static BigInteger [] tabFactorial;
static BigInteger [] tabEulerFi;
static int    nMax;


public static void main (String [] args) {

int n, d, h = 0, h2, k, i, l, iWk, rmd4;

BigInteger [] nrCoset = new BigInteger [5];
final int idxN = 0,  idxNw2 = 1,  idxNw = 2,
idxNx = 3,  idxNwx = 4;

BigInteger wk1, wk2, wk3, big2, bigH, bigN,
product, numerator, denominator,
subSum, allCoset, Tn;

BigInteger [] divisionResult;  // Array containing result
and remainder
// given by method divideAndRemainder

nMax = 30;
prepareValues ();

big2 = BigInteger.ONE.add (BigInteger.ONE);

for (n = 1; n < nMax; n++) {
//Diag System.out.print ("Diag: starting for " + n + "\n");
// Program follows the cases in the paper (Page 396).

rmd4 = n % 4; k = n / 4;
switch (rmd4) {
case 0:  h = n / 2;   break;
case 1:  h = (n-1)/2;  break;
case 2:  h = n / 2;   break;
case 3:  h = (n-1)/2;  break;
}
bigN = makeBig (n);
bigH = makeBig (h);

// First all cosets without reflection
// Compute nrCoset [idxN],    i.e. subgroup itself.
wk1 = BigInteger.ZERO;
for (d = 1; d <= n; d++) {
if (n % d == 0) {
iWk = n / d;
wk2 = eulerFi (iWk).pow (2);
wk3 = makeBig (iWk).pow (d);
wk2 = wk2.multiply (wk3);

wk3 = factorial (d);
wk2 = wk2.multiply (wk3);

wk1 = wk1.add (wk2);
}
}         // all factors d of n
nrCoset [idxN] = wk1;

// Compute nrCoset [idxNw2],   i.e. 180� rotation.
wk3 = big2.pow (h);
if (rmd4 == 0 || rmd4 == 2) {
wk2 = factorial (h + 1);
wk1 = wk3.multiply (wk2).multiply (bigH);
} else {
wk2 = factorial (h);
wk1 = bigH.multiply (big2).add (BigInteger.ONE).pow (2);
wk1 = wk1.multiply (wk2).multiply (wk3);
}
nrCoset [idxNw2] = wk1;

// Compute nrCoset [idxNw],   i.e. 90� rotation (or 270�).
// For that, compute the product first
product = BigInteger.ONE;
for (l = 1; l <= k; l++) {
wk3 = makeBig (2 * l - 1);
product = product.multiply (wk3);
}

switch (rmd4) {
case 0:
wk3 = makeBig (k).pow (2);
wk2 = big2.pow (k + 3);
wk1 = wk2.multiply (wk3).multiply (product);
break;
case 1:
wk3 = makeBig (n).pow (2);
wk2 = big2.pow (k);
wk1 = wk2.multiply (wk3).multiply (product);
break;
case 2:
wk3 = makeBig (2 * k + 1).pow (2);
wk2 = big2.pow (k + 2);
wk1 = wk2.multiply (wk3).multiply (product);
break;
case 3:
wk1 = BigInteger.ZERO;
break;
}
nrCoset [idxNw] = wk1;

// Now cosets with reflection (x means reflection)

// Compute nrCoset [idxNx],   i.e. only reflection.
// This part is the most complex part (complex not in
mathematical sense)
wk1 = BigInteger.ZERO;
for (d = 1; d <= n; d++) {
if (n % d == 0) {   // Take only factors of n
numerator = factorial (d).multiply (makeBig (n)).multiply
(eulerFi (n/d));
if ((n / d) % 2 == 1) {
subSum = BigInteger.ZERO;
for (l = 0; 2 * l <= d; l ++) {
denominator = makeBig (2 * d).pow (l);
denominator = denominator.multiply (factorial (l));
denominator = denominator.multiply (factorial (d - 2 *
l));

wk3 = numerator.multiply (bigN.pow (l));
divisionResult = wk3.divideAndRemainder (denominator);
if (! divisionResult [1].equals (BigInteger.ZERO)) {
System.out.print ("Error: division results in a
remainder;"
+ " that is unexpected! ! !\n");
System.exit (7);
}
subSum = subSum.add (divisionResult [0]);
}
wk1 = wk1.add (subSum);
} else {
if (rmd4 == 0) {
if (d % 2 == 0) {
h2 = d / 2;
denominator = makeBig (2 * d).pow (h2);
denominator = denominator.multiply (factorial (h2));

wk3 = numerator.multiply (bigN.pow (h2));
divisionResult = wk3.divideAndRemainder (denominator);
if (! divisionResult [1].equals (BigInteger.ZERO)) {
System.out.print ("Error: division results in a
remainder;"
+ " that is unexpected! ! !\n");
System.exit (7);
}
wk1 = wk1.add (divisionResult [0]);
}     // even d
}      // n divisible by 4
}       // even n / d
}        // d factor of n
}         // all factors d of n
nrCoset [idxNx] = wk1;


// No nrCoset [idxNw2x],    i.e. reflection plus 180�
rotation:
//            it is the same as reflection only!

// Compute nrCoset [idxNwx],   i.e. reflection plus 90�
rotation.
wk3 = big2.pow (h);
if (rmd4 == 0 || rmd4 == 2) {
wk2 = factorial (h);
nrCoset [idxNwx] = wk3.multiply (wk2).multiply (bigH);
} else {
nrCoset [idxNwx] = BigInteger.ZERO;
}

// Corrections for small values:
if (n == 1) {
nrCoset [idxNx]  = BigInteger.ONE;
nrCoset [idxNwx] = BigInteger.ONE;
}
if (n == 2) {
nrCoset [idxNx]  = makeBig (4);
nrCoset [idxNwx] = makeBig (4);
}

// Now take the sum; for that, I sum up first the three
parts that occur twice.
allCoset = nrCoset [idxNw].add (nrCoset [idxNx]).add
(nrCoset [idxNwx]);
allCoset = allCoset.add (allCoset);  // double it
allCoset = allCoset.add (nrCoset [idxN]).add (nrCoset
[idxNw2]);

Tn = allCoset.divide (makeBig (8 * n * n));
System.out.print (" " + Tn.toString () + ",\n");
}          // all my n
}          // End of method main

/** Convert variable of simple type int to a BigInteger.
* Class BigInteger has no constructor for that (I think that
would
* be more appropriate). Therefore, I take the "deviation"
via String.
*/
public static BigInteger makeBig (int intIn) {
return new BigInteger (new Integer (intIn).toString ());
}

/** Prepare values for factorials and Euler's Fi function.
*/
static public void prepareValues () {
int i;
int n, d, h, k;
BigInteger wk1, wk2, wk3;

tabFactorial = new BigInteger [nMax + 1];
tabEulerFi = new BigInteger [nMax + 1];

tabFactorial [0] = BigInteger.ONE;
tabFactorial [1] = BigInteger.ONE;
tabEulerFi [1] = BigInteger.ONE;

for (n = 2; n < nMax; n++) {
wk1 = new BigInteger (new Integer (n).toString ());
tabFactorial [n] = tabFactorial [n-1].multiply (wk1);

tabEulerFi [n] = BigInteger.ZERO;
for (k = 1; k < n; k++) {
wk2 = new BigInteger (new Integer (k).toString ());
// Test if k is relativ prime to n, i.e. greatest common
divisor is one.
if (wk2.gcd (wk1).equals (BigInteger.ONE)) {
tabEulerFi [n] = tabEulerFi [n].add (BigInteger.ONE);
}
}
}          // prepare all needed values
}          // End prepareValues

static public BigInteger eulerFi (int i) {
if (i > nMax) {
System.out.print ("eulerFi: argument " + i + " exceeds
generated size\n");
System.exit (9);
}
return (tabEulerFi [i]);
}          // eulerFi

static public BigInteger factorial (int i) {
if (i > nMax) {
System.out.print ("factorial: argument " + i + " exceeds
generated size\n");
System.exit (9);
}
return (tabFactorial [i]);
}          // factorial

}          // End of class SquareGrid