六边形幻方

wayne

20# hujunhua
的确是7个参数，我对增广矩阵做行变换，发现 只有12行，也就是说有7个独立参数，代码如下：

1. b1={{1,2,3},{4,5,6,7},{8,9,10,11,12},{13,14,15,16},{17,18,19}};
   
2. b2={{1,4,8},{2,5,9,13},{3,6,10,14,17},{7,11,15,18},{12,16,19}};
   
3. b3={{3,7,12},{2,6,11,16},{1,5,10,15,19},{4,9,14,18},{8,13,17}};
   
4. cc=Append[ReplacePart[PadLeft[{38},20],List/@#->1]&/@Join[b1,b2,b3],PadLeft[{5*38},20,1]];
   
5. RowReduce[cc]
   

复制代码

wayne

即便是7个参数，也只有C(19,7) = 50388 种情况，这个，枚举一下也是非常快的
@chyanog，写一个程序呗

这7个独立参数可以从上面约简后的矩阵中提取

wayne

20# hujunhua
看到老大把7改成8了，这是要考虑一般情形啊

wayne

有一种最懒的方法，预先设定任意两个数的值(总共有19*18种)，然后解线性方程的整数解，

很不幸，把机子扔在一边穷搜大概五分钟的时间，得到的12个答案其实是一个：
{3,17,18,19,7,1,11,16,2,5,6,9,12,4,8,14,10,13,15}
{3,19,16,17,7,2,12,18,1,5,4,10,11,6,8,13,9,14,15}
{9,11,18,14,6,1,17,15,8,5,7,3,13,4,2,19,10,12,16}
{9,14,15,11,6,8,13,18,1,5,4,10,17,7,2,12,3,19,16}
{10,12,16,13,4,2,19,15,8,5,7,3,14,6,1,17,9,11,18}
{10,13,15,12,4,8,14,16,2,5,6,9,19,7,1,11,3,17,18}
{15,13,10,14,8,4,12,9,6,5,2,16,11,1,7,19,18,17,3}
{15,14,9,13,8,6,11,10,4,5,1,18,12,2,7,17,16,19,3}
{16,12,10,19,2,4,13,3,7,5,8,15,17,1,6,14,18,11,9}
{16,19,3,12,2,7,17,10,4,5,1,18,13,8,6,11,15,14,9}
{18,11,9,17,1,6,14,3,7,5,8,15,19,2,4,13,16,12,10}
{18,17,3,11,1,7,19,9,6,5,2,16,14,8,4,12,15,13,10}

chyanog

22# wayne
用Mathematica我还写不出快一点的代码，还有计算Permutations[Range[19], {7}]内存不够

hujunhua

容易回避Permutations[Range[19], {7}]的计算。方法如下：
1、容易证明中心数不大于8,  将中心数取1~8分别计算
2、内圈6数中最大的小于16，从中取5个数，只需要计算Permutations[Range[16], {5}]
3、角格6数中最大的大于15，取一个角格数从16~19分别计算.

hujunhua

除12个对称解，通过编程搜索，还找到了7个非对称解，包括2个零中心的和5个非零中心的. 12个对称解如下图：

对称解的内圈6格与对应的6个角格可以互换，所以12个对称解可以分为6对，第1列-第2列，第3列-第4列。
第一行的4个图，将每个数乘以3后模19取绝对值最小剩余（正负对称），可以得到第二行的对应图。真是奇异的变换关系。

wayne

27# hujunhua
不会吧，我已经搜完了，好想只有12个对称解。

hujunhua

2个零中心非对称解

5个非零中心解

hujunhua

我选择的7个参数格是旋转对称的，即下图中的 a,b,c,d,e,f,o. 其它各格都可用这7个参数表达。


编程时，这7个参数分成3组：
1、中心格取非负数0~9，因为当中心格为负数时可以全体取相反数，将中心格转为正数。
2、{a,b,c}=x, 记除去中心数o后的18个数表为tab, x扫描从tab中取3个不同数的组合, 组合数为C(18,3)=816.
3、{d,e,f}=y, 记从tab中除去x的剩余集为tab1,  y扫描从tab1中取3个不同数的排列，排列数为P(15,3)=2730.
这样就避免了产生P(19,7)的大表。

Mathematica程序如下：

1. term[x_List,y_List,z_]:=Module[{a,b,c,d,e,f,o=z},{a,b,c}=x;{d,e,f}=y;(a+b+c-o-d-e-f)/3];
   
2. Hexagon[x_List,y_List,z_]:=Module[{a,b,c,d,e,f,o=z},{a,b,c}=x;{d,e,f}=y;s=(a+b+c-o-d-e-f)/3;t=-s-o;{e,c-e-s,-c+s,a-e-s,b,-f+t,c-d-s,-a+s,-d+t,a,d,a-f-s,c,-e+t,b-d-s,f,b-f-s,-b+s}];
   
3. Hexagono[x_List,y_List,z_]:=Module[{a,b,c,d,e,f,o=z},{a,b,c}=x;{d,e,f}=y;s=(a+b+c-o-d-e-f)/3;t=-s-o;{{e,c-e-s,-c+s},{a-e-s,b,-f+t,c-d-s},{-a+s,-d+t,o,a,d},{a-f-s,c,-e+t,b-d-s},{f,b-f-s,-b+s}}];
   
4. Do[tab=Cases[Range[-9,9],Except[z]];
   
5.   Do[tab1=Complement[tab,x] ;
   
6.     Do[If[IntegerQ[term[x,y,z]&&Union@Hexagon[x,y,z]==tab],Print[MatrixForm@Hexagono[x,y,z]]],
   
7.     {y,Permutations[tab1,{3}]}],
   
8.   {x,Subsets[tab,{3}]}],
   
9. {z,0,9}]

复制代码

不过实际上我是把z=0(零中心)的情况单独计算的，所以上述程序中{z,0,9}是写成{z,9}的。