《射鵰英雄傳》裡面有一個情節，郭靖帶著受傷的黃蓉四處求高人療傷，遇見瑛姑。瑛姑也愛好各種奇門術數，但是花了好多年卻解不出一個三階幻方。這個三階幻方也就是“洛書”，它有三行三列，九個空格分別填上一到九這九個數字，使得每行、每列、每條對角線上三個數的和都相等。

黃蓉是黃老邪的女兒，古靈精怪，自然也精通此道，很快告訴了瑛姑答案。於是瑛姑就告訴郭靖黃蓉可以找段皇爺療傷。當然瑛姑其實也是為了找段皇爺尋仇。這是後話。其實三階的幻方太簡單了。我倒覺得金庸可以可以把這一段情節改成瑛姑解的是一個五階幻方，就是五行五列的幻方，甚至是七行七列的幻方，這樣花十幾年解不出也情有可原，難度大些，瑛姑也不至於顯得那麼笨。不知道金庸是不是為了襯托黃蓉的聰明？

瑛姑只是在幻方上花了十幾年而已。接下來主要要討論一個六角幻方，也稱幻六邊形，卻花費了一個人前後52年的時間！用盡一生的心血就為解一個幻方。無論如何這種精神很讓人欽佩，畢竟那是一個計算機還沒有普及的年代。

 

一百年前的1910年，一位叫阿當斯的青年人，對六角幻方產生了濃厚興趣。他趁著在鐵路公司閱覽室當職員之便，利用一些空閒時間，去擺弄從1到19這19個數。冬去春來，度過了漫長的47個年頭。經過了無數次的挫折、失敗，使他由一個英俊少年，變成了白髮蒼蒼的老頭，但是他仍然不甘心失敗，這就是興趣的魔力。

1957年的一天，病中的阿當斯，在病床上無意中將六角幻方排列成功了。他驚喜萬分，連忙找紙記錄下來，了卻了他多年的宿願。幾天後，他病癒出院。到家後卻不幸地發現，他填的寶圖不見了。

真是好事多磨，可是阿當斯沒有灰心，他又繼續奮鬥了5年，終於在1962年12月的一天，有志者事竟成，阿當斯又重新填出了他盼望已久的寶圖。

下面就是這個耗費了他52年心血的來之不易的六角幻方。

阿當斯隨即將他的寶圖拿給當時美國的幻方專家馬丁·加德納鑑定。面對這無與倫比的珍奇寶圖，馬丁博士欣喜萬分，當即寫信給才華橫溢的數學遊戲專家特里格。

特里格手捧寶圖敬佩不已。這位專家也一頭扎進了六角幻方，想在層數上作出突破。又耗費了不知多少心血，他才驚奇地發現，兩層以上的六角幻方根本不存在。

1969年，滑鐵盧大學二年級學生阿萊爾，對特里格的結論做出了嚴格的證明，並排除了一些不可能的情況，加上一些條件，利用電子計算機進行測試。僅用了

17秒的時間，就得出了與阿當斯完全相同的結果。電子計算機向人類宣告：雖然普通幻方有千萬種排法，但是，六角幻方卻只有這一個，難怪阿當斯為之奮鬥了52年。

當然阿萊爾用的方法是經過了分析，排除了絕大多數不可能的情況，所以在上個世紀六十年代計算機運算速度還不夠快的情況下，只花17秒的時間就完成了求解。如果不想動腦筋，用暴力窮舉法，這個解的嘗試需要19的階乘這麼多次的嘗試，這個數量是巨大的，也就是19*18*17*16*........*1這麼多次嘗試，就是普通的個人電腦來計算也需要很多年的時間吧。

我不想研究圖形幾何結構，於是想用最笨的窮舉法進行嘗試，畢竟這是計算機的強項。於是我設計了個方法，嘗試的次數只需要19*18*17*16*15*14*13約等於兩億五千萬次，在我電腦上大概運行了5個多小時。估計巨型計算機，只需要秒級的計算時間。我得到12個解，這12個解本質其實就是一個解。因為這一個解可以對稱反轉和旋轉，得到12個位置本質一樣的解。

這個六角幻方可以作為線性代數課程老師的教學案例，完全可以激發學生的挑戰欲和學習興趣，因為它的求解和矩陣的秩，初等變換，線性方程組的解的理論，線性空間等等有重大關係。只要解這麼一道題，很多線性代數的內容就會自然搞懂了。我想如果線性代數老師是這麼教學的，很多同學對抽象的線代至少也會多那麼一些興趣。

這裡附帶了一些六角幻方的問題。有興趣的朋友可以嘗試求解。明天我把我的解答貼出，也當做一個備份。

六角幻方的相關求解及程式碼如下，做個備份。程式採用窮舉法，並不考慮空間結構的約束條件，進行了兩億五千多萬次嘗試（迴圈），實際耗時六個多小時。得到12個解，本質上是一種解的12種不同形式。同時下面還有其他各種問題都挺好的。

下圖的每一列就是一個解。

解答：

1. 需要19個未知數15個方程。

本題的解法就是，給每個小格子設定未知數，從x1~x19，建立方程。如下圖所示

圖1

圖2

AX=0，即為齊次線性方程組。

下圖就代表A矩陣，裡面的資料填的地方都是1，不填的都為0.

圖3

  2. 計算前有4個空格可以填數。單從方程數量和未知數個數對比來看，未知數19個大於方程個數15個，若未經過求解，我們肯定會想，即使係數矩陣A的秩為15（最多也只是15），那麼也多了4個未知數，這四個未知數作為自由引數，可以取任意值。從而對應六邊形圖中的某4個位置，可以任意填數。當然，這四個位置不固定，但也不是隨便都可以（比如，六邊形最外面的六條邊，就不能隨便填數，那樣加起來不為0）

   3. 計算後有7個位置可以任意填數。當把A矩陣輸入MATLAB，計算其秩rank（A）=12，這說明A矩陣真正線性無關的就12行，一定可以化到最簡，解是7個向量的線性組合。7個係數即為7個對應的未知數。（詳細請查閱線性代數中齊次方程組內容）具體如下：

圖4

如上圖，利用MATLAB rref（A）函式對  矩陣A進行 行化簡， 其實也就是相當於高斯迭代法求解AX=0 這一齊次線性方程組，圖中可以驗證rank（A）=12，從圖中看出，選擇x10，x11，x14，x15，x16，x18，x19.這七個變數作為可任填的自由係數。

即方程組改寫為更容易的理解的形式如下，並把上述七個變數新增進去，得到19個解的通解形式：

 

圖5

解用向量的形式表示如下：

圖6

上面的x10，x11，x14，x15，x16，x18，x19.這七個變數就是可以任意填的數，同時它們給出了位置資訊。按照我前面給的x1~x19對應的位置資訊，上面七個位置就是下面圖中的圈出的紅色六邊形的位置，包含那七個變數。你可以先任意填那七個數，然後，其他位置的都可以在圖中利用 和為0直接都手算出來。另外也可以用上圖的通解代入資料驗算，兩種方法結果一致！矩陣的不同化簡方式會導致不同的自由變數選取，這會引起任意填的數的位置的改變，但是一定是像如圖那樣的位於角上的六邊形位置。就好像把下面圖不斷旋轉60，位置就改變不斷改變，有六個這樣的位置可以任意填數。

圖7

4. 數字3和零空間的關係？首先說明下零空間，就是指滿足AX=0這個齊次方程組的x的全體子空間。顯然該零空間的維數為7（包含7個任意的的自由係數，所以維數為7）。那麼3和7的關係就是3=7-4，4是第二問中未經計算的可任填的空格數。

即 零空間維數（7）-（19-15）=3.

5. 顯然是要考線性代數的理論，用的是暴力窮舉法，所以最樸素的想法就是嘗試19！這麼多的次數。現在其實通過上述分析就知道，rank（A）=12，當現在每行每列的和為38時，我們一樣將此時非齊次方程組AX=b的通解求出。 b為全是38的列向量，長度為15. 通過行化簡（利用MATLAB rref函式計算），可得如下圖8所示的通解形式。改變x10，x11，x14，x15，x16，x18，x19這些值，從1~19選取，互不相同，通過下面通解式子計算得到x1~x19，看看這十九個結果是否全部落在1~19，互不重複，是的話，滿足要求，同時得到了一種填法！

需要嘗試A(19,7)=19*18*17*16*15*14*13 這麼多次，效率已經大大提高！

圖8

如下圖9所示就是程式計算出來的一組解，把它旋轉和對稱，有12種情況。本質就是一種解而已。以上所有均在MATLAB程式中實現。

 圖9

現在解釋MATLAB程式碼，MATLAB檔案另附，檔案中不含下面註釋，下面紅色部分為註釋，特此說明：

以下是和為0的計算

a=zeros(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 1014 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 1015 19])=1;

a(15,[2 6 1116])=1;

a(14,[3 712])=1;  以上實現A矩陣，大小15*19。注意MATLAB變數區分大小寫，我這裡先用小寫a表示A矩陣

rank(a);   計算A矩陣的秩，執行完得到A的秩為12

A=rref(a);  對A矩陣進行行化簡，得到如下：

圖10

aEx=[aones(15,1)*38]; AX=b 非齊次方程組對應的增廣矩陣

AEx=rref(aEx);  同上，對增廣矩陣進行行變換，從執行結果可知，

rank(A) =rank(AEx)，所以非齊次方程組有解。AEx矩陣的內容就是A和下面的b，把b拼在A後面就得到AEx

 B=[-1 -1 -1 -2 -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 0 38 38 0 0 0 38 0 0]';

b是從上面計算所得的AEx變數的最右邊一列。因為我寫程式碼是按順序寫的，前面先運行了，看AEx結果了，然後直接輸入這裡的。

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

x=[3 3 7 4-10 -12 21]'; 這個x就是x10，x11，x14，x15，x16，x18，x19這七個任填的資料構成的列向量。這裡給的兩個x的驗證資料就是我前面圖7驗證一驗證二圖中的資料。

XZero=B*x; 這裡就是要驗證“用任何軟體MATLAB，SYMPY，JAVA實現填入實數使得每行或列都為0 填入的數字可以重複。” B*x 實現了圖6的計算，結果儲存在XZero變數中，一下子就得到了全部要填的數。。。可以更改x的值試試。

下面是和為38的計算。

temp=1:19;

template=round(temp');根據前面我的設計思路，就是要利用窮舉法計算出滿足非齊次線性方程組一種可能結果，如果是正確的結果的話，解一定是從1到19的整數，當然順序是可以打亂的，所以到時候把計算的結果排序後同這裡的template做對比，template的值是從1到19的升序的整數。

tic   開始計時，按照前面的解答，需要嘗試A(19,7)=19*18*17*16*15*14*13 這麼多次，約是兩億五千萬次迴圈，所以這裡蠻計時一下，大概需要5~10個小時。視電腦配置和程式執行方式不同而不同。

hits=0          檢測一個結果，加1

Counts=0;        迴圈次數的統計累加

XTable=zeros(19,20);  符合要求的填數的結果變數。每個結果以列向量的形式存放，蠻假設有20個結果，弄個20列。。其實本質的解只有1個，總共有12種形式。道理很簡單，你想象做好了一個解，你可以把整張直面翻過來，填數方式就變了，這叫做對稱。也可以不斷旋轉60度，所以結合起來共12種解。

下面其實最關鍵的就是針對x10，x11，x14，x15，x16，x18，x19的所有取值情況（任意填的1~19的不重複整數）利用圖8計算出結果，並加以判斷合理的解。方法就是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= nchoosek(1:19,7);構造1~19 裡面任選七的所有組合，把結果以行向量的形式構成矩陣，注意是行向量!

[rowCountsCom,columnCountsCom]=size(CombinationX);取得組合結果矩陣的行數，rowCountsCom代表的是所有組合數，後面有個所有排列數，道理一樣。

for i=1:rowCountsCom 遍歷每一種組合，直到完成所有行數，即所有組合數，

    PermutationX=perms(CombinationX(i,:));取當前第i個組合，計算其所有排列結果，同樣結果以行向量形式存放在結果矩陣中

   [rowCountsPerm,columnCountsPerm]=size(PermutationX);

for j=1:rowCountsPerm 同樣遍歷所有的排列情況

x=PermutationX(j,:);取其中的一種排列

  X=B*x'+b; 計算全解x1~x19，注意這裡使用x’，因為上面x是行向量，要轉置為列向量。

            XSorted=round(sort(X));

相關推薦