/* -*- MaTX -*-
 *
 *  NAME
 *      magicsq() - Magic Square
 *
 *  SYNOPSIS
 *      M = magicsq(n)
 *         Matrix M;
 *         Integer n;
 *
 *  DESCRIPTION
 *      magicsq(n) computs the n-by-n Magic Square matrix.
 */

Func Matrix magicsq(n)
	Integer n;
{
	Integer i, j, k;
	Matrix A, M;

	A = M = Z(n);
	if (rem(n,2) == 0) {
		// Make a seed matrix of magic square of even-order
		for (i = 1; i <= n/2; i++) {
			for (j = 1; j <= n; j++) {
				if (rem(i,2) == 1) {
					A(i,j) = j;
					A(n-i+1,j) = j;
				} else {
					A(i,j) = n - j + 1;
					A(n-i+1,j) = n - j + 1;
				}
			}
		}

		// 4m + 2 order
		if (rem(n-2,4) == 0) {
			A(n/2,:) = fliplr(A(n/2,:));
			A(n/2,Index([n/2,n/2+1])) = A(n/2,Index([n/2+1,n/2]));
			A(n,Index([n/2,n/2+1])) = A(n,Index([n/2+1,n/2]));
		}

		for (i = 1; i <= n; i++) {
			for (j = 1; j <= n; j++) {
				M(i,j) = A(i,j) + n*(A(j,i) - 1);
			}
		}
	} else {
		// odd order
		k = 0;
		for (i = -n/2; i <= n/2; i++) {
			for (j = 0; j < n; j++) {
				M(rem(j-i+n,n)+1,rem(j+i+n,n)+1) = ++k;
			}
		}
	}

	return M;
}