算額（その6）新版
書き直した。
香川県高松市宮脇町 石清尾八幡宮 文化三年(1806)
http://www.wasan.jp/kagawa/iwaseo.html
千葉県富津市西大和田 吾妻神社 明治十年(1877)
http://www.wasan.jp/chiba/azuma.html
山口正義：やまぶき 第57号
https://yamabukiwasan.sakura.ne.jp/ymbk57.pdf
キーワード：円10個，外円，正三角形
#Julia, #SymPy, #算額, #和算
大円の中に，正三角形 1 個，甲円 3 個，乙円 6 個を容れる。大円の直径が 103 寸 5 分のとき，正三角形の一辺の長さを求めよ。
吾妻神社のものは，外円と甲円の直径を与えて乙円の直径を問う問題にしている。
正三角形の一辺の長さを a
大円の半径と中心座標を R, (0,0)
甲円の半径と中心座標を r2, (0, R - r2)
乙円の半径と中心座標を r1, (x1, y1)
とおき，以下の連立方程式を解く。
include("julia-source.txt")
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R, r1, x1, y1, r2
R = 1035//20
eq1 = x1^2 + (R - r2 - y1)^2 - (r1 + r2)^2
eq2 = x1^2 + y1^2 - (R - r1)^2
eq3 = dist2(0, 0, √Sym(3), 1, x1, y1, r1)
eq4 = dist2(0, R - 2r2, (R - 2r2)*√Sym(3)/2, -(R - 2r2)/2, x1, y1, r1);
まず x1, y1 を求める。
res = solve([eq1, eq2], (x1, y1))[1];
# x1
res[1] |> println
-6*sqrt(23)*sqrt(-r1*r2*(4*r1 + 4*r2 - 207))/(4*r2 - 207)
# y1
res[2] |> println
(16*r1*r2 + 828*r1 + 828*r2 - 42849)/(4*(4*r2 - 207))
eq4 の x1, y1 に代入し，r1 を求める。
eq14 = eq4(x1 => res[1], y1 => res[2]) |> simplify |> numerator;
ans_r1 = solve(eq14, r1)[1]
ans_r1 |> println
(32*r2^3 - 72*sqrt(23)*r2^2*sqrt(621 - 8*r2) - 11592*r2^2 + 3726*sqrt(23)*r2*sqrt(621 - 8*r2) + 514188*r2)/(16*r2^2 + 4968*r2 + 385641)
eq3 の x1, y1, r1 に代入し，r2 を求める。
eq13 = eq3(x1 => res[1], y1 => res[2])(r1 => ans_r1) |> simplify |> numerator;
ans_r2 = solve(eq13, r2)[6] # 6 of 7
ans_r2 |> println
207*CRootOf(1024*x^4 + 14080*x^3 - 7664*x^2 + 936*x - 27, 1)
r2 は 1024x^4 + 14080x^3 - 7664x^2 + 936x - 27 = 0 の実数解の 207 倍ということである。
実際の値は，ans_r2.evalf() で得られる。r2 = 8.78541777524225 ということである。
正三角形の一辺の長さ a は，a = (R - 2r2)*√3 なので，R = 103.5/2, r2 = 8.78541777524225 のとき，a = 59.2000493868128 である。
((R - 2r2)*√3)(R => 103.5/2, r2 => 8.78541777524225) |> println
59.2000493868128
「術」は以下のように計算を進める。
大円径 = 103.5
103.5
乾 = √3
1.7320508075688772
坤 = 57乾 + 102
200.726896031426
(坤 - sqrt(11364乾 + 19695))*大円径/4
59.20004938681292
まとめると以下のように，大円径の (57√3 + 102 - sqrt(11364√3 + 19695))/4 = 0.571981153495777 倍ということである。
大円径が 103.5 寸のとき，三角形の一辺の長さは 59.20004938681292 であり，上述の解と一致する。
(57√3 + 102 - sqrt(11364√3 + 19695))/4*大円径
59.20004938681292
function triangle(ox, oy, r; color=:blue)
θ = [90, 210, 330, 90]
x = r .* cosd.(θ)
y = r .* sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, lw=0.5)
end;
function draw(R, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 2R*0.08488326352891062
r1 = (32*r2^3 - 72*sqrt(23)*r2^2*sqrt(621 - 8*r2) - 11592*r2^2 + 3726*sqrt(23)*r2*sqrt(621 - 8*r2) + 514188*r2)/(16*r2^2 + 4968*r2 + 385641)
x1 = -6*sqrt(23)*sqrt(-r1*r2*(4*r1 + 4*r2 - 207))/(4*r2 - 207)
y1 = (16*r1*r2 + 828*r1 + 828*r2 - 42849)/(4*(4*r2 - 207))
a = √3(R - 2r2)
@printf("大円の直径が %g のとき，三角形の一辺の長さは %g である。\nまた，甲円，乙円の直径はそれぞれ %g, %g である。\n", 2R, a, 2r2, 2r1)
plot()
circle(0, 0, R, :magenta)
triangle(0, 0, R - 2r2)
rotate(0, R - r2, r2, :green)
rotate(x1, y1, r1, :red)
rotate(-r1, -(R - 2r2)/2 - r1, r1, :red)
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)
point(0, R, "R", :magenta, :center, :bottom, delta=delta)
point(0, R - r2, "甲円:r2\n(0,R-r2)", :green, :center, :bottom, delta=delta)
point(x1, y1, "乙円:r1,(x1,y1)", :red, :center, delta=-delta)
end
end
draw(103.5/2, true)
