A magic square of odd order can be efficiently constructed with a simple method brought from the travels of the French mathematician Simon de la LoubĂ¨re to what was then Thailand (Siam), in the 17th century. I like to think of a time when you brought things like that from overseas trips, and not only a gazillion photos and some cheap souvenirs.

The method is easier to undertand with a visual example:

A Python implementation is also very easy, especially when using Numpy, with which the `pos`

pointer variable, being a `numpy.array`

, can be both updated with one-liner arithmetic operators, and used as a two-dimensional index into the matrix (when cast as a `tuple`

):

import numpy as np def de_la_loubere(n): m = np.zeros((n, n), dtype=int) pos = np.asarray([0, n/2]) # initial m[tuple(pos)] = 1 i = 2 while i <= n ** 2: if m[tuple((pos + [-1, 1]) % n)] != 0: pos += [1, 0] # blocked: go down else: pos += [-1, 1] # go up/right pos %= n # wrap m[tuple(pos)] = i i += 1 return m n = 15 m = de_la_loubere(n) # verify that it's really magic magic = n * (n ** 2 + 1) / 2 # 65 assert (np.sum(m, axis=0) == np.repeat(magic, n)).all() assert (np.sum(m, axis=1) == np.repeat(magic, n)).all() assert m.trace() == magic assert np.fliplr(m).trace() == magic print m # [[122 139 156 173 190 207 224 1 18 35 52 69 86 103 120] # [138 155 172 189 206 223 15 17 34 51 68 85 102 119 121] # [154 171 188 205 222 14 16 33 50 67 84 101 118 135 137] # [170 187 204 221 13 30 32 49 66 83 100 117 134 136 153] # [186 203 220 12 29 31 48 65 82 99 116 133 150 152 169] # [202 219 11 28 45 47 64 81 98 115 132 149 151 168 185] # [218 10 27 44 46 63 80 97 114 131 148 165 167 184 201] # [ 9 26 43 60 62 79 96 113 130 147 164 166 183 200 217] # [ 25 42 59 61 78 95 112 129 146 163 180 182 199 216 8] # [ 41 58 75 77 94 111 128 145 162 179 181 198 215 7 24] # [ 57 74 76 93 110 127 144 161 178 195 197 214 6 23 40] # [ 73 90 92 109 126 143 160 177 194 196 213 5 22 39 56] # [ 89 91 108 125 142 159 176 193 210 212 4 21 38 55 72] # [105 107 124 141 158 175 192 209 211 3 20 37 54 71 88] # [106 123 140 157 174 191 208 225 2 19 36 53 70 87 104]]

## No comments:

## Post a Comment

Note: Only a member of this blog may post a comment.