幻方问题
幻方(Magic Square)是一种将数字安排在正方形格子中，使每行、列和对角线上的数字和都相等的方法。
幻方包括：奇数幻方、单偶数幻方和双偶数幻方
奇数幻方(2n+1类): 如: 3, 5, 7
单偶数幻方(即4n+2类)：幻方阶数n为不能被4整除的偶数。如: 6、10、14
双偶数幻方(即4n类)：幻方行列数n能被4整除的幻方。如:4、8、12
奇数幻方实现
import numpy as np
def oddMagic(rank: int):
if (rank < 3) or (rank % 2) == 0:
print("Rank数据必须大于2,且是奇数")
return None
if rank > 1000:
print("rank数值太大")
return None
curRow = 0
curCol = rank // 2
magic = np.zeros((rank, rank), dtype=np.uint32)
magic[curRow, curCol] = 1
i = 2
endNum = rank * rank
while i <= endNum:
tmpRow = curRow - 1
tmpCol = curCol + 1
# row,col正常数值范围
if (0 <= tmpRow <= rank-1) and (0 <= tmpCol <= rank-1):
# 前面有数字(row-1,col不变)
if magic[tmpRow, tmpCol] > 0:
curRow = curRow + 1
else:
curRow = tmpRow
curCol = tmpCol
magic[curRow, curCol] = i
# row越界
if (tmpRow < 0) and (0 <= tmpCol <= rank-1):
# row-> rank - 1
curRow = rank - 1
curCol = tmpCol
magic[curRow, curCol] = i
# col越界
if (0 <= tmpRow <= rank - 1) and (tmpCol > rank - 1):
curRow = tmpRow
curCol = 0
magic[curRow, curCol] = i
# row,col都越界
if (tmpRow < 0) and (tmpCol > rank - 1):
curRow = curRow + 1
magic[curRow, curCol] = i
i = i + 1
return magic
11阶奇数幻方结果：
sumTotal=7381, perValue=671
validate->: True
[[ 68 81 94 107 120 1 14 27 40 53 66]
[ 80 93 106 119 11 13 26 39 52 65 67]
[ 92 105 118 10 12 25 38 51 64 77 79]
[104 117 9 22 24 37 50 63 76 78 91]
[116 8 21 23 36 49 62 75 88 90 103]
[ 7 20 33 35 48 61 74 87 89 102 115]
[ 19 32 34 47 60 73 86 99 101 114 6]
[ 31 44 46 59 72 85 98 100 113 5 18]
[ 43 45 58 71 84 97 110 112 4 17 30]
[ 55 57 70 83 96 109 111 3 16 29 42]
[ 56 69 82 95 108 121 2 15 28 41 54]]
双偶数幻方实现
def evenMagicDouble(rank: int):
if (rank < 4) or (rank % 4 != 0):
print("rank值必须是4的倍数,且大于4!")
return None
if rank > 1000:
print("rank数值太大")
return None
# 对每个元素按行列依次赋值
# refer: https://hanspub.org/journal/PaperInformation?paperID=64877
magic = np.zeros((rank, rank), dtype=np.int32)
for tmpRow in range(rank):
for tmpCol in range(rank):
magic[tmpRow, tmpCol] = rank * tmpRow + (tmpCol + 1)
print("step1:\n", magic)
# 分成n*n个4x4的小方阵(对角线上的元素标记为0)
maskMagic = np.ones((rank, rank), dtype=np.int32)
nDim = rank // 4
nSpan = 4
for tmpY in range(nDim):
for tmpX in range(nDim):
offsetY = nSpan * tmpY
offsetX = nSpan * tmpX
for index in range(nSpan):
# 主对角线
tmpRow = index + offsetY
tmpCol = index + offsetX
maskMagic[tmpRow, tmpCol] = 0
# 副对象线
tmpRow = nSpan - 1 - index + offsetY
tmpCol = index + offsetX
maskMagic[tmpRow, tmpCol] = 0
#print("step2:\n", maskMagic)
# 将n*n平方个小方阵的对角线元素按(rank*rank方阵的中点)对称交换
# 元素交换的互补值
compleValue = rank * rank + 1
for tmpRow in range(rank):
for tmpCol in range(rank):
if maskMagic[tmpRow, tmpCol] == 0:
magic[tmpRow, tmpCol] = compleValue - magic[tmpRow, tmpCol]
#print("step3:\n", magic)
return magic
16阶双偶数幻方结果：
validate->: True
[[256 2 3 253 252 6 7 249 248 10 11 245 244 14 15 241]
[ 17 239 238 20 21 235 234 24 25 231 230 28 29 227 226 32]
[ 33 223 222 36 37 219 218 40 41 215 214 44 45 211 210 48]
[208 50 51 205 204 54 55 201 200 58 59 197 196 62 63 193]
[192 66 67 189 188 70 71 185 184 74 75 181 180 78 79 177]
[ 81 175 174 84 85 171 170 88 89 167 166 92 93 163 162 96]
[ 97 159 158 100 101 155 154 104 105 151 150 108 109 147 146 112]
[144 114 115 141 140 118 119 137 136 122 123 133 132 126 127 129]
[128 130 131 125 124 134 135 121 120 138 139 117 116 142 143 113]
[145 111 110 148 149 107 106 152 153 103 102 156 157 99 98 160]
[161 95 94 164 165 91 90 168 169 87 86 172 173 83 82 176]
[ 80 178 179 77 76 182 183 73 72 186 187 69 68 190 191 65]
[ 64 194 195 61 60 198 199 57 56 202 203 53 52 206 207 49]
[209 47 46 212 213 43 42 216 217 39 38 220 221 35 34 224]
[225 31 30 228 229 27 26 232 233 23 22 236 237 19 18 240]
[ 16 242 243 13 12 246 247 9 8 250 251 5 4 254 255 1]]
单偶数幻方实现
def evenMagicSingle(rank: int):
if not (rank >= 6 and ((rank % 2 == 0) and (rank % 4 != 0))):
print("rank值必须是2的奇数倍,且大于5!")
return None
if rank > 1000:
print("rank数值太大")
return None
# 罗伯法
divRank = rank // 2
oriMagic = oddMagic(divRank)
offsetValue = divRank * divRank
mags = []
for index in range(4):
tmpMagic = oriMagic + offsetValue * index
mags.append(tmpMagic)
print("step1:\n")
for i in range(len(mags)):
print(f"[{i}]:\n", mags[i])
"""
4个小方阵按ACDB(原来的顺序是ABCD)方式合成一个大的方阵
A B => A C
C D D B
"""
mH1 = np.hstack((mags[0], mags[2]))
mH2 = np.hstack((mags[3], mags[1]))
magic = np.vstack((mH1, mH2))
print("step2:\n", magic)
"""
refer1: https://zhidao.baidu.com/question/1608589126045729347.html
A象限和D象限对换数据(相对应的位置的元素)
在A每行取m个小格（中心格及一侧对角线格为必换格，其余m-1格只要不是另一侧对角线格即可），
简单地说，就是说在A中间一行取包括中心格在内的m个小格，其他行左侧（或右侧）边缘取m个小格，将其与D相应方格内交换；
B与C在任取m-1列相互交换（6阶幻方m=1，m-1=0，B与C列不用相互交换）
"""
mGrid = divRank // 2
# 记录对角线元素的状态值(1)
maskMagic = np.zeros((divRank, divRank))
for index in range(divRank):
maskMagic[index, index] = 1
maskMagic[divRank-1-index, index] = 1
print("maskMagic:\n", maskMagic)
# 每一行需要交换mGrid个元素(一侧对角线格为必换格)
for tmpRow in range(divRank):
# 交换对象线元素
for tmpCol in range(divRank):
if maskMagic[tmpRow, tmpCol] == 1:
# 小方阵的元素行列索引值要相对应
curRowA = tmpRow
curColA = tmpCol
curRowD = curRowA + divRank
curColD = curColA
tmpValue = magic[curRowA, curColA]
magic[curRowA, curColA] = magic[curRowD, curColD]
magic[curRowD, curColD] = tmpValue
break
# 其它非对角线上元素交换
rmnCnt = 0
while rmnCnt < mGrid - 1:
for tmpCol in range(divRank):
if maskMagic[tmpRow, tmpCol] != 1:
curRowA = tmpRow
curColA = tmpCol
curRowD = curRowA + divRank
curColD = curColA
tmpValue = magic[curRowA, curColA]
magic[curRowA, curColA] = magic[curRowD, curColD]
magic[curRowD, curColD] = tmpValue
# 标记已交换的状态
maskMagic[tmpRow, tmpCol] = 1
break
rmnCnt += 1
# B与C再任取m-1列相互交换
tmpColCnt = 0
while tmpColCnt < mGrid - 1:
for tmpRow in range(divRank):
curRowC = tmpRow
curColC = tmpColCnt + divRank
curRowB = curRowC + divRank
curColB = curColC
tmpValue = magic[curRowC, curColC]
magic[curRowC, curColC] = magic[curRowB, curColB]
magic[curRowB, curColB] = tmpValue
tmpColCnt += 1
print("step3:\n", magic)
return magic
14阶单偶数幻方结果：
sumTotal=19306, perValue=1379
validate->: True
[[177 186 195 1 10 19 28 79 88 146 99 108 117 126]
[185 194 154 9 18 27 29 87 96 105 107 116 125 127]
[193 153 155 17 26 35 37 95 55 106 115 124 133 135]
[152 161 16 172 34 36 45 54 63 114 123 132 134 143]
[160 162 171 33 42 44 4 62 64 122 131 140 142 102]
[168 170 179 41 43 3 12 70 72 130 139 141 101 110]
[169 178 187 49 2 11 20 71 80 138 147 100 109 118]
[ 30 39 48 148 157 166 175 128 137 97 50 59 68 77]
[ 38 47 7 156 165 174 176 136 145 56 58 67 76 78]
[ 46 6 8 164 173 182 184 144 104 57 66 75 84 86]
[ 5 14 163 25 181 183 192 103 112 65 74 83 85 94]
[ 13 15 24 180 189 191 151 111 113 73 82 91 93 53]
[ 21 23 32 188 190 150 159 119 121 81 90 92 52 61]
[ 22 31 40 196 149 158 167 120 129 89 98 51 60 69]]
验证幻方是否正确
def validateMagic(magic: np.ndarray, rank: int):
maxNum = rank * rank
sumTotal = (1 + maxNum) * maxNum // 2
perValue = sumTotal // rank
print(f"sumTotal={sumTotal}, perValue={perValue}")
# validate row
for tmpRow in range(rank):
tmpSum = np.sum(magic[tmpRow,0:rank])
if tmpSum != perValue:
print(f"error row:{tmpSum} != {perValue}", magic[tmpRow,0:rank])
return False
# validate col
for tmpCol in range(rank):
tmpSum = np.sum(magic[0:rank, tmpCol])
if tmpSum != perValue:
print(f"error col:{tmpSum} != {perValue}", magic[0:rank, tmpCol])
return False
# validate hypotenuse 1
tmpSum = 0
for index in range(rank):
tmpSum += magic[index,index]
if tmpSum != perValue:
print(f"error hypotenuse1: col:{tmpSum} != {perValue}")
return False
# validate hypotenuse 2
tmpSum = 0
for index in range(rank):
tmpSum += magic[rank-1-index,index]
if tmpSum != perValue:
print(f"error hypotenuse1: col:{tmpSum} != {perValue}")
return False
return True
全部代码
"""
幻方（Magic Square）是一种将数字安排在正方形格子中，使每行、列和对角线上的数字和都相等的方法。
幻方包括：奇数幻方、单偶数幻方和双偶数幻方
单偶数幻方（即4n+2式）：幻方阶数n为不能被4整除的偶数，比如6、10、14
双偶数幻方（即4n式）：幻方行列数n能被4整除的幻方，比如4、8、12
"""
import numpy as np
def oddMagic(rank: int):
if (rank < 3) or (rank % 2) == 0:
print("Rank数据必须大于2,且是奇数")
return None
if rank > 1000:
print("rank数值太大")
return None
curRow = 0
curCol = rank // 2
magic = np.zeros((rank, rank), dtype=np.uint32)
magic[curRow, curCol] = 1
i = 2
endNum = rank * rank
while i <= endNum:
tmpRow = curRow - 1
tmpCol = curCol + 1
# row,col正常数值范围
if (0 <= tmpRow <= rank-1) and (0 <= tmpCol <= rank-1):
# 前面有数字(row-1,col不变)
if magic[tmpRow, tmpCol] > 0:
curRow = curRow + 1
else:
curRow = tmpRow
curCol = tmpCol
magic[curRow, curCol] = i
# row越界
if (tmpRow < 0) and (0 <= tmpCol <= rank-1):
# row-> rank - 1
curRow = rank - 1
curCol = tmpCol
magic[curRow, curCol] = i
# col越界
if (0 <= tmpRow <= rank - 1) and (tmpCol > rank - 1):
curRow = tmpRow
curCol = 0
magic[curRow, curCol] = i
# row,col都越界
if (tmpRow < 0) and (tmpCol > rank - 1):
curRow = curRow + 1
magic[curRow, curCol] = i
i = i + 1
return magic
# 交换方阵的对角元素
def swapHypotenuseEle(magic: np.ndarray):
tmpRank = magic.shape[0]
semiRank = tmpRank // 2
semiHigh = tmpRank - 1
centerPntY = int((semiHigh / 2 + 0) * 2)
centerPntX = int((semiHigh / 2 + 0) * 2)
for index in range(semiRank):
tmpRowA = index
tmpColA = index
tmpRowB = centerPntY - tmpRowA
tmpColB = centerPntX - tmpColA
tmpValue = magic[tmpRowA, tmpColA]
magic[tmpRowA, tmpColA] = magic[tmpRowB, tmpColB]
magic[tmpRowB, tmpColB] = tmpValue
tmpRowA = semiHigh - index
tmpColA = index
tmpRowB = centerPntY - tmpRowA
tmpColB = centerPntX - tmpColA
tmpValue = magic[tmpRowA, tmpColA]
magic[tmpRowA, tmpColA] = magic[tmpRowB, tmpColB]
magic[tmpRowB, tmpColB] = tmpValue
return magic
def evenMagicDouble(rank: int):
if (rank < 4) or (rank % 4 != 0):
print("rank值必须是4的倍数,且大于4!")
return None
if rank > 1000:
print("rank数值太大")
return None
# 对每个元素按行列依次赋值
# refer: https://hanspub.org/journal/PaperInformation?paperID=64877
magic = np.zeros((rank, rank), dtype=np.int32)
for tmpRow in range(rank):
for tmpCol in range(rank):
magic[tmpRow, tmpCol] = rank * tmpRow + (tmpCol + 1)
print("step1:\n", magic)
# 分成n*n个4x4的小方阵(对角线上的元素标记为0)
maskMagic = np.ones((rank, rank), dtype=np.int32)
nDim = rank // 4
nSpan = 4
for tmpY in range(nDim):
for tmpX in range(nDim):
offsetY = nSpan * tmpY
offsetX = nSpan * tmpX
for index in range(nSpan):
# 主对角线
tmpRow = index + offsetY
tmpCol = index + offsetX
maskMagic[tmpRow, tmpCol] = 0
# 副对象线
tmpRow = nSpan - 1 - index + offsetY
tmpCol = index + offsetX
maskMagic[tmpRow, tmpCol] = 0
#print("step2:\n", maskMagic)
# 将n*n平方个小方阵的对角线元素按(rank*rank方阵的中点)对称交换
# 元素交换的互补值
compleValue = rank * rank + 1
for tmpRow in range(rank):
for tmpCol in range(rank):
if maskMagic[tmpRow, tmpCol] == 0:
magic[tmpRow, tmpCol] = compleValue - magic[tmpRow, tmpCol]
#print("step3:\n", magic)
return magic
def evenMagicSingle(rank: int):
if not (rank >= 6 and ((rank % 2 == 0) and (rank % 4 != 0))):
print("rank值必须是2的奇数倍,且大于5!")
return None
if rank > 1000:
print("rank数值太大")
return None
# 罗伯法
divRank = rank // 2
oriMagic = oddMagic(divRank)
offsetValue = divRank * divRank
mags = []
for index in range(4):
tmpMagic = oriMagic + offsetValue * index
mags.append(tmpMagic)
print("step1:\n")
for i in range(len(mags)):
print(f"[{i}]:\n", mags[i])
"""
4个小方阵按ACDB(原来的顺序是ABCD)方式合成一个大的方阵
A B => A C
C D D B
"""
mH1 = np.hstack((mags[0], mags[2]))
mH2 = np.hstack((mags[3], mags[1]))
magic = np.vstack((mH1, mH2))
print("step2:\n", magic)
"""
refer1: https://zhidao.baidu.com/question/1608589126045729347.html
A象限和D象限对换数据(相对应的位置的元素)
在A每行取m个小格（中心格及一侧对角线格为必换格，其余m-1格只要不是另一侧对角线格即可），
简单地说，就是说在A中间一行取包括中心格在内的m个小格，其他行左侧（或右侧）边缘取m个小格，将其与D相应方格内交换；
B与C在任取m-1列相互交换（6阶幻方m=1，m-1=0，B与C列不用相互交换）
"""
mGrid = divRank // 2
# 记录对角线元素的状态值(1)
maskMagic = np.zeros((divRank, divRank))
for index in range(divRank):
maskMagic[index, index] = 1
maskMagic[divRank-1-index, index] = 1
print("maskMagic:\n", maskMagic)
# 每一行需要交换mGrid个元素(一侧对角线格为必换格)
for tmpRow in range(divRank):
# 交换对象线元素
for tmpCol in range(divRank):
if maskMagic[tmpRow, tmpCol] == 1:
# 小方阵的元素行列索引值要相对应
curRowA = tmpRow
curColA = tmpCol
curRowD = curRowA + divRank
curColD = curColA
tmpValue = magic[curRowA, curColA]
magic[curRowA, curColA] = magic[curRowD, curColD]
magic[curRowD, curColD] = tmpValue
break
# 其它非对角线上元素交换
rmnCnt = 0
while rmnCnt < mGrid - 1:
for tmpCol in range(divRank):
if maskMagic[tmpRow, tmpCol] != 1:
curRowA = tmpRow
curColA = tmpCol
curRowD = curRowA + divRank
curColD = curColA
tmpValue = magic[curRowA, curColA]
magic[curRowA, curColA] = magic[curRowD, curColD]
magic[curRowD, curColD] = tmpValue
# 标记已交换的状态
maskMagic[tmpRow, tmpCol] = 1
break
rmnCnt += 1
# B与C再任取m-1列相互交换
tmpColCnt = 0
while tmpColCnt < mGrid - 1:
for tmpRow in range(divRank):
curRowC = tmpRow
curColC = tmpColCnt + divRank
curRowB = curRowC + divRank
curColB = curColC
tmpValue = magic[curRowC, curColC]
magic[curRowC, curColC] = magic[curRowB, curColB]
magic[curRowB, curColB] = tmpValue
tmpColCnt += 1
print("step3:\n", magic)
return magic
def validateMagic(magic: np.ndarray, rank: int):
maxNum = rank * rank
sumTotal = (1 + maxNum) * maxNum // 2
perValue = sumTotal // rank
print(f"sumTotal={sumTotal}, perValue={perValue}")
# validate row
for tmpRow in range(rank):
tmpSum = np.sum(magic[tmpRow,0:rank])
if tmpSum != perValue:
print(f"error row:{tmpSum} != {perValue}", magic[tmpRow,0:rank])
return False
# validate col
for tmpCol in range(rank):
tmpSum = np.sum(magic[0:rank, tmpCol])
if tmpSum != perValue:
print(f"error col:{tmpSum} != {perValue}", magic[0:rank, tmpCol])
return False
# validate hypotenuse 1
tmpSum = 0
for index in range(rank):
tmpSum += magic[index,index]
if tmpSum != perValue:
print(f"error hypotenuse1: col:{tmpSum} != {perValue}")
return False
# validate hypotenuse 2
tmpSum = 0
for index in range(rank):
tmpSum += magic[rank-1-index,index]
if tmpSum != perValue:
print(f"error hypotenuse1: col:{tmpSum} != {perValue}")
return False
return True
def testOddMagic():
rank = 15
tmpMagic = oddMagic(rank)
boSuccess = validateMagic(tmpMagic, rank)
print("validate->:", boSuccess)
print(tmpMagic)
def testEvenMagicDouble():
rank = 16
tmpMagic = evenMagicDouble(rank)
boSuccess = validateMagic(tmpMagic, rank)
print("validate->:", boSuccess)
print(tmpMagic)
def testEvenMagicSingle():
rank = 14
tmpMagic = evenMagicSingle(rank)
boSuccess = validateMagic(tmpMagic, rank)
print("validate->:", boSuccess)
print(tmpMagic)
if __name__ == "__main__":
# testOddMagic()
# testEvenMagicDouble()
testEvenMagicSingle()
结论
1.奇数幻方填写容易
2.双偶数幻方网上很多参考资料都有误(包括维基百科上描述的都不全)
建议参考如下论文:
https://hanspub.org/journal/PaperInformation?paperID=64877
3.单偶数幻方填写相对复杂一点
如上代码给出三类幻方解决的一种参考方案(3-1000阶,扩展可到任意阶)
幻方解有很多种，如果要求出全部解(需要强大的计算资源，可能也计算不出来)
922
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
