算額（その673）

算額（その673）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問題. 100　19✕19魔方陣

関孝和が考案した方式があるらしいが，多分追跡できないので，ヒンズー方式，バシェー方式の二通りで書いてみた。

function generate_odd_magic_square(n)
   # ヒンズーの連続方式
   m = zeros(Int, n, n)
   (i, j) = (div(n + 1, 2), 1)
   m[j, i] = 1
   for num in 2:n^2
       (next_i, next_j) = (i + 1, j - 1)
       next_i > n && (next_i = 1)
       next_j < 1 && (next_j = n)
       m[next_j, next_i] != 0 && ((next_i, next_j) = (i, j + 1))
       m[next_j, next_i] = num
       (i, j) = (next_i, next_j)
   end
   return m
end;

function generate_odd_magic_square2(n)
   m = zeros(Int, n, n)
   (i, j) = (n - 1, div(n - 1, 2))
   for num in 1:n^2
       m[i + 1, j + 1] = num
       (next_i, next_j) = ((i + 1) % n + 1, (j + 1) % n + 1)
       if m[next_i, next_j] == 0
           (i, j) = (next_i - 1, next_j - 1)
       else
           i -= 1
       end
   end
   return m
end;

function generate_odd_magic_square3(n)
   # バシェー方式
   a = zeros(Int, n, n)
   m = n ÷ 2
   for i = 1:n
       (sx, sy) = (i - 1 - m, n - i - m)
       for j = 1:n
           (x, y) = (sx + j, sy + j)
           x <= 0 && (x += n)
           x >  n && (x -= n)
           y <= 0 && (y += n)
           y >  n && (y -= n)
           a[x, y] = (i - 1) * n + j
       end
   end
   return a
end;

function print_magic_square(m)
   n = size(m, 1)
   for i in 1:n
       for j in 1:n
           print(m[i, j], "\t")
       end
       println()
   end
end;

function check_magic_square(magic_square)
   n = size(magic_square, 1)
   sums = n*(n^2 + 1)/2
   allequal(sum(magic_square, dims=1) .== sums) &&
   allequal(sum(magic_square, dims=2) .== sums) &&
   sum(magic_square[i, i] for i in 1:n) == sums &&
   sum(magic_square[i, n + 1 -i] for i in 1:n) == sums
end;

n = 19
magic_square = generate_odd_magic_square(n);
magic_square2 = generate_odd_magic_square2(n);
magic_square3 = generate_odd_magic_square3(n);

print_magic_square(magic_square);
println("--")
print_magic_square(magic_square2);
println("--")
print_magic_square(magic_square3);
println("--")

check_magic_square(magic_square) |> println
check_magic_square(magic_square2) |> println
check_magic_square(magic_square3) |> println

   192 213 234 255 276 297 318 339 360 1 22 43 64 85 106 127 148 169 190 
   212 233 254 275 296 317 338 359 19 21 42 63 84 105 126 147 168 189 191 
   232 253 274 295 316 337 358 18 20 41 62 83 104 125 146 167 188 209 211 
   252 273 294 315 336 357 17 38 40 61 82 103 124 145 166 187 208 210 231 
   272 293 314 335 356 16 37 39 60 81 102 123 144 165 186 207 228 230 251 
   292 313 334 355 15 36 57 59 80 101 122 143 164 185 206 227 229 250 271 
   312 333 354 14 35 56 58 79 100 121 142 163 184 205 226 247 249 270 291 
   332 353 13 34 55 76 78 99 120 141 162 183 204 225 246 248 269 290 311 
   352 12 33 54 75 77 98 119 140 161 182 203 224 245 266 268 289 310 331 
   11 32 53 74 95 97 118 139 160 181 202 223 244 265 267 288 309 330 351 
   31 52 73 94 96 117 138 159 180 201 222 243 264 285 287 308 329 350 10 
   51 72 93 114 116 137 158 179 200 221 242 263 284 286 307 328 349 9 30 
   71 92 113 115 136 157 178 199 220 241 262 283 304 306 327 348 8 29 50 
   91 112 133 135 156 177 198 219 240 261 282 303 305 326 347 7 28 49 70 
   111 132 134 155 176 197 218 239 260 281 302 323 325 346 6 27 48 69 90 
   131 152 154 175 196 217 238 259 280 301 322 324 345 5 26 47 68 89 110 
   151 153 174 195 216 237 258 279 300 321 342 344 4 25 46 67 88 109 130 
   171 173 194 215 236 257 278 299 320 341 343 3 24 45 66 87 108 129 150 
   172 193 214 235 256 277 298 319 340 361 2 23 44 65 86 107 128 149 170 
   --
   172 193 214 235 256 277 298 319 340 361 2 23 44 65 86 107 128 149 170 
   171 173 194 215 236 257 278 299 320 341 343 3 24 45 66 87 108 129 150 
   151 153 174 195 216 237 258 279 300 321 342 344 4 25 46 67 88 109 130 
   131 152 154 175 196 217 238 259 280 301 322 324 345 5 26 47 68 89 110 
   111 132 134 155 176 197 218 239 260 281 302 323 325 346 6 27 48 69 90 
   91 112 133 135 156 177 198 219 240 261 282 303 305 326 347 7 28 49 70 
   71 92 113 115 136 157 178 199 220 241 262 283 304 306 327 348 8 29 50 
   51 72 93 114 116 137 158 179 200 221 242 263 284 286 307 328 349 9 30 
   31 52 73 94 96 117 138 159 180 201 222 243 264 285 287 308 329 350 10 
   11 32 53 74 95 97 118 139 160 181 202 223 244 265 267 288 309 330 351 
   352 12 33 54 75 77 98 119 140 161 182 203 224 245 266 268 289 310 331 
   332 353 13 34 55 76 78 99 120 141 162 183 204 225 246 248 269 290 311 
   312 333 354 14 35 56 58 79 100 121 142 163 184 205 226 247 249 270 291 
   292 313 334 355 15 36 57 59 80 101 122 143 164 185 206 227 229 250 271 
   272 293 314 335 356 16 37 39 60 81 102 123 144 165 186 207 228 230 251 
   252 273 294 315 336 357 17 38 40 61 82 103 124 145 166 187 208 210 231 
   232 253 274 295 316 337 358 18 20 41 62 83 104 125 146 167 188 209 211 
   212 233 254 275 296 317 338 359 19 21 42 63 84 105 126 147 168 189 191 
   192 213 234 255 276 297 318 339 360 1 22 43 64 85 106 127 148 169 190 
   --
   172 353 154 335 136 317 118 299 100 281 82 263 64 245 46 227 28 209 10 
   11 173 354 155 336 137 318 119 300 101 282 83 264 65 246 47 228 29 191 
   192 12 174 355 156 337 138 319 120 301 102 283 84 265 66 247 48 210 30 
   31 193 13 175 356 157 338 139 320 121 302 103 284 85 266 67 229 49 211 
   212 32 194 14 176 357 158 339 140 321 122 303 104 285 86 248 68 230 50 
   51 213 33 195 15 177 358 159 340 141 322 123 304 105 267 87 249 69 231 
   232 52 214 34 196 16 178 359 160 341 142 323 124 286 106 268 88 250 70 
   71 233 53 215 35 197 17 179 360 161 342 143 305 125 287 107 269 89 251 
   252 72 234 54 216 36 198 18 180 361 162 324 144 306 126 288 108 270 90 
   91 253 73 235 55 217 37 199 19 181 343 163 325 145 307 127 289 109 271 
   272 92 254 74 236 56 218 38 200 1 182 344 164 326 146 308 128 290 110 
   111 273 93 255 75 237 57 219 20 201 2 183 345 165 327 147 309 129 291 
   292 112 274 94 256 76 238 39 220 21 202 3 184 346 166 328 148 310 130 
   131 293 113 275 95 257 58 239 40 221 22 203 4 185 347 167 329 149 311 
   312 132 294 114 276 77 258 59 240 41 222 23 204 5 186 348 168 330 150 
   151 313 133 295 96 277 78 259 60 241 42 223 24 205 6 187 349 169 331 
   332 152 314 115 296 97 278 79 260 61 242 43 224 25 206 7 188 350 170 
   171 333 134 315 116 297 98 279 80 261 62 243 44 225 26 207 8 189 351 
   352 153 334 135 316 117 298 99 280 81 262 63 244 45 226 27 208 9 190 
   --
   true
   true
   true

 

算額（その672）

算額（その672）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問. 61　大円の内外に中円，小円を 19 個配置する。大円の直径が 1 尺のとき，中円，小円の直径はいかほどか。

プログラム的には外円に内接する小円を求めるのが簡単である。
外円（半径 R，中心座標 (0, 0)）に外接する中円の半径，内接する小円の半径をそれぞれ r1, r2 とする。以下の方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms θ::positive, R::positive, r1::positive, r2::positive, r3::positive, n::integer

θ = 360/n
eq1 = (R + r1)*sind(θ/2) - r1
eq2 = (R - r2)*sind(θ/2) - r2
solve([eq1, eq2], (r1, r2)) |> println

   Dict{Any, Any}(r2 => R*sin(pi/n)/(sin(pi/n) + 1), r1 => -R*sin(pi/n)/(sin(pi/n) - 1))

r1 = R*sin(pi/n)/(1 - sin(pi/n))
r2 = R*sin(pi/n)/(1 + sin(pi/n))
である。
外円の直径が = 1 尺 = 10 寸のとき，中円の直径は 1.9702361077126025 寸，小円の直径は 1.413320924331764 寸である。

R = 10/2
n = 19
r1 = R*sin(pi/n)/(1 - sin(pi/n))
r2 = R*sin(pi/n)/(1 + sin(pi/n))
(r1, r2) .* 2 |> println

   (1.9702361077126025, 1.413320924331764)

同じプログラムで，外円の直径が R, 円の数が n の図形を順次縮小して図を描くことができる(maxn=2 を指定する）。

function draw(R =10, n=19, maxn=2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   θ = 360/n
   x = cosd.(range(0, 360, n + 1))
   y = sind.(range(0, 360, n + 1))
   r1 = -R*sin(pi/n)/(sin(pi/n) - 1)
   R += 2r1
   plot()
   for j = 1:maxn
       if maxn == 2 || j != 1
           println("外円の直径 = $R,  外円に内接する円の直径 = $r1")
           for i in 1:n
               circle((R - r1)*x[i], (R - r1)*y[i], r1, :blue)
           end
       end
       maxn == 2 && j == 2 && break
       R -= 2r1
       r1 = R*sin(pi/n)/(sin(pi/n) + 1)
       circle(0, 0, R)
   end
   if more && maxn == 2
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       #vline!([0], color=:black, lw=0.5)
       #hline!([0], color=:black, lw=0.5)
       point(R, 0, " R", :red, :left, :vcenter)
       point(R + 1.97, 0, " R+r1", :red, :left, :vcenter)
       point(R - 1.41, 0, "R-r2 ", :red, :right, :vcenter)
   end
end;

算額（その671）

算額（その671）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問. 25　正三角形に内接する大小 2 個の円がある。小円の直径が 1 寸のとき，正三角形の一辺の長さはいかほどか。

正三角形の一辺の長さを 2a とする。
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (0, 2r1 + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, a::positive
s3 = sqrt(Sym(3))
eq1 = r2/(s3*a - (2r1 + r2)) - 1//2
eq2 = r1/(s3*a - r1) - 1//2
res = solve([eq1, eq2], (r1, a))

   Dict{Any, Any} with 2 entries:
     r1 => 3*r2
     a  => 3*sqrt(3)*r2

正三角形の一辺の長さは小円の半径の 6√3 倍である。
小円の直径が 1 寸のとき，三角形の一辺の長さは 1/2 * 6√3 = 5.196152422706632 寸である。
ちなみに，そのときの大円の直径は小円の直径の 3 倍なので，3 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (r1, a) = r2 .* (3, 3√3)
   @printf("正三角形の一辺の長さ = %g;  大円の直径 = %g\n", 2a, 2r1)
   @printf("r1 = %g;  a = %g\n", r1, a)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, 2r1 + r2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, r1, " 大円:r1,(0,r1)", :red, :left, :vcenter)
       point(0, 2r1 + r2, " 小円:r2,(0,2r1+r2)", :blue, :left, :vcenter)
       point(0, √3a, " √3a", :black, :left, :vcenter)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
   end
end;