算額（その1608）
福島県田村市船引町堀越明神前 明石神社(堀越明宮) 明治11年(1878)
http://www.wasan.jp/fukusima/horikosiakasi.html
～落書き帳「○△□」～ 696. 第10回「街角の問題」―堀越村明宮神社算額（その３）
http://streetwasan.web.fc2.com/math19.9.27.html
キーワード：円7個，直線上
#Julia, #SymPy, #算額, #和算, #数学
直線上に，東円 1 個，西円 2 個，南円 1 個，北円 3 個を載せる。東円の直径が 5 寸のとき北円の直径はいかほどか。
注：horikosiakasi.html の写真では，「北円径一寸四分□東円径幾何」とあり，「答曰東円径五寸」となっている。
東円の半径と中心座標を r1, (0, 2r3 - r1)
西円の半径と中心座標を r2, (r2, r2)
南円の半径と中心座標を r3, (0, r3)
北円の半径と中心座標を r4, (x4, y4), (0, r4)
とおき，以下の連立方程式を解く。
まず，右にある西円と中央の 北円が外接する条件 eq4 で，r2 を求めると，r2 = 4r4 であることがわかる。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive, r3::positive, r4::positive, x4::positive, y4::positive
eq4 = r2^2 + (r2 - r4)^2 - (r2 + r4)^2;
ans_r2 = solve(eq4, r2)[1]
ans_r2 |> println
4*r4
r2 = 2r4 である。また，y4 = 2r4 でもある。
これを条件として加え，残りの 3 方程式を解く。
horikosiakasi.html の「北円径一寸四分□東円径幾何」に基づき（r4 が与えられたとき），r1, r3, x4 を求める。
r2 = 4r4
y4 = 2r4
eq1 = r2^2 + (2r3 - r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = x4^2 + (r3 - y4)^2 - (r3 + r4)^2 |> expand
eq3 = (x4 - r2)^2 + (y4 - r2)^2 - (r4 - r2)^2 |> expand
res = solve([eq1, eq2, eq3], (r1, r3, x4))[2] # 2 of 2
(r4*(7*sqrt(5) + 27)/12, 4*r4*(sqrt(5) + 3)/3, r4*(sqrt(5) + 4))
# r1 東円
res[1] |> println
r4*(7*sqrt(5) + 27)/12
東円の半径 r1 は，北円の半径 r4 の (7√5 + 27)/12 = 3.5543729868748777 倍である。
北円の直径が 1 寸 4 分のとき，東円の直径は 1.4*(7√5 + 27)/12 = 4.976122181624828 寸である。
「術」は「北径を 0.28 で割れ」とあり，1.4/0.28 = 5 である。
実にきれいな結果になるが，0.28 は，r4 に掛ける倍数とすれば 1/0.28 = 3.571428571428571 なので，(7√5 + 27)/12 = 3.5543729868748777 の近似値ではなく，答えが 5 になるように選んだものであろう。
ということで，ちょっと汚いので，「落書き帳」ではこれを逆にしたのであろう。
「連立方程式を解く」というメソッドを取っているので，方程式はそのままで，何を求めるかを変えればよいだけである。
r1 が与えられるとして，r3, r4, x4 を求めればよい。
res2 = solve([eq1, eq2, eq3], (r3, r4, x4))[2] # 2 of 2
(24*sqrt(5)*r1/121 + 184*r1/121, -21*sqrt(5)*r1/121 + 81*r1/121, 6*r1*(73/242 - sqrt(5)/242))
# r3
res[1] |> simplify |> println
r4*(7*sqrt(5) + 27)/12
# r4
res[2] |> simplify |> println
4*r4*(sqrt(5) + 3)/3
# x4
res[3] |> simplify |> println
r4*(sqrt(5) + 4)
北円の半径 r4 は，東円の半径 r1 の 3(27 - 7√5)/121 である。
東円の直径が 5 寸のとき，北円の直径は 1.4067178707646453 寸である。
これならば，「小数点以下 1 桁まで取って，北円の直径は 1 寸 4 分である」としても，さほどの違和感はない。（しかし，算額の原文を変更するのはあまり好ましいものとは思えない）
5*3(27 - 7√5)/121
1.4067178707646453
function draw(r1, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r3 = r1*8(3√5 + 23)/121
r4 = r1*3(27 - 7√5)/121
x4 = r1*3(73 - √5)/121
r2 = 4r4
y4 = 2r4
@printf("r1 = %g; r2 = %g; r3 = %g; r4 = %g; x4 = %g; y4 = %g\n", r1, r2, r3, r4, x4, y4)
plot()
circle(0, r3, r3)
circle2(r2, r2, r2, :blue)
circle(0, 2r3 - r1, r1, :green)
circle2(x4, y4, r4, :magenta)
circle(0, r4, r4, :magenta)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
hline!([r4, 2r4, 3r4, 4r4], color=:orange, lw=0.5, style=:dot)
point(0, 2r3 - r1, "東円:r1,(0,2r3-r1)", :green, :center, delta=-delta/2)
point(r2, r2, "西円:r2,(r2,r2)", :blue, :center, delta=-delta/2)
point(0, r3, "南円:r3,(0,r3)", :red, :center, delta=-delta/2)
point(x4, y4, "北円:r4,(x4,y4)", :magenta, :center, delta=-delta/2)
point(0, r4, "北円:r4,(0,r3)", :magenta, :center, delta=-delta/2)
end
end;
draw(5/2, true)
