幻方，也有人称之为魔方，就是横竖斜和都相等的那种方阵。

小学6年级在奥数课上老师告诉我奇数阶幻方的通解，不知道的小朋友们，我在这科普一下吧：

比如说三阶幻方，先向外翻折扩展，然后按上图左二的规律，按顺序写上1-9的数字，接下来幻方之外的数，按左往右仍，右往左仍，上往下扔，下往上扔的规律填进幻方，将其余的删去，就得到一个横竖斜都等于15的幻方了！

下图是五阶幻方的解法，方法相同，只是规模大了点。
七阶幻方如下：（唉，上面那种做图太累，后面的图就来自于互联网了。。）

只要按照这个方法，无论多少阶，只要是个奇数，都可以画得出来，至少一个！你可以奸诈一点，比如说画好菱形后，1的起始位置是可以换的，写的方向也是可以换的，但是最后出来的幻方本质上是一样的。。

对于偶数呢，最小是4阶的，四阶的幻方老师也讲了一个解法，就是大对角线换，小对角线也换。步骤如下：

先按顺序写出1-16的数在4阶幻方里面，如下：

接下来所谓的大对角线换，小对角线换就是1和16换，4和13换，6和11，7和10，换完就出来了：

横竖斜都是34。

然后问题就来了，有没有办法可以解出任意高偶数阶的幻方的方法呢？

我曾经很傻很天真的试图把4阶这种换对角线的方法推广到6阶，但是怎么弄都未果，估计这种方法对于4阶只是种巧合吧。

后来大学玩matlab后，发现matlab里面函数magic可以输出任意阶的幻方，哦，soga，原来真的有的啊。

后来我就对着matlab里面magic的源文件写出了这个C++版本，只是为了巩固自己对四阶的理解罢了。

然后下面整理一下一般的偶数阶幻方的解法，解法来源于互联网。

首先一般的偶数阶解法都是把偶数分成两种，4,8,12,16这种4m的双偶数和6,10,14这种4m+2的单偶数，一般的解法都是分开来两类的，包括matlab里面的magic函数，不过查了一下也有很多大牛研究出了统一解法，更有大神把奇偶阶全部同意了，膜拜ing。。。

双偶数解法：偶数阶下面先讲简单的双偶数解法，看了很多解法，但是最后发现了一个通解，网上看到的大部分解法都是这个通解的特例。

首先呢，如下图所示，先把n阶幻方分成4个小块，对于左上角那个你任意的把一半放个填成灰色，但是有一个约束条件，就是左上角这个小块中每一行每一列都要只有n/4个灰色的。然后呢，右上的那个小块的填色方案就是左上填色方案的左右镜像对称，左下的就是左上天色方案的上下镜像对称，自然，右下就是左上的中心对称了。如下图所示：

然后呢，你把1-n2这么多个数按顺序填进白色的格子里去，灰色的部分要留着。如下面左图所示：

之后呢，把剩下的没填的数反过来填进去，也就是从右下到左上的顺序，填完双偶数阶幻方就出来了。

现在我们来讨论一下这种方法，首先看我们原本的四阶幻方的解法，有没有发现其实和这种方法是一个东西。

然后再看看双偶数阶的另一种解法，比如说下面这个8阶幻方：

这里的解法呢，就是把整个幻方分成2×2个4×4的小块，按顺序填好1-64个数，然后每个4×4小块的对角线上的数不变，其余的数做中心对称。

再看看下面这个：

12阶，分成3×3个4×4的小块，和之前一样，按顺序填好数，然后每个4×4小块的对角线上的数不变，其余的数做中心对称。

虽然和我最开始的那种分法不一样，但是你仔细一想，其实是完全一样的，只是他的填色方案是固定的一种模式而已。

还有一种说法是每个小块对角线上的数换成互补的那个数，其实本质还是一样嘛。

下面是一个双偶数的matlab程序，我填色方案用时是国际象棋棋盘那种黑白相间。

function a = hf_4m(n)
flag = zeros(n/2,n/2);
flag(1:2:n/2,1:2:n/2) = 1;
flag(2:2:n/2,2:2:n/2) = 1;
flag = [flag fliplr(flag);flipud(flag) flipud(fliplr(flag))];
a = reshape(1:n^2,n,n)';
a = a .* flag;
a = reshape(a',1,n^2);
blank_idx = find(a==0);
number_left = (1:n^2) .* (a==0);
number_left = fliplr(setdiff(number_left,0));
a(blank_idx) = number_left;
a = reshape(a,n,n)';

单偶数解法：下面来看看单偶数的解法，这种现在主要有两种方法，分区法和易位法。其中呢，分区法也有两种。

先说分区法，首先呢就是把方阵划分成下面A,B,C,D四块，因为是单偶数，所以每一块必然是个奇数幻方。
然后把1~n2/4这些数组成的奇数阶幻方算出来，填进A里面，然后接下来的n2/4的幻方填进D里面，（其实有个很简单的方法，就是把A里面的每个数加上n2/4就可以了），再把D里面的加上n2/4放到B里面，最后那些放到C里面。下面是10阶幻方的一个例子：

然后下面的东西有点拗口，但是细细读就会明白了：

先假设阶数是4k+2，那么k=(n-2)/4，然后下面是第一种方法：

从A小块的中间行中间格开始（上图中的13），向右找k个数（包括中间行中间格那个），和C小块的相应位置的数换位。A小块的其他行（也就是除了最中间那一行）从最左开始数出k个数，和C中相应位置的数换。

B小块中间列开始，向左数K-1列出来（当然也包括B小块中间那一列），然后这些列和D小块中相应未知的数换位。（6阶k-1=0，就不用了）

然后就完成了。

14阶幻方的换位方式如下：

这种方法的matlab函数如下：

function a = hf_4m_2(n)
a = zeros(n,n);
a(1:n/2,1:n/2) = magic(n/2);
a(n/2+1:n,n/2+1:n) = a(1:n/2,1:n/2) + n^2/4;
a(1:n/2,n/2+1:n) = a(1:n/2,1:n/2) + n^2/4*2;
a(n/2+1:n,1:n/2) = a(1:n/2,1:n/2) + n^2/4*3;
m = (n/2-1)/2;

temp = a((n/2+1)/2,(n/2+1)/2:(n/2+1)/2+m-1);
a((n/2+1)/2,(n/2+1)/2:(n/2+1)/2+m-1) = a((n/2+1)/2+n/2,(n/2+1)/2:(n/2+1)/2+m-1);
a((n/2+1)/2+n/2,(n/2+1)/2:(n/2+1)/2+m-1) = temp;

temp = a(setdiff(1:n/2,(n/2+1)/2),1:1+m-1);
a(setdiff(1:n/2,(n/2+1)/2),1:1+m-1) = a(setdiff(1:n/2,(n/2+1)/2)+n/2,1:1+m-1);
a(setdiff(1:n/2,(n/2+1)/2)+n/2,1:1+m-1) = temp;

if(m>1)
    temp = a(1:n/2,n*3/4+1/2-m+2:n*3/4+1/2);
    a(1:n/2,n*3/4+1/2-m+2:n*3/4+1/2) = a(1+n/2:n,n*3/4+1/2-m+2:n*3/4+1/2);
    a(1+n/2:n,n*3/4+1/2-m+2:n*3/4+1/2) = temp;
end

然后还有一种换位方法，A小块中间那一行第2列开始往右数k个数，和C小块中相应位置的数换位，A小块中其余行都从最左开始向右数k列，这些数也和C小块的做交换。B小块中，从最右开始向左数k-1个列，与D中相应位置的数换位，结果也是一样的。这种方法matlab代码如下：

function a = hf_4m_2(n)
a = zeros(n,n);
a(1:n/2,1:n/2) = magic(n/2);
a(n/2+1:n,n/2+1:n) = a(1:n/2,1:n/2) + n^2/4;
a(1:n/2,n/2+1:n) = a(1:n/2,1:n/2) + n^2/4*2;
a(n/2+1:n,1:n/2) = a(1:n/2,1:n/2) + n^2/4*3;
m = (n/2-1)/2;

temp = a((n/2+1)/2,2:2+m-1);
a((n/2+1)/2,2:2+m-1) = a((n/2+1)/2+n/2,2:2+m-1);
a((n/2+1)/2+n/2,2:2+m-1) = temp;

temp = a(setdiff(1:n/2,(n/2+1)/2),1:1+m-1);
a(setdiff(1:n/2,(n/2+1)/2),1:1+m-1) = a(setdiff(1:n/2,(n/2+1)/2)+n/2,1:1+m-1);
a(setdiff(1:n/2,(n/2+1)/2)+n/2,1:1+m-1) = temp;

if(m>1)
    temp = a(1:n/2,n:-1:n-m+2);
    a(1:n/2,n:-1:n-m+2) = a((1:n/2)+n/2,n:-1:n-m+2);
    a((1:n/2)+n/2,n:-1:n-m+2) = temp;
end

这两中方法为什么可行我还没仔细研究，但是刚刚编程发现一个很神奇的现象，就是第一种方法的BD小块交换规则配上第二种方法的AC小块交换规则，也是可以的。。囧。。

对于单偶数的幻方，还有一种杨辉创造的二阶方阵易位法（我发现杨辉老兄很喜欢玩数阵）。。
对于n = 4m+2阶幻方，先用奇数阶的方法做出一个2m+1阶幻方来，然后把1~n2那么多个数4个一组，分成(2m+1)2个组，{1,2,3,4}{5,6,7,8}{9,10,11,12}…分别称为第1组，第2组，第3组…第(2m+1)2组。

接下来那每一组四个数按下面的方法放入2×2的方阵中：

然后把之前那个2m+1阶幻方，每个位置上的数如果是i，那么就换成第i组2×2方阵，这样就有了一个n×n的方阵了，但是这个方阵还不是幻方，需要再修正。

我们继续讨论刚才那一个2m+1阶幻方，假设我们n=14，那么2m+1 = 7，对于下图中这个7阶的方阵，我们把倒数第二行染绿，然后从中间那一行开始向下知道倒数第三行为止全部染蓝，如果中间那一行就是倒数第二行，那么不染蓝。

接下来在把最左和最右两列的染色向下拉一格。如下图：

我们知道对应于刚刚做出来的那个n×n的方阵，每2×2方阵，四个数对应于上图的一个格。我们现在做如下操作，如果是绿色的格子，那么2×2方阵的最下面两个数交换，如果是蓝色格子的话，2×2方阵不仅下面两个数交换，而且上面两个数也要交换。

下面举个例子：

对于14阶幻方，先生成一个7阶幻方

然后把相应位置填上相应的2阶方阵：
填好色：

绿色格子下面两个数换位，蓝色的上下都换，就得到结果啦~~~

杨辉易位法代码如下：

function a = hf_4m_2_yiwei(n)
h = magic(n/2);
a = zeros(n,n);
for i = 1 : n/2
    for j = 1 : n/2
        a(i*2-1:i*2,j*2-1:j*2) = [2 3;4 1] + (h(i,j)-1)*4;
    end
end
flag = zeros(n/2,n/2);
flag(n/2-1,:) = 1; %%下面两个互换
if(n > 6)
    flag((n/2+1)/2:n/2-2,:) = 2;%%上面下面都要换
end
flag(2:n/2,[1 n/2]) = flag(1:n/2-1,[1 n/2]);
for i = 1 : n/2
    for j = 1 : n/2
        if(flag(i,j) >0)
            temp = a(i*2,j*2-1);
            a(i*2,j*2-1) = a(i*2,j*2);
            a(i*2,j*2) = temp;
        end
        if(flag(i,j) == 2)
            temp = a(i*2-1,j*2-1);
            a(i*2-1,j*2-1) = a(i*2-1,j*2);
            a(i*2-1,j*2) = temp;
        end
    end
end

哟西，好了，终于写的差不多了。自个研究了一下，收获颇多。不过幻方可不仅仅是构造那么简单，以前看的一本书里面有各种变态的幻方，什么切尾幻方什么的。。还有很多数学上的东西，下面提问，请证明：偶数阶幻方行列式值一定是0！

下面附上一个很多年前改写matlab的magic函数的C++代码：

#include<iostream .h>
#include<iomanip .h>
class magic
{
public:
    int **m;
    magic();
    magic(int n);
    magic(magic &c);
    void show();
    int size;
};

magic::magic(magic &c)
{
    int n=c.size;
    size=n;
    m=new int*[n];
    for(int l=0;l<n ;l++)
    {
        m[l]=new int[n];
        for(int k=0;k<n;k++)
        {
            m[l][k]=c.m[l][k];
        }
    }

}

magic::magic(int n)
{
    size=n;
    m=new int*[n];
    for(int l=0;l<n;l++)
    {
        m[l]=new int[n];
        for(int k=0;k<n;k++)
            m[l][k]=0;
    }
}

void magic::show()
{
    for(int l=0;l<size;l++)
    {
        for(int k=0;k<size;k++)
            cout<<setw(3)<<m[l][k]<<" ";
        cout<<endl;
    }
    cout<<endl;
}

magic creatmagic(int n)
{
    int l,k;
    if(n%2==1)
    {
        magic i(n),j(n),m(n),a(n),b(n);
        for(l=0;l<n;l++)
            for(k=0;k<n;k++)
            {
                j.m[l][k]=k+1;
                i.m[l][k]=l+1;
            }
            for(l=0;l<n;l++)
                for(k=0;k<n;k++)
                {
                    a.m[l][k]=(i.m[l][k]+j.m[l][k]-(n+3)/2)%n;
                    if(a.m[l][k]&lt;0) a.m[l][k]+=n;
                    b.m[l][k]=(i.m[l][k]+j.m[l][k]*2-2)%n;
                    if(b.m[l][k]&lt;0) b.m[l][k]+=n;
                    m.m[l][k]=n*a.m[l][k]+b.m[l][k]+1;
                }
                return m;
    }
    else if(n%4==0)
    {
        magic i(n),j(n),d(n),m(n);
        for(l=0;l<n;l++)
            for(k=0;k<n;k++)
            {
                j.m[l][k]=k+1;
                i.m[l][k]=l+1;
                d.m[l][k]=(   ((i.m[l][k]%4)/2)  ==  ((j.m[l][k]%4)/2)  );
                m.m[l][k]=l*n+k+1;
                if(d.m[l][k]==1) m.m[l][k]=n*n+1-m.m[l][k];
            }
            return m;

    }
    else
    {
        int p=n/2;
        magic t(creatmagic(p)),m(n);
        for(l=0;l<n;l++)
            for(k=0;k<n;k++)
            {if(l<n/2&&k<n/2)
            m.m[l][k]=t.m[l][k];
            else if(l<n/2&&k>=n/2)
                m.m[l][k]=t.m[l][k-n/2]+2*p*p;
            else if(l>=n/2&&k</n><n /2)
                m.m[l][k]=t.m[l-n/2][k]+3*p*p;
            else
                m.m[l][k]=t.m[l-n/2][k-n/2]+p*p;
            }
            int e=(n-2)/4;
            //cout<<e;
            for(l=0;l<n;l++)
            {
                if( ( l>=0 && l< =(e-1) )  ||  (  l>=(n+1-e)  ))
                {
                    for(k=0;k</n><n /2;k++)
                    {
                        int temp=m.m[k][l];
                        m.m[k][l]=m.m[k+n/2][l];
                        m.m[k+n/2][l]=temp;
                    }
                }

            }
            int temp=m.m[e][0];
            m.m[e][0]=m.m[e+n/2][0];
            m.m[e+n/2][0]=temp;

            temp=m.m[e][e];
            m.m[e][e]=m.m[e+n/2][e];
            m.m[e+n/2][e]=temp;
            return m;
    }
}

void main()
{
    float m;
    int n;
    while(1)
    {
        cout<<"input the size:";
        cin>>m;
        n=int(m);
        if(n==2||n< =0)
            continue;
        magic result(creatmagic(n));
        result.show();

        cout<<"每一列和为：\n";
        for(int l=0;l<n;l++)
        {
            cout<<l+1<<":";
            int sum=0;
            for(int k=0;k<n;k++)
            {
                sum+=result.m[l][k];
                cout<<result.m[l][k];
                if(k!=n-1)
                    cout<<"+";
            }
            cout<<"="<<sum<<endl;
        }

        cout<<"每一行和为：\n";
        for( l=0;l<n;l++)
        {
            int  sum=0;
            cout<<l+1<<":";
            for(int k=0;k<n;k++)
            {
                sum+=result.m[k][l];
                cout<<result.m[k][l];
                if(k!=n-1)
                    cout<<"+";}
            cout<<"="<<sum<<endl;
        }

        int sum=0;
        cout<<"主对角线和为：";
        for( l=0;l<n;l++)
        {
            sum+=result.m[l][l];
            if(l==n-1)
                cout<<result.m[l][l];
            else
                cout<<result.m[l][l]<<"+";
        }
        cout<<"="<<sum<<endl;

        sum=0;
        cout<<"次对角线和为：";
        for( l=0;l<n;l++)
        {
            sum+=result.m[l][n-1-l];
            if(l==n-1)
                cout<<result.m[l][n-1-l];
            else cout<<result.m[l][n-1-l]<<"+";
        }
        cout<<"="<<sum<<endl;
    }
}

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com