n-阶幻方问题

前言

这是很久以前做的一个实验的内容，觉得特别有意思，所以一直想发布出来，没想到拖着拖着就到现在了。

问题描述

一个 n 阶幻方是把从 1 到 n^2 的整数赶往一个 n 阶方阵，每一个数只出现一次，每一行、主副对角线的和都相等。

分析和本文基本概念

分析

据了解，4 阶幻方个数的基本型就有 880 个，通过旋转和反射总共可有 7040 个不同的形式的，5 阶幻方基本型有 275 305 224 个，6 阶幻方的个数非常之多，皮恩和维茨考夫斯基利用蒙特卡洛模拟和统计学方法，也只能获得一个大概的估计数字，其数量在 1.7743 x 10^19 ~ 1.7766 x 10^19 之间，可见随着幻方阶数的增加其数量将会快速地增长。同时 6 阶幻方的个数仅仅是理论值，其中的幻方未必就能利用现有的方法构造出来。但我们依旧希望能够找到尽可能多的幻方。

本文基本概念

为了能够在下文中方便地进行描述，这里先说明几个基本的概念：

1. 不存在 2 阶幻方；
2. 奇数阶幻方：阶数 n 为奇数的幻方；
3. 双偶数阶幻方：阶数 n 可以被 4 整除的偶数阶幻方，即 4K 阶幻方（如 n =4，8，12 ……的幻方）。
4. 单偶数阶幻方：阶数 n 不可以被 4 整除的偶数阶幻方, 即 4K+ 2 阶幻方（如 n =6，10，14……的幻方）。
5. 幻和：幻方行、列、主副对角线的和，为 n (nn+1)/2；
6. 互补数：在 n 阶幻方中，和为 n^2+ 1 的两个数为互补数；
7. 同构幻方：可以通过矩阵行列交换相互得到的幻方；
8. 映射：幻方左右或上下两侧关于对称轴对称；
9. 映射行、列：关于幻方水平对称轴对称的行，关于垂直对称轴对称的的列；
10. 映射操作：对映射列或映射行操作。

在通过查找大量的资料之后我决定采用罗伯法（解决奇数阶幻方，也叫连续摆数法），海尔法（解决双偶数阶幻方，也称象限对称交换法），斯特拉兹法（解决单偶数阶幻方），以及我通过研究罗伯法得到的素数阶幻方构造法，基本原理也是连续摆数法。此外，为了得到更多的幻方，我对以上四种方法多进行研究拓展。

奇数阶幻方：罗伯法（连续摆数法）

填写方法

把 1（或最小的数）放在第一行正中，并按以下规律排列剩下的 (n×n-1) 个数：

1. 每一个数放在前一个数的右上一格；
2. 如果这个数所要放的格已经超出了顶行那么就把它放在底行，仍然要放在右一列；
3. 如果这个数所要放的格已经超出了最右列，那么就把它放在最左列，仍然要放在上一行；
4. 如果这个数所要放的格已经超出了顶行且超出了最右列，那么就把它放在底行且最左列；
5. 如果这个数所要放的格已经有数填入，那么就把它放在前一个数的下一行同一列的格内。

图解说明如下：

方法的研究

经研究，发现罗伯法有一下规律：
将两个按 1~n^2 顺序填写如下的两个矩阵放在一起，并按副对角线方向取出五列，如图 1。按照 45123 排列可以得到图 2。很容易就可以发现，图 2 的每一列和主副对角线的和都是 5 *(5^2+1)/2，也就是幻和。同时也可以很容易的发现每两个相对中间的数字“13”中心对称的数字都是互补数，和为 5^2+1。按照罗伯法，“1”右边的每一列都相对左边的一列往上移动一行，相应的，为了保持对角线的和保持不变，又因为每两个相对中间的数字“13”中心对称的数字都是互补数，所以只需要对左边的列做映射操作即可，即分别下移动相应的行数。利用超出方框的数字补全缺少的位置即可得到一个幻方，如图 3 中的中间的一个矩阵。

得到这个结果之后我们开始设想，如果往上移一行可以使得每一行的和为幻和，那么往上移动两行是否也可以使得每一行的和也为幻和呢？在对这个 5 阶幻方进行往上移动 2~4 行得到的幻方进行计算得知当移动 2~3 行时都可以构成幻方。当移动 4 行时相当于用图 3 往下移动一行得到如图 4。分析可知向下移动一行是不可行的，由此可以得出往上移动 1~n- 2 行都可以成立，其实往上移动 n - 2 行，相当于往下移动 2 行就是已有的“马步法”构造素数幻方。于是利用计算机对 3~500 中的所有奇数进行往上移动 1~n- 2 行构成幻方的测试，发现 3~500 中的所有合数中仅有少数的情况可以构成幻方，而素数均可以构成。为了确认素数是否均可符合这个规律，我对 3~5000 以内的所有素数都进行以上的测试，均可构成幻方。由此推测所有大于 2 的素数均可构成幻方。这就是我所发现的 素数阶幻方构造法，通过这个方法一个素数 n 可以得到 n - 2 个幻方，非素数奇数可以得到 1 个幻方。

自创素数阶幻方构成法

该方法是基于罗伯法的基础上而得出的，也是连续摆数法。对于素数 n，通过这个方法可得到 n - 2 个幻方。该方法在罗伯法中进行了详细的介绍，这里就不在重复说明。

双偶数阶幻方：海尔法

所谓双偶阶幻方就是当 n 可以被 4 整除时的偶阶幻方，即 4K 阶幻方。以 8 阶幻方为例：

1. 先把数字按顺序填。然后，按 4×4 把它分割成 4 块
2. 每个小方阵对角线上的数字，换成和它互补的数。

单偶数阶幻方：斯特雷奇法

单偶数阶幻方就是当 n 不可以被 4 整除时的偶数阶幻方，即 4k+ 2 阶幻方（如 n =6，10，14……的幻方）。半阶数 m = n/2。

填写方法：
以 10 阶幻方为例。这时，k=2，m = 5。

1. 把 n 阶矩阵分为 A，B，C，D 四个象限（分别对应象限 2,1,3,4），这样每一个象限肯定是奇数阶，阶数为 m。按照奇数阶的方法用 1~(n/2)^2 在 A 象限填写数字，并用 A 象限初始化其他象限，其中 B = A+2m^2，C = A + 3m^2，D = A + m^2。如图 1。
2. （2）在 A 象限的中间行、中间格开始，按自左向右的方向，标出 k 格。A 象限的其它行则标出最左边的 k 格。将这些格，和 C 象限相对位置上的数，互换位置。如图 2。
3. （3）在 B 象限任一行的中间格，自右向左，标出 k - 1 列。(注：6 阶幻方由于 k -1=0，所以不用再作 B、D 象限的数据交换) 将 B 象限标出的这些数，和 D 象限相对位置上的数进行交换，就形成幻方。

方法的研究

对该方法生成幻方的过程，我们可以利用矩阵的加法及数乘进行实现（如下图）。其中由 0,1,2,3 组成的矩阵我们称之为叠加矩阵，另一个由奇数幻方生成的因数矩阵我们称之为源矩阵，最右边的矩阵我们称之为生成矩阵。

易知对这三个矩阵同时进行步骤 2,3 的操作等式仍是成立的，其中由于源矩阵的 A,B,C,D 象限都是一样的，在进行象限对称交换时，源矩阵没有发生变化，可知引起生成矩阵变化的是叠加矩阵。叠加矩阵变换如下图。

对叠加矩阵进行分析。对于数字 2 和数字 1 的矩阵（B 象限和 D 象限），由于是映射交换，所以无论交换哪 k - 1 列对每一行，每一列的和的作用是一样的。同时可以发现，交换 k - 1 列对主副对角线的影响都是交换了 k - 1 个数。因此，我们只需要对应地交换 B,D 象限的任意 k - 1 列都可以使得生成矩阵成为幻方。其中叠加矩阵的系数 m^2 可以使得得到生成矩阵的数字由 1~n^2 组成。

再对源矩阵进行分析。源矩阵的四个象限都是 m 阶的奇数阶矩阵，所有构成的源矩阵每一行、每一列、主副对角线的和都是 2 (m(mm+1)/2)。如果我们独立地构造这四个象限的 m 阶矩阵，所得到的源矩阵的每一行、每一列、主副对角线的和仍是 2 (m(mm+1)/2)，与数乘后叠加矩阵相加仍可以构成幻方。

通过这些变换我们可以得到的 n 阶幻方个数为 ((m – 2)^isPrimeNum(m))^4*C(m, k-1) 个，其中 m =n/2，k=n/4，isPrimeNum(m)是一个判断是否为素数的函数，如果是素数，函数值为 1，否则为 0。m- 2 是由素数法得到的奇数阶幻方个数，由于由四个独立的幻方，所以就有 4 次方，C(m, k-1) 组合是表示在 B，D 象限中选择 k - 1 列进行交换可以得到的种类数。由于这些是后来进一步拓展的，所有在之前实现的代码是分别生成 A,B 两个象限的幻方，并利用这两个幻方分别生成 C,D 两个换，同时 B,D 象限我们交换的是自右向左的 k - 1 列。这样我们得到的幻方个数偏少，为 ((m – 2)^isPrimeNum(m))^2 个。

注：C(n, m) 为组合公式，n > m，表示从 n 中取 m 个的组合。C(n,m)=n!/(m!(n-m)!) 即 C(n,m)=A(n,m)/m!。*

通过矩阵行列交换、旋转实现幻方个数的暴增

通过对这三种方法生成的幻方进行进一步的分析，我们发现连续摆数法（包含罗伯法和素数阶幻方构成法）生成的奇数阶幻方和海尔法生成的双偶数阶幻方都有以下规律：幻方中的每两个关于幻方中心中心对称的两个数均为互补数。因而关于幻方中心中心对称的两列或者两行上的对应位置数字也是互补数。因此关于幻方对称轴映射交换任意行列都是可以构成幻方的（如下图）。

行列交换

行列映射交换的方式又可以有多种，总的来说可以概括为以下的四种（由于行列的交换都是一样的，这里以列交换进行分析）：左右对称交换，左右交叉交换，单侧映射交换，混合交换（即左右对称交换、左右交叉交换和单侧交换混合）。

当阶数很大的时候，混合交换的情况就会变得非常复杂，似乎是无法找到所有的情况。但是经过我们分析，所有的情况都可以通过以下两个步骤得到：

第一步：确定左右交换的列（行），进行左右对称交换；
第二步：对对称交换后的幻方进行单侧映射交换。

例子：对行 12345 操作，使得到的结果 3614725
第一步：选出需要交换的列，由交换后的结果知，左侧需要和右侧交换的有：2（由映射关系 2,6 交换）

第二步：左侧单侧交换为 3,6,1（由映射关系得，右侧单侧交换为 7,2,5）。

这样就可以得到最终的结果了。
这两个步骤其实就是组合排列的过程。第一步，利用组合对左侧找到所有组合情况，即所有左右交换的情况，包含交换 0 列，1 列……n/ 2 列，对组合结果进行左右对称交换，可得到 2n/ 2 种。第二步，对左侧进行全排列，找到所有排列的情况，即所有单侧映射交换的情况，进行单侧映射交换，可得到(n/2)! 种。由于组合排列得到的结果均是不同的，所有这里得到的所有结果都是唯一的。这样，一个奇数阶幻方或一个双偶数阶幻方就可以拓展得到 2n/2*(n/2)! 种。

旋转

这里我们以一个 2 阶的矩阵来分析。按照逆时针旋转，每次旋转 90 度可以得到下列四个矩阵。

通过行列交换可以由下图中的 a 图得到 c 图，但无法得到 b，c 图，即无法使得同列关系的数字 1,3 变为同行关系。但是可以通过 a，c 图顺时针旋转分别得到 b，d 图。一次我们可以通过一个矩阵行列交换，并用行列交换后的结果顺时针旋转 90 度得到它的四个旋转后的矩阵，这样得到的矩阵不会重复。那么，一个奇数阶幻方或一个双偶数阶幻方就可以拓展得到 2n/2(n/2)!2 种。由于单偶数阶幻方没有行列交换，因此，可以对一个单偶数阶幻方进行三次旋转得到 4 个不同幻方（包括自身）。
由于这个结论也是后来才得到，所以在我们的代码实现总没有进行旋转。

行列交换的实现

要实现 5 中的幻方个数暴增，必须解决的问题就是行列交换问题，也就是组合排列问题。

组合问题

这里我们将组合问题转化为图的路径遍历问题，在 n 个数字中选出 m 个数的所有组合，相当于在一个这样的图（下面以从 1,2,3,4 中选出 3 个数字为例说明），求从 [1,1] 位置出发到达 [m,x] (m<=x<=n) 位置的所有路径：

上图是截取 n×n 右上对角矩阵的 m 行构成，我们要求的所有组合就相当于从第一行的第一列元素 [1,1] 出发，到第三行的任意一列元素作为结束的所有路径，规定每走一步需跨越一行，并且从上一行的任何一个元素到其下一行中列处于其右面的任何一个元素均有一路径相连，显然任一路径经过的数字序列就对应一个符合要求的组合。由于组合的数组不一定是按照 1~n 这样升序的，但我们可以将数字存放到一个一维数组中，并对其下标进行组合即可。
这里用回溯的方法实现。由 [1,1] 开始，[1,1] -> [2,2] -> [3,3]，得到一组，然后 [2,2] 位置指向[3,4], 这得到另一组，第三行没法再往右移动，回溯到第二行往右移动，路径变成[1,1] -> [2.3] -> [3,4]，按照这个规律下去，知道第一行往右移到最后一列。这样就可以得到所有合要求的组合了。

排列问题

这里用递归实现，下面以 1,2,3,4,5 的全排列为例说明。

1. 首先看最后两个数 4, 5。它们的全排列为 4 5 和 5 4, 即以 4 开头的 5 的全排列和以 5 开头的 4 的全排列。由于 4 或者 5 的全排列就是其本身，从而得到以上结果。
2. 再看后三个数 3, 4, 5。它们的全排列为 3 4 5、3 5 4、4 3 5、4 5 3、5 3 4、5 4 3 六组数。即以 3 开头的和 4,5 的全排列的组合、以 4 开头的和 3,5 的全排列的组合和以 5 开头的和 3,4 的全排列的组合。
3. 从而可以推断，设一组数 p = {r1, r2, r3, … ,rn}, 全排列为分别以各个数字为第一个数字的全排列的总和。

实现上其实就是，每一次都将所有的元素和第一个元素交换即可

如图 1，每一次递归都会 start 指向的位置都会移动，数组的长度也会跟着变化。再看图 2，由于每个数字都可以作为当前递归层次的开始位置，因此可以得到如图 2 的情况。在第一层递归中，开始的位置可以是 1~5,，上一层为 1 时候，第二层中，开始的位置有 2~5，按照这个规律下去就可以得到所有的全排列。
由以上的描述可以知道每生成一个 n 个数字的全排列需要进行 n 次递归，也就是有 n 层递归。对于 n 阶的幻方，由于我们排列时只需要对左侧排列即可，也就是需要对 n / 2 个数字进行排列，这样就会有 n / 2 层的递归，对于阶数很大的幻方这无疑会成为一个限制。所以下一步的优化就是需要将递归实现全排列替换为非递归的全排列方法，这个暂时还没有实现。

代码最终结构及幻方个数计算

代码最终结构

已实现的代码的结构：

幻方个数计算

通过以上这些过程，可以得到幻方个数如下：

如果对得到的幻方都进行顺时针旋转 90 生成另一个幻方，同时利用四个奇数阶幻方对单偶数阶幻方进行优化，这样得到的新的幻方个数计算公式如下：

代码实现已经上传到了 github：MagicSquare