算額（その951）
一六 武州 一之宮（武州赤坂氷川明神） 文化12年(1815)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
大円の中に甲円，乙円，丁円がそれぞれ 2 個，丙円が 4 個入っている。
黒積（図の灰色で示した面積の 8 倍）が 950625 歩であるとき，丁円の直径はいかほどか。
この類の問題は，黒積の算出に用いた円周率が何なのかとか問題中の数値の誤記などで，算額の答えと合わないことがママあるので，ここでは連立方程式の次数を下げる目的と合わせ，黒積の数値を与えずに丁円の直径（半径）を未知数のまま他のパラメータを求め，しかる後に丁円の直径の確定値がいくつのときに与えられた黒積の数値になるかという順に進める。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (R - r2, 0)
丙円の半径と中心座標を r3, (x2, y3); x3 = R - r2, y3 = R - r1
丁円の半径と中心座標を r4, (r4, 0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
@syms R::positive, r1::positive, r2::positive,
r3::positive, x3::positive, y3::positive,
r4::positive, 円周率
x3 = R - r2
y3 = R - r1
eq1 = x3^2 + y3^2 - (R - r3)^2
eq2 = (R - r2)^2 + (R -r1)^2 - (r1 + r2)^2
eq3 = r4^2 + (R - r1)^2 - (r1 + r4)^2
eq4 = 2r4 + 2r2 - R
res = solve([eq1, eq2, eq3, eq4], (R, r1, r2, r3))
2-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
(2*r4*(2 + sqrt(5)), 2*r4*(2*sqrt(5) + 5)/5, r4*(1 + sqrt(5)), r4*(7 + 19*sqrt(5)/5))
(2*r4*(2 + sqrt(5)), 2*r4*(2*sqrt(5) + 5)/5, r4*(1 + sqrt(5)), r4*(sqrt(5)/5 + 1))
丁円の半径を与えれば，外円，甲円，乙円，丙円の半径が求まる。
丁円の半径を 1 とすれば，
外円の半径 R = 2*r4*(2 + sqrt(5))
甲円の半径 r1 = 2*r4*(2*sqrt(5) + 5)/5
乙円の半径 r2 = r4*(1 + sqrt(5))
丙円の半径 r3 = r4*(sqrt(5)/5 + 1)
黒積は，図の灰色で示した面積の 8 倍である。灰色の面積は，外円の面積の(中心角θ/360)倍の扇形の面積から，直角三角形と乙円の面積の 1/4 と丙円の面積の (θ + 90)/360 を除いたものである。
tan(θ) = y3/x3
r4 = 1
R = 2*r4*(2 + sqrt(5))
r1 = 2*r4*(2*sqrt(5) + 5)/5
r2 = r4*(1 + sqrt(5))
r3 = r4*(sqrt(5)/5 + 1)
x3 = R - r2
y3 = R - r1
黒積1 = 8(pi*R^2*atand(y3/x3)/360 - x3*y3/2 - pi*r2^2/4 - pi*r3^2*(atand(y3/x3) + 90)/360)
26.35149033973981
問での図形は，「黒積は 950625」とあるので，r4 = 1 のときの黒積と比較すると図形は √(950625/26.35149033973981) = 189.9336982367789 倍である。
すなわち，丁円の半径は 189.9336982367789 寸，直径は 379.8673964735578 である。
その他のパラメータは以下のとおりである。
R = 1609.14; r1 = 719.631; r2 = 614.638; r3 = 274.875; x3 = 994.506; y3 = 889.513; r4 = 189.934
このとき，確かに黒積2は 950625 になる。
算額の答えは「丁円径 358 寸有奇」で，食い違っている。
r4 = 189.9336982367789
R = 2*r4*(2 + sqrt(5))
r1 = 2*r4*(2*sqrt(5) + 5)/5
r2 = r4*(1 + sqrt(5))
r3 = r4*(sqrt(5)/5 + 1)
x3 = R - r2
y3 = R - r1
黒積2 = 8(pi*R^2*atand(y3/x3)/360 - x3*y3/2 - pi*r2^2/4 - pi*r3^2*(atand(y3/x3) + 90)/360)
950625.0000000008
「答」も「黒積」も正しいとすると，用いた円周率が 3.22534682526815 ということになってしまう。
@syms R, r1, r2, r3, x3, y3, r4
r4 = 358//2 # 「答」とされた値
R = 2*r4*(2 + sqrt(5))
r1 = 2*r4*(2*sqrt(5) + 5)/5
r2 = r4*(1 + sqrt(5))
r3 = r4*(sqrt(5)/5 + 1)
x3 = R - r2
y3 = R - r1
eq99 = 8(円周率*R^2*atand(y3/x3)/360 - x3*y3/2 - 円周率*r2^2/4 - 円周率*r3^2*(atand(y3/x3) + 90)/360) - 950625
solve(eq99, 円周率)[1] |> println
3.22534682526815
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r4 = 189.9336982367789
(R, r1, r2, r3) = (2*r4*(2 + sqrt(5)), 2*r4*(2*sqrt(5) + 5)/5, r4*(1 + sqrt(5)), r4*(sqrt(5)/5 + 1))
x3 = R - r2
y3 = R - r1
println(4*R^2*atan((R - r1)/(R - r2)) - 2*pi*r2^2 - pi*r3^2*(180*atan((R - r1)/(R - r2))/pi + 90)/45 - 4*(R - r1)*(R - r2))
println(2r4)
@printf("R = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g; y3 = %g; r4 = %g\n", R, r1, r2, r3, x3, y3, r4)
plot()
circle(0, 0, R, :black)
circle22(0, R - r1, r1)
circle2(R - r2, 0, r2, :blue)
circle4(x3, y3, r3, :green)
circle2(r4, 0, r4, :magenta)
segment(x3, 0, x3, y3, :black)
θ = atand(y3, x3)
#segment(x3, y3, x3 + r3*cosd(θ), y3 + r3*sind(θ), :black)
segment(0, 0, x3 + r3*cosd(θ), y3 + r3*sind(θ), :black)
circle(0, 0, 0.9r4, :black, beginangle= 0, endangle=θ)
circle(0, 0, 0.8r4, :black, beginangle= 0, endangle=θ)
x = []
y = []
for theta = 0:0.05:θ
append!(x, R*cosd(theta))
append!(y, R*sind(theta))
end
for theta = θ:-0.05:-90
append!(x, x3 + r3*cosd(theta))
append!(y, y3 + r3*sind(theta))
end
for theta = 90:-0.05:0
append!(x, x3 + r2*cosd(theta))
append!(y, r2*sind(theta))
end
plot!(x, y, seriestype=:shape, fillcolor=:gray80)
segment(0, 0, R, 0, :black)
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.9r4, 0.7r4, "θ", :black, mark=false)
point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
point(R - r2, 0, "乙円:r2,(R-r2,0)", :blue, :center, delta=-delta/2)
point(x3, y3, "丙円:r3,(x3,y3)", :green, :right, :bottom, delta=1.5delta, deltax=3delta)
point(r4, 0, "丁円:r4,(r4,0)", :magenta, :center, delta=-delta/2, deltax=-4delta)
point(x3, r2, "(x3,r2)", :black, :right, delta=-2delta)
end
end;
