现在要求将$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$ 原来还定义了三个对角线方向，每个方向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$
由于约束条件太多了，存在解的可能性很小 (%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)
好像还有12个自由变量 应该无解，代码效率不是很高，但是还是可以穷举完毕
#include <stdio.h>
#include <math.h>
#include <algorithm>
using namespace std;
double v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,v17,v18;
double v19,v20,v21,v22,v23,v24,v25,v26,v27,v28,v29,v30,v31,v32,v33,v34,v35;
double r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13;
void setv()
{
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;
}
int used;
#define R(x)r##x=a;
#define V(x)if(fabs(round(v##x)-v##x)>0.1)continue;\
c=(int)round(v##x);if(c<0||c>=36)continue;\
if(used)continue;\
used=1;
#define P(x)printf("v%d=%d\n",x,(int)round(v##x));
void output()
{
printf("Find one:\n");
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)
P(13)P(14)P(15)P(16)P(17)P(18)P(19)P(20)P(21)P(22)P(23)P(24)
P(25)P(26)P(27)P(28)P(29)P(30)P(31)P(32)P(33)P(34)P(35)
}
int main()
{
int a;
int i,c;
for(i=0;i<12;i++)a=i;
do{
for(i=0;i<36;i++)used=0;used=1;
R(2)R(3)R(4)R(5)R(6)R(7)R(8)R(9)R(10)R(11)R(12)R(13)
setv();
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)
V(13)V(14)V(15)V(16)V(17)V(18)V(19)V(20)V(21)V(22)V(23)V(24)
V(25)V(26)V(27)V(28)V(29)V(30)V(31)V(32)V(33)V(34)V(35)
output();
}while(next_permutation(a,a+12));
} 上面的代码错了，因为其中对参数r2,r3,....,r13代码中只让它们选择0~11之间的值了。
重新修改并运行代码，花费了将近4小时才结束，还是没有找到合法的解。
另外发现我说的方法hujunhua在5#就已经使用了。 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.
一个数据松散的6阶模幻方实例
既然用±{1, 2, ..., 18}也填不成，那咱把填入的数组再放松一些，好歹填一个出来再说。以下是一个松散的实例，36个数总和为零，负向最大数为-109，正的最大数为110，宽度219. 数组也不是正负对称的。
图中红色为负数，excel自动变色的。
可以用此图验证前述各导出律。
这是我取-24至24间的12个随机整数做自由变量，1000次试验所得到的最小宽度。
肯定可以更紧凑啦。 改进程序(加循环）后重试10万多次，运气还不错，猫到一个-42~40的。
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 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