6 阶模幻方的搜索

mathe

应该是指，总共有36个位置，分别可以用$a_{i,j}$表示，其中$0\le i\lt 6, 0\le j\lt 6$.
现在要求将$0~17,19~36$填入这36个位置，使得
$R_i=\sum_{j=0}^5 a_{i,j}=18\times 6$
$C_j=\sum_{i=0}^5 a_{i,j}=18\times 6$
$D_k=\sum_{i=0}^5 a_{i, (i+k)(\mod 6)}=18\times 6$

mathe

原来还定义了三个对角线方向，每个方向6条，所以还需要添加额外3*6-3个关系式。
也就是36个变量，6*6-5=31个约束，只余下5个自由变量了。
先解方程得到31个非自由变量的表达式，然后依次试验5个自由变量的取值组合，判断余下的31个导出变量是否还在范围之内。
而对角线可以使用如下坐标定义
$H_{2,j}=\sum_{i=0}^5 a_{i,(2i+j)(\mod 6)}=18\times6$
$H_{5,j}=\sum_{i=0}^5 a_{i,(j-i)(\mod 6)}=18\times6$
$H_{4,j}=\sum_{i=0}^5 a_{(2i+j)(\mod 6),i}=18\times6$
由于约束条件太多了，存在解的可能性很小

mathe

(%i4)        linsolve([+v0+v1+v2+v3+v4+v5=108,+v0+v6+v12+v18+v24+v30=108,+v0+v7+v14+v21+v28+v35=108,+v0+v8+v16+v18+v26+v34=108,+v0+v11+v16+v21+v26+v31=108,+v0+v13+v26+v3+v16+v29=108,+v6+v7+v8+v9+v10+v11=108,+v1+v7+v13+v19+v25+v31=108,+v1+v8+v15+v22+v29+v30=108,+v1+v9+v17+v19+v27+v35=108,+v1+v6+v17+v22+v27+v32=108,+v6+v19+v32+v9+v22+v35=108,+v12+v13+v14+v15+v16+v17=108,+v2+v8+v14+v20+v26+v32=108,+v2+v9+v16+v23+v24+v31=108,+v2+v10+v12+v20+v28+v30=108,+v2+v7+v12+v23+v28+v33=108,+v12+v25+v2+v15+v28+v5=108,+v18+v19+v20+v21+v22+v23=108,+v3+v9+v15+v21+v27+v33=108,+v3+v10+v17+v18+v25+v32=108,+v3+v11+v13+v21+v29+v31=108,+v3+v8+v13+v18+v29+v34=108,+v18+v31+v8+v21+v34+v11=108,+v24+v25+v26+v27+v28+v29=108,+v4+v10+v16+v22+v28+v34=108,+v4+v11+v12+v19+v26+v33=108,+v4+v6+v14+v22+v24+v32=108,+v4+v9+v14+v19+v24+v35=108,+v24+v1+v14+v27+v4+v17=108,+v30+v31+v32+v33+v34+v35=108,+v5+v11+v17+v23+v29+v35=108,+v5+v6+v13+v20+v27+v34=108,+v5+v7+v15+v23+v25+v33=108,+v5+v10+v15+v20+v25+v30=108,+v30+v7+v20+v33+v10+v23=108],[v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16
                ,v17,v18,v19,v20,v21,v22,v23,v24,v25,v26,v27,v28,v29,v30,v31,v32,v33,v34,v35]);
       
solve: dependent equations eliminated: (34 35 28 11 31 32 5 6 2 13 27 21)

(%o4)        [v0=-%r9+%r8-%r7+%r6+%r4-2*%r3+2*%r13+%r12-%r11+%r10-18,v1=(%r9-%r8+2*%r7-2*%r6-%r5-%r4+%r3+%r2-2*%r13-%r12-%r10+108)/2,v2=%r5,v3=-%r9-2*%r7+%r6-%r5-%r4-%r3-%r2+%r13-%r11+126,v4=%r8,v5=(3*%r9-3*%r8+4*%r7-2*%r6+%r5+%r4+5*%r3+%r2-4*%r13-%r12+4*%r11-%r10-108)/2,v6=-%r4-%r13-%r10+72,v7=%r3,v8=%r11,v9=%r4,v10=%r7,v11=-%r7-%r3+%r13-%r11+%r10+36,v12=-%r9-%r5+54,v13=-(-%r9-%r8-%r5-%r4+%r3+%r2-2*%r13+%r12-%r10+36)/2,v14=-%r8-%r13+54,v15=-(%r9-3*%r8+2*%r7-2*%r6+%r5+%r4+3*%r3+%r2-4*%r13-%r12+4*%r11-%r10-72)/2,v16=(3*%r9-3*%r8+4*%r7-4*%r6+%r5-%r4+5*%r3+3*%r2-6*%r13-%r12+4*%r11-%r10-36)/2,v17=-(%r9-%r8+2*%r7-2*%r6-%r5-%r4+%r3+%r2-2*%r13-%r12+%r10)/2,v18=2*%r9-%r8+2*%r7-%r6+%r5+2*%r3+%r2-2*%r13-%r12+%r11-54,v19=%r12,v20=%r2,v21=%r7-%r6+%r3-%r13+%r11-%r10+18,v22=-(%r9+%r8+2*%r7-2*%r6-%r5-%r4+%r3+%r2-2*%r13+%r12-%r10-36)/2,v23=-(3*%r9-3*%r8+4*%r7-2*%r6+3*%r5+%r4+5*%r3+3*%r2-4*%r13-%r12+4*%r11-%r10-252)/2,v24=%r13,v25=-%r9-%r7-%r3+72,v26=-(%r9-%r8+2*%r7-2*%r6+%r5+%r4+%r3+3*%r2-2*%r13+%r12+2*%r11+%r10-180)/2,v27=%r10,v28=%r9,v29=(%r9-%r8+4*%r7-2*%r6+%r5+%r4+3*%r3+3*%r2-4*%r13+%r12+2*%r11-%r10-108)/2,v30=-%r7-%r2+54,v31=%r6,v32=(%r9+%r8+2*%r7-2*%r6-%r5+%r4+%r3+%r2+%r12+%r10-72)/2,v33=(3*%r9-3*%r8+4*%r7-2*%r6+3*%r5+%r4+3*%r3+3*%r2-4*%r13-%r12+4*%r11-%r10-144)/2,v34=-2*%r9+%r8-2*%r7+%r6-%r5-2*%r3-%r2+2*%r13+%r12-2*%r11+108,v35=-%r4-%r12+54]
好像还有12个自由变量

mathe

应该无解，代码效率不是很高，但是还是可以穷举完毕

1. #include <stdio.h>
   
2. #include <math.h>
   
3. #include <algorithm>
   
4. using namespace std;
   
5. 
   
6. double v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,v17,v18;
   
7. double v19,v20,v21,v22,v23,v24,v25,v26,v27,v28,v29,v30,v31,v32,v33,v34,v35;
   
8. double r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13;
   
9. 
   
10. void setv()
    
11. {
    
12. v0=-r9+r8-r7+r6+r4-2*r3+2*r13+r12-r11+r10-18;
    
13. v1=(r9-r8+2*r7-2*r6-r5-r4+r3+r2-2*r13-r12-r10+108)/2;
    
14. v2=r5;
    
15. v3=-r9-2*r7+r6-r5-r4-r3-r2+r13-r11+126;
    
16. v4=r8;
    
17. v5=(3*r9-3*r8+4*r7-2*r6+r5+r4+5*r3+r2-4*r13-r12+4*r11-r10-108)/2;
    
18. v6=-r4-r13-r10+72;
    
19. v7=r3;
    
20. v8=r11;
    
21. v9=r4;
    
22. v10=r7;
    
23. v11=-r7-r3+r13-r11+r10+36;
    
24. v12=-r9-r5+54;
    
25. v13=-(-r9-r8-r5-r4+r3+r2-2*r13+r12-r10+36)/2;
    
26. v14=-r8-r13+54;
    
27. v15=-(r9-3*r8+2*r7-2*r6+r5+r4+3*r3+r2-4*r13-r12+4*r11-r10-72)/2;
    
28. v16=(3*r9-3*r8+4*r7-4*r6+r5-r4+5*r3+3*r2-6*r13-r12+4*r11-r10-36)/2;
    
29. v17=-(r9-r8+2*r7-2*r6-r5-r4+r3+r2-2*r13-r12+r10)/2;
    
30. v18=2*r9-r8+2*r7-r6+r5+2*r3+r2-2*r13-r12+r11-54;
    
31. v19=r12;
    
32. v20=r2;
    
33. v21=r7-r6+r3-r13+r11-r10+18;
    
34. v22=-(r9+r8+2*r7-2*r6-r5-r4+r3+r2-2*r13+r12-r10-36)/2;
    
35. v23=-(3*r9-3*r8+4*r7-2*r6+3*r5+r4+5*r3+3*r2-4*r13-r12+4*r11-r10-252)/2;
    
36. v24=r13;
    
37. v25=-r9-r7-r3+72;
    
38. v26=-(r9-r8+2*r7-2*r6+r5+r4+r3+3*r2-2*r13+r12+2*r11+r10-180)/2;
    
39. v27=r10;
    
40. v28=r9;
    
41. v29=(r9-r8+4*r7-2*r6+r5+r4+3*r3+3*r2-4*r13+r12+2*r11-r10-108)/2;
    
42. v30=-r7-r2+54;
    
43. v31=r6;
    
44. v32=(r9+r8+2*r7-2*r6-r5+r4+r3+r2+r12+r10-72)/2;
    
45. v33=(3*r9-3*r8+4*r7-2*r6+3*r5+r4+3*r3+3*r2-4*r13-r12+4*r11-r10-144)/2;
    
46. v34=-2*r9+r8-2*r7+r6-r5-2*r3-r2+2*r13+r12-2*r11+108;
    
47. v35=-r4-r12+54;
    
48. }
    
49. 
    
50. int used[36];
    
51. #define R(x)  r##x=a[x-2];
    
52. #define V(x)  if(fabs(round(v##x)-v##x)>0.1)continue;\
    
53.               c=(int)round(v##x);if(c<0||c>=36)continue;\
    
54.               if(used[c])continue;\
    
55.               used[c]=1;
    
56. 
    
57. #define P(x)  printf("v%d=%d\n",x,(int)round(v##x));
    
58. void output()
    
59. {
    
60.     printf("Find one:\n");
    
61.     P(0)P(1)P(2)P(3)P(4)P(5)P(6)P(7)P(8)P(9)P(10)P(11)P(12)
    
62.     P(13)P(14)P(15)P(16)P(17)P(18)P(19)P(20)P(21)P(22)P(23)P(24)
    
63.     P(25)P(26)P(27)P(28)P(29)P(30)P(31)P(32)P(33)P(34)P(35)
    
64. }
    
65. 
    
66. int main()
    
67. {
    
68.     int a[12];
    
69.     int i,c;
    
70.     for(i=0;i<12;i++)a[i]=i;
    
71.     do{
    
72.         for(i=0;i<36;i++)used[i]=0;used[18]=1;
    
73.         R(2)R(3)R(4)R(5)R(6)R(7)R(8)R(9)R(10)R(11)R(12)R(13)
    
74.         setv();
    
75.         V(0)V(1)V(2)V(3)V(4)V(5)V(6)V(7)V(8)V(9)V(10)V(11)V(12)
    
76.         V(13)V(14)V(15)V(16)V(17)V(18)V(19)V(20)V(21)V(22)V(23)V(24)
    
77.         V(25)V(26)V(27)V(28)V(29)V(30)V(31)V(32)V(33)V(34)V(35)
    
78.         output();
    
79.     }while(next_permutation(a,a+12));
    
80. }

复制代码

mathe

上面的代码错了，因为其中对参数r2,r3,....,r13代码中只让它们选择0~11之间的值了。
重新修改并运行代码，花费了将近4小时才结束，还是没有找到合法的解。
另外发现我说的方法hujunhua在5#就已经使用了。

hujunhua

chyanog 发表于 2019-12-15 17:15
这句话怎么理解？

由于没有构造出6阶模幻方的实例，谨以一个7阶模幻方的实例来加以说明。

如图所示为一个7阶模幻方，其中填入了0, 1, 2, ..., 48（用二位七进数00, 01, ..., 66表示）。
取任意一个7X7的菱形视野，均正好包含0~48 。菱形可拉正为正方形，以合其名称中的“方”字。
为了表现其旋转对称性，我们取了一个正六方形视野，将其中49个数涂红。
三条黑实线方向为密排方向，三条橙色虚线方向为次密排方向。
在任一方向的直线上的格位以最小正周期7重复，视为同一格。
即任一方向的一条直线实际上只有7个数。
这个7阶幻方的约束是：
1、任意的连续7个密排数之和为7X24.
2、任意的连续7个次密排数之和为7X24.

这样定义的7阶模幻方的导出律是：
1、密排方向的任意正六边形，六角加中心数之和为7X24.
2、次密排方向的任意正六边形，六角加中心数之和为7X24.

hujunhua

既然用±{1, 2, ..., 18}也填不成，那咱把填入的数组再放松一些，好歹填一个出来再说。
以下是一个松散的实例，36个数总和为零，负向最大数为-109，正的最大数为110，宽度219. 数组也不是正负对称的。
图中红色为负数，excel自动变色的。

可以用此图验证前述各导出律。
这是我取-24至24间的12个随机整数做自由变量，1000次试验所得到的最小宽度。
肯定可以更紧凑啦。

hujunhua

改进程序(加循环）后重试10万多次，运气还不错，猫到一个-42~40的。

mathe

7阶的解好像很多很多，hujunhua帮忙验证一下，我的理解是否正确的，比如:
Find:
        a_(0,0)=6               a_(0,1)=25              a_(0,2)=7               a_(0,3)=30              a_(0,4)=47              a_(0,5)=15              a_(0,6)=38
        a_(1,0)=12              a_(1,1)=29              a_(1,2)=45              a_(1,3)=20              a_(1,4)=39              a_(1,5)=0               a_(1,6)=23
        a_(2,0)=46              a_(2,1)=14              a_(2,2)=37              a_(2,3)=5               a_(2,4)=22              a_(2,5)=10              a_(2,6)=34
        a_(3,0)=36              a_(3,1)=3               a_(3,2)=27              a_(3,3)=11              a_(3,4)=28              a_(3,5)=44              a_(3,6)=19
        a_(4,0)=21              a_(4,1)=9               a_(4,2)=33              a_(4,3)=43              a_(4,4)=17              a_(4,5)=41              a_(4,6)=4
        a_(5,0)=31              a_(5,1)=48              a_(5,2)=18              a_(5,3)=35              a_(5,4)=2               a_(5,5)=26              a_(5,6)=8
        a_(6,0)=16              a_(6,1)=40              a_(6,2)=1               a_(6,3)=24              a_(6,4)=13              a_(6,5)=32              a_(6,6)=42
Find:
        a_(0,0)=6               a_(0,1)=25              a_(0,2)=7               a_(0,3)=30              a_(0,4)=40              a_(0,5)=15              a_(0,6)=45
        a_(1,0)=12              a_(1,1)=29              a_(1,2)=38              a_(1,3)=20              a_(1,4)=46              a_(1,5)=0               a_(1,6)=23
        a_(2,0)=39              a_(2,1)=14              a_(2,2)=44              a_(2,3)=5               a_(2,4)=22              a_(2,5)=10              a_(2,6)=34
        a_(3,0)=43              a_(3,1)=3               a_(3,2)=27              a_(3,3)=11              a_(3,4)=28              a_(3,5)=37              a_(3,6)=19
        a_(4,0)=21              a_(4,1)=9               a_(4,2)=33              a_(4,3)=36              a_(4,4)=17              a_(4,5)=48              a_(4,6)=4
        a_(5,0)=31              a_(5,1)=41              a_(5,2)=18              a_(5,3)=42              a_(5,4)=2               a_(5,5)=26              a_(5,6)=8
        a_(6,0)=16              a_(6,1)=47              a_(6,2)=1               a_(6,3)=24              a_(6,4)=13              a_(6,5)=32              a_(6,6)=35

mathe

6阶放大范围到45：
              a_(0,0)=8                a_(0,1)=24              a_(0,2)=25              a_(0,3)=40              a_(0,4)=15              a_(0,5)=26
              a_(1,0)=14              a_(1,1)=31              a_(1,2)=21              a_(1,3)=23              a_(1,4)=19              a_(1,5)=30
              a_(2,0)=44              a_(2,1)=17              a_(2,2)=38              a_(2,3)=1               a_(2,4)=32              a_(2,5)=4
              a_(3,0)=11              a_(3,1)=18              a_(3,2)=5               a_(3,3)=33              a_(3,4)=35              a_(3,5)=36
              a_(4,0)=16              a_(4,1)=42              a_(4,2)=29              a_(4,3)=39              a_(4,4)=0               a_(4,5)=12
              a_(5,0)=45              a_(5,1)=6                a_(5,2)=20              a_(5,3)=2               a_(5,4)=37              a_(5,5)=28