幻方 幻方的故事与求解

图9

现在解释一下MATLAB代码，MATLAB文件附后。以下注释不包含在文件中，下面红色部分为注释，特此说明:

以下是总和为0的计算

a =零(15，19)；

a(1，1:3)= 1；

a(2，4:7)= 1；

a(3，8:12)= 1；

a(4，13:16)= 1；

a(5，17:19)= 1；

a(6，[1 48)]= 1；

a(7，[2 5 913)]= 1；

a(8，[3 6 10 14 17)]= 1；

a(9，[7 11 1518)]= 1；

a(10，[12 1619])= 1；

a(11，[8 1317])= 1；

a(12，[4 9 1418)]= 1；

a(13，[1 5 10 15 19)]= 1；

a(15，[2 6 1116])= 1；

a(14，[3 712])= 1；上面实现了一个矩阵，大小15*19。注意MATLAB变量是区分大小写的。这里我用小写的A来表示A矩阵

排名(a)；计算A矩阵的秩，运行后A的秩为12

a = rref(a)；矩阵按行简化，得到以下结果:

图10

aEx=[aones(15，1)* 38]；AX=b非齐次方程对应的增广矩阵

AEx = rref(AEx)；同上，对增广矩阵进行行变换。根据运行结果，

秩(A)=秩(AEx)，所以非齐次方程有解。AEx矩阵的内容如下A和B。把B放在A后面才能得到AEx

b =[-1-1-1-2-1-1-1-1；

1 0 1 1 0 0 0;

0 1 0 1 1 1 1;

1 1 0 1 0 0 0;

0 1 1 1 1 1 0;

-1 -1 -1 -1 -1 0 0;

0 -1 0 -1 0 -1 0;

0 0 1 1 1 1 1;

-1 -1 -1 -1 0 -1 0;

1 0 0 0 0 0 0;

0 1 0 0 0 0 0;

0 0 0 0 -1 0-1;

0 0-1 -1 -1 0 0;

0 0 1 0 0 0 0;

0 0 0 1 0 0 0;

0 0 0 0 1 0 0;

0 0 0 0 0 -1-1;

0 0 0 0 0 1 0;

0 0 0 0 0 0 1;];

B矩阵是手工输入的，即把上图简化矩阵A的第10、11、14、15、16、18、19列向右移动，改变它们的个数得到的列向量构成矩阵。其实B是基本解。

b =[76 0-38 0-38 38 38-38 38 0 38 38 0 0 38 0 0 38 0 0]'；

b是从上面计算出来的AEx变量最右边一列。因为是按顺序写的代码，所以先运行，看了一下AEx结果，然后直接在这里输入。

x =[1 2 5 710 3 8]'；

x =[3 3 7 4-10-12 21]'；这个x是由7个任意数据x10，x11，x14，x15，x16，x18，x19组成的列向量。这里给出的两个X验证数据是我前面图7验证-验证-2中的数据。

XZero = B * x；这里是为了验证“使用任何软件MATLAB、SYMPY、JAVA填充实数，使得每行或每列用0填充的数字是可重复的。”B*x实现了图6中的计算，结果保存在XZero变量中，一次得到所有需要填充的数字。。。尝试更改x的值。

下面是总和为38的计算。

temp = 1:19；

template = round(temp’)；根据我以前的设计思想，有必要用穷举法来计算一个满足非齐次线性方程组的可能结果。如果结果正确，解必须是1到19的整数。当然顺序是可以打乱的，所以到时候会把计算结果排序，和这里的模板进行对比。模板的值是从1到19的升序整数。

Tic开始计时。根据前面的回答，需要尝试A(19，7)=19*18*17*16*15*14*13这么多次，大概是2.5亿个周期，所以这里大概需要5~10个小时才能计时。取决于计算机配置和程序运行方式。

命中=0检测结果，加1

counts = 0；周期时间的统计累积

XTable =零(19，20)；所需填充编号的结果变量。每个结果都以列向量的形式存储。假设有20个结果，得到20列。。其实本质解决方案只有一个，一共12种形式。原因很简单。如果你想象你做了一个解决方案，你可以把整张脸翻过来，填充数字的方式就会改变。这叫对称。它还可以连续旋转60度，所以总共有12种解决方案。

其实最重要的是用图8对x10，x11，x14，x15，x16，x18，x19(1-19的非重复整数任意填充)的所有值计算结果，判断合理的解。方法是X = B * X+B；x是存储x1至x19的列向量，b，b是前述矩阵，x10、x11、x14、x15、x16、x18和x19依次构成列向量x。所谓的x10，x11，x14，x15，x16，x18，x19都是用nchoosek()组合函数和perms()置换函数。结合循环。关键是要搞清楚以上两个功能的含义

CombinationX = nchosek(1:19，7)；构造1~19中七的所有组合，并将结果以行向量的形式形成矩阵，注意行向量！

[rowCountsCom，column counts com]= size(Combinationx)；获取组合结果矩阵的行数。rowCountsCom表示所有组合数，后面跟所有排列数，也是这个道理。

Fori=1:rowCountsCom遍历每个组合，直到完成所有行，即所有组合。

PermutationX = perms(Combinationx(I，))；取当前第I个组合，计算所有排列结果。相同的结果以行向量的形式存储在结果矩阵中

[行计数精子，列计数精子]=大小(PermutationX)；

Forj = 1: rowcountspell也遍历所有排列

x=PermutationX(j，)；采取其中一种安排

x = B * x '+B；计算满解x1~x19，注意这里用的是x’，因为上面的x是行向量，应该转置成列向量。

xs sorted = round(sort(X))；计算结果应按照升序排序，四舍五入应避免计算过程中出现微小误差，导致程序误判结果

如果(isequal(template，XSorted))结果相等，

点击=点击+1找到一个，添加一个

XTable(:，hits)=X最终结果被添加到解决方案的记录中，并存储为列向量。

结束

计数=计数+1；为每个内部循环添加一个计数。

fprintf('总计253955520次尝试。现在%d...点击%dn '，计数，点击)；这一行是为了避免MATLAB运行的枯燥感，让命令窗口不断显示统计周期次数，相当于界面程序的进度条。。当然缺点是会大大降低程序的运行速度。所以不需要的可以注释掉。

结束

结束

toc计时结束。Tic开始，toc结束，非常生动。按表计时tic，按toc结束计时。命令窗口显示程序运行的总时间。如前所述，这个程序运行时间很长，因为是穷举法，需要循环2亿次以上。下面是运行界面，还没有完成。到时候可以自己在XTable里查结果。图9的结果是一组到位的运行解，其他解只是位置变换。运行程序10分钟左右就会出来两套解决方案，也就是此时显示Hits=2。如果此时不想等所有结果出来，可以按ctrl+C停止MATLAB运行，然后查看XTable结果。

-

我们会在这里等你