幻方的故事与求解
《射雕英雄傳》里面有一個情節，郭靖帶著受傷的黃蓉四處求高人療傷，遇見瑛姑。瑛姑也愛好各種奇門術數，但是花了好多年卻解不出一個三階幻方。這個三階幻方也就是“洛書”，它有三行三列，九個空格分別填上一到九這九個數字，使得每行、每列、每條對角線上三個數的和都相等。
?
黃蓉是黃老邪的女兒，古靈精怪，自然也精通此道，很快告訴了瑛姑答案。于是瑛姑就告訴郭靖黃蓉可以找段皇爺療傷。當然瑛姑其實也是為了找段皇爺尋仇。這是后話。其實三階的幻方太簡單了。我倒覺得金庸可以可以把這一段情節改成瑛姑解的是一個五階幻方，就是五行五列的幻方，甚至是七行七列的幻方，這樣花十幾年解不出也情有可原，難度大些，瑛姑也不至于顯得那么笨。不知道金庸是不是為了襯托黃蓉的聰明？
瑛姑只是在幻方上花了十幾年而已。接下來主要要討論一個六角幻方，也稱幻六邊形，卻花費了一個人前后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！這么多的次數?，F在其實通過上述分析就知道，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 填入的數字可以重復?！?/span> B*x 實現了圖6的計算，結果保存在XZero變量中，一下子就得到了全部要填的數。。?？梢愿?/span>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));算出來的結果要升序排序，同時round四舍五入避免因為計算過程導致的微小誤差，引起程序誤判結果不等
???????????if?(isequal(template,XSorted))結果相等的話
??????????????? hits=hits+1??????? 找中了一個，加一
??????????????? XTable(:,hits)=X?最終結果添加到解的記錄中，以列向量存放。
??????????? end
??????????? Counts=Counts+1;每一個內層循環加一計數。
??????????? fprintf('Total 253955520 tries. Now %d ...Hits %d\n',Counts,hits);這一行是為了避免MATLAB運行感覺太枯燥，讓command window不斷顯示出統計的循環次數，相當于界面程序的進度條。。當然缺點就是會大大降低程序運行速度。因此不想要可以注釋掉。
??????? end????????
end
toc?? 計時結束。Tic開始，toc結束，很形象吧，按下表計時tic，按下toc結束計時。command window 會顯示程序運行的總時間。按前面說的，本程序運行時間長，因為是窮舉法，要進行兩億多次循環。下面是運行界面，還沒運行完。到時候你可以自己去XTable里面查看結果。圖9的結果就是運行出來的一組位置上的解，其他解就是位置變換而已。運行程序大概10分鐘左右就會出來兩組解，也就是此時顯示Hits=2，這時候要是不想等所有的結果出來，可以按 ctrl+C終止MATLAB運行，然后查看XTable結果。
來源：黃島主的世界
--------------------
明明共同關注公眾號，彼此卻互不認識；
明明具有相同的愛好，卻無緣相識；
有沒有覺得這就是上帝給我們的一個bug！
想不想認識更多寫程序的小伙伴？
C++，Java，VB……應有盡有。
還等什么？趕快上車加入我們吧！
(･??･?)っ算法與數學之美-計算機粉絲群
我們在這里等你喲
算法數學之美微信公眾號歡迎賜稿
稿件涉及數學、物理、算法、計算機、編程等相關領域。
稿件一經采用，我們將奉上稿酬。
投稿郵箱：math_alg@163.com
總結
- 上一篇： mysql命令行导入url_Mysql
- 下一篇： MSN-LDL论文修改（B-Y Rong