周小安+赵宇

摘要： 关键词： 中图分类号： 文献标志码： A文章编号： 2095-2163

Abstract： Magic square originated thousands of years ago in China，which is one of the research objects of combinatorial designs. In a square matrix，if the sums of every row，every column and two diagonals are equal， this square matrix is called magic square. Currently， the achievements in the research of magic square have been widely used in image encryption，image scrambling，etc. This paper introduces a new continuous numbering method to structure arbitrary odd order magic square，which is more flexible than traditional construction methods.The new method can not only increase the forms of odd order magic square， but also parallelly structure magic square and reduce the time of structuring it.

0引言

幻方具有非常悠久的历史，虽然古老但却是最流行的数字游戏之一。一个n阶幻方的构造方法是将整数1、2、3、…、n2共n2个连续自然数，填入到n × n的方阵中，使得每行、每列以及每条对角线上的数字之和都等于同一个常数。幻方是组合设计的重要研究对象之一[1]。

《易经》中记载的洛书是世界上最早的幻方，之后幻方逐渐传入世界各地，引起了广泛关注，并在许多方面取得了研究成果。幻方可应用于哲理思想的研究、美术设计、智力开发、科学技术等方面[2]。尤其是在科学技术中，随着计算机的快速发展，幻方在信息隐藏[3]、数字图像[4]、认证与加密[5-6]、量子信息[7]等方面应用前景十分广阔。

1使用改进的连续摆数法构造奇数阶幻方

1.1连续摆数法

连续摆数法是较为古老的幻方构造方法，适用于奇数阶幻方的构造，一般把数字1放在幻方中第一行最中间的方格中，然后从这里开始，按对角线方向（比如说按从左下到右上的方向）把其余各数按从小到大的顺序依次放入其它格子中，如果到达顶部，则转向底部；如果到达右侧，则转向左侧；如果需要放入的格子中已经有数存在或者到达了右上角，则需要退至前一格的下方[1]。

按照上述法则建立的5阶幻方如图1所示。

这个方法可以推广到一般情况，即起始数1不一定要放在第一行的中间，而下一个数字也不一定要放在上一个数的右上方格中，为此本文提出了改进的连续摆数法来构造奇数阶幻方。

1.2改进的连续摆数法

1.2.1算法设计步骤

改进的连续摆数法构造n阶幻方的步骤如下：

步骤1以奇数阶幻方中心位置为坐标原点，建立坐标系；

步骤2将幻方中任意一格的位置用坐标来表示，坐标范围从（-n-12，-n-12）到（n-12， n-12）；

步骤3定义“起始向量”，表示1的位置；

步骤4定义“偏移向量”，表示一个数的位置到下一个数位置所指的方向，依次按顺序填写其它数字；

步骤5当遇到要填入的格子中已经被其它数字占据，用“中断向量”重新计算该数字的坐标并填入计算后的空格；

步骤6重复步骤4和步骤5，直到所有数字均填入方格中。

1.2.2算法设计过程

这里，下面将详述给出改进连续摆数法构造5阶幻方的过程。先以5阶幻方中心位置为坐标原点，建立坐标系，每一点的坐标如图2所示。

接着确定1的位置，连续摆数法是基于对称性的方法，即以5阶幻方为例，1～25的中间数字是13，所以假设要让13固定在幻方的中心点，那么1就只能选择除中心点外的其它位置。为此，研究定义一个起始向量（x0，y0），表示1的起始位置，x0、y0满足条件（x0，y0） ≠ （0，0）。

确定好1的位置后，还需要定义一个偏移向量（u，v），表示一个数的位置到下一个数位置所指的方向。当然，偏移向量的选择在理论上可以是任意的，但是却需要满足如下条件：

首先，u和v都不能等于0，否则无法产生偏移；其次，必须要保证（1+n2）/2固定在幻方的中心点，那么选取偏移向量时就要避免让（1+n2）/2之前的数字占据中心位置，就需要据此原理求解数个特征方程，这里暂时不考虑这个问题。至此，研究将针对起始向量和偏移向量展开全面的分析论述。

假设阶数n=5，起始向量（x0，y0） = （0，-1），偏移向量（u，v） = （2，1），那么1的位置如图3所示，2的坐标由1的坐标加上偏移向量所得，即（0，-1） + （2，1） = （2，0），将2填入到图3中对应位置，3、4、5采用相同的方法依次填入。需要說明的是，3的坐标由2的坐标加上偏移向量所得，即（2，0） + （2，1） = （4，1），在幻方中显然没有对应的坐标存在，这里需要把超出范围的坐标进行“模n加”运算使该坐标在该幻方中寻得对应位置。具体来说，就是坐标（4，1）的横坐标超出范围，将（4，1） - （5，0） = （-1，1），于是找到3在幻方中的位置，填入图3中。对4和5的处理也采取和3类似的方式。endprint

下一步就是确定6的位置。无论是从计算6的坐标位置还是从图3直接来看，研究发现不管起始向量和偏移向量的取值是多少，理论上计算得到的6的位置上已经有1存在了，这和连续摆数法的周期性有关。那么要想确定6的位置，必须引入一个中断向量（xt，yt），表示发生冲突时的偏置量。具体来说，就是将数字正准备填入某一方格时，却遭遇该方格中已存在数字的情形，这时即需改变该数字坐标的计算方法，将中断向量代替偏移向量，重新计算该数字的坐标。

继续用图3中的例子来说明，当用偏移向量计算得到的6的位置里已经存在1了，这时需要将中断向量代替偏移向量，重新计算6的坐标，即用5的坐标加上中断向量。那么如何确定中断向量的大小？

经过推导可知，中断向量的值与起始向量有关，中断向量（xt，yt）满足方程（1）：（xt，yt）=2×（x0，y0）（1）在该例中，（x0，y0） = （0，-1），所以（xt，yt） = （0，-2），这时计算得到的6的坐标为5的坐标加上中断向量，即（-2，-2） + （0，-2） = （-2，-4），再利用“模5加”运算得到6的坐标为（-2，1），填入幻方中。

紧接着，继续用偏移向量完成其后顺承数字的填写，直到遇到下次被“占位”时才使用中断向量。于是可以将7、8、9、10依次填入幻方中，如图4所示。

计算该5阶幻方每行、每列及对角线数字的和均为65，符合幻方的规律。

经过尝试和推导，可以将改进的连续摆数法推广到任意奇数阶幻方，用来构造阶数更多的幻方。由于起始向量和偏移向量是随机选取的，两者组合后可以构造出许多种完全不同的幻方图样，如果事先不知道這两个向量，就极难复制出一样的幻方。

2改进连续摆数法的通项公式

如果从1到n2逐个去填入幻方的话，显然可得时间复杂度即为O（n2），而且需要知道一个数的坐标才能推算出下一个数的坐标，所以有必要推导出通项公式，可以直接计算出一个数的坐标，直接填入幻方中，提高构造幻方的效率。

参考文献：

[1]吴鹤龄. 幻方及其他-娱乐数学经典名题[M]. 北京：科学出版社，2005.

[2] 高治源. 幻方的应用前景[J]. 延安教育学院学报，2000（1）：44-46.

[3] 丁玮，齐东旭. 数字图像变换及信息隐藏与伪装技术[J]. 计算机学报，1998，21（9）：838-843.

[4] ARONOV B， ASANO T， KIKUCHI Y， et al. A generalization of magic squares with applications to digital halftoning[M]//FLEISCHER R， TRIPPEN G：ISAAC 2004， LNCS 3341.Berlin/Heidelberg： Springer-Verlag， 2004：89-100.

[5] XIE T， CHEN H， KANG L. A unified method of magic square-based two-way certificate authentication and key transferring：China，02114288.2 [P].2002.

[6] XIE T. A magic square based signing method for identification and authentication：China， 200410046922.2[P].2004.

[7] RAMZAN M， KHAN M K. Distinguishing quantum channels via magic squares game[J]. Quantum Information Processing， 2010， 9（6）： 667-679.

[8] ZHANG Jie，LU Yang， YANG Shu， et al. NHPPbased software reliability model considering testing effort and multivariate fault detection rate[J]. Journal of Systems Engineering and Electronics， 2016，27（1）： 260-270.

[9] 马飒飒，王光平，赵守伟. 基于时间序列的软件可靠性预测模型研究[J]. 计算机工程与设计，2007，28（11）：2520-2523.

[10]贾治宇，康锐. 软件可靠性预测的ARIMA方法研究[J]. 计算机工程与应用，2008，44（35）：17-19.

[11]何书元. 应用时间序列分析[M]. 北京： 北京大学出版社，2003.

[12]黄雄波. 多周期时序数据的傅氏级数拟合算法[J]. 计算机系统应用，2015，24（7）：142-148.

[13]黄雄波. 基于自相关函数的非平稳时序数据的辨识改进[J]. 微型机与应用，2016，35（13）：10-14，18.

[14]刘思峰，杨英杰. 灰色系统研究进展（2004-2014）[J]. 南京航空航天大学学报，2015，47（1）：1-18.

[15]党耀国，王俊杰，康文芳. 灰色预测技术研究进展综述[J]. 上海电机学院学报，2015，18（1）：1-7，18.

[16]郑咸义，姚仰新，雷秀仁，等. 应用数值分析[M]. 广州： 华南理工大学出版社，2008.

[17]李庆扬，关治，白峰彬. 数值计算原理[M]. 北京： 清华大学出版社，2005.endprint