算額（その981）
一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
正方形内に四分円 3 個，甲円 1 個，乙円 1 個，丙円 2 個を入れる。乙円の直径が 3.1 寸（注）のとき，丙円の直径はいかほどか。
注：後述の通り，2.1 寸は 3.1 寸の誤り。
正方形の一辺の長さを 2a
甲円の半径と中心座標を r1, (a, a)
乙円の半径と中心座標を r2, (x2, 2a - x2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。
正方形，甲円，乙円を決定する方程式 eq3, eq4, eq5 は独立なので，まず a, r1, x2 を求める。
include("julia-source.txt");
using SymPy
@syms a::positive, r1::positive,
r2::positive, x2::positive,
r3::positive, x3::positive, y3::positive
# eq1 = (a - x3)^2 + (y3 - a)^2 - (r1 + r3)^2 |> expand
# eq2 = (2a - x3)^2 + (2a - y3)^2 - (2a - r3)^2 |> expand
eq3 = 2a^2 - (2a - r1)^2 |> expand
eq4 = (2a - x2)^2 + x2^2 - (2a - r2)^2 |> expand
eq5 = 2(x2 - a)^2 - (r1 + r2)^2 |> expand;
# eq6 = dist2(2a, 0, 0, 2a, x3, y3, r3)
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (a, r1, x2, r3, x3, y3))
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
-2*a^2 + 4*a*r1 - r1^2, # eq3
4*a*r2 - 4*a*x2 - r2^2 + 2*x2^2, # eq4
2*a^2 - 4*a*x2 - r1^2 - 2*r1*r2 - r2^2 + 2*x2^2, # eq5
res = solve([eq3, eq4, eq5], (a, r1, x2))[2] # 2 of 2
(r2*(2 + 3*sqrt(2))/2, r2*(-1 + 2*sqrt(2)), r2*(3*sqrt(2)/2 + 3))
(res[1])(r2 => 3.1/2).evalf() |> println
(res[2])(r2 => 3.1/2).evalf() |> println
(res[3])(r2 => 3.1/2).evalf() |> println
4.83804653251745
2.83406204335659
7.93804653251745
a, r1, x2 が既知として，eq1, eq2, eq6 から 最終的に r3 を求める（r2 を含む式）。
include("julia-source.txt");
using SymPy
@syms a::positive, r1::positive,
r2::positive, x2::positive,
r3::positive, x3::positive, y3::positive
a = r2*(2 + 3*sqrt(Sym(2)))/2
r1 = r2*(-1 + 2*sqrt(Sym(2)))
x2 = r2*(3*sqrt(Sym(2))/2 + 3)
eq1 = (a - x3)^2 + (y3 - a)^2 - (r1 + r3)^2 |> expand
eq2 = (2a - x3)^2 + (2a - y3)^2 - (2a - r3)^2 |> expand
eq6 = dist2(2a, 0, 0, 2a, x3, y3, r3);
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq6, ", # eq6")
2*r2^2 + 10*sqrt(2)*r2^2 - 4*sqrt(2)*r2*r3 + 2*r2*r3 - 3*sqrt(2)*r2*x3 - 2*r2*x3 - 3*sqrt(2)*r2*y3 - 2*r2*y3 - r3^2 + x3^2 + y3^2, # eq1
12*sqrt(2)*r2^2 + 22*r2^2 + 4*r2*r3 + 6*sqrt(2)*r2*r3 - 6*sqrt(2)*r2*x3 - 4*r2*x3 - 6*sqrt(2)*r2*y3 - 4*r2*y3 - r3^2 + x3^2 + y3^2, # eq2
6*sqrt(2)*r2^2 + 11*r2^2 - 3*sqrt(2)*r2*x3 - 2*r2*x3 - 3*sqrt(2)*r2*y3 - 2*r2*y3 - r3^2 + x3^2/2 + x3*y3 + y3^2/2, # eq6
eq1, eq2 から x3, y3 を求める。
res = solve([eq1, eq2], (x3, y3))[1] # 1 of 2
(-r2 + 2*sqrt(2)*r2 - 2*sqrt(r3)*sqrt(6733*sqrt(2)*r2 + 9608*r2 - 2194*r3 - 1513*sqrt(2)*r3)/(58 + 45*sqrt(2)) - sqrt(2)*r3/2 + 2*r3, -r2 + 2*sqrt(2)*r2 + 2*sqrt(r3)*sqrt(6733*sqrt(2)*r2 + 9608*r2 - 2194*r3 - 1513*sqrt(2)*r3)/(58 + 45*sqrt(2)) - sqrt(2)*r3/2 + 2*r3)
res[1](r2=>3.1/2, r3 => 1.0019922445804257).evalf() |> println
res[2](r2=>3.1/2, r3 => 1.0019922445804257).evalf() |> println
1.51119910595402
6.74786293740257
eq6 に x3, y3 を代入する。
eq6 = eq6(x3 => res[1], y3 => res[2]) |> simplify;
ans_r3 = solve(eq6, r3)[2]
ans_r3 |> println
-r2*sqrt(9 - 4*sqrt(2))/4 - sqrt(2)*r2*sqrt(9 - 4*sqrt(2))/8 + sqrt(2)*r2/8 + 5*r2/4
簡約化する。
ans_r3 = ans_r3 |> sympy.sqrtdenest |> simplify
ans_r3 |> println
r2*(4 - sqrt(2))/4
正しい答えが得られることを確認する。
ans_r3(r2 => 3.1/2).evalf() |> println
1.00199224458043
丙円の半径は，乙円の半径の (4 - sqrt(2))/4 倍である。
乙円の直径が 3.1 寸のとき，丙円の直径は 3.1 * (4 - sqrt(2))/4 = 2.0039844891608514 である。
逆にたどって，全てのパラメータの解を求める。
r2 = 3.1/2
r3 = r2*(4 - √2)/4
x3 = -r2 + 2√2*r2 - 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
y3 = -r2 + 2√2*r2 + 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
a = r2*(2 + 3√2)/2
r1 = r2*(2√2 - 1)
x2 = r2*(3√2/2 + 3)
r2 = 1.55; a = 4.83805; r1 = 2.83406; x2 = 7.93805; r3 = 1.00199; x3 = 1.5112; y3 = 6.74786
丙円の半径は 0.6787689398770625（直径は 1.357537879754125）
「答」は「丙円径は2寸00有奇」となっている。
「術」は「(1-sqrt(0.5)/2)*乙円径」と言っているのに。計算すると 1.3575378797541253 になる。
乙円の径が 3.1 寸のとき, 丙円の径が 2.0039844891608514 になる。
結局，「答」，「術」は正しく，「問」の「只云乙円径二寸一分」が「只云乙円径三寸一分」の誤記ということだろう。算額が間違えているのか，出版時の校正漏れなのかはわからない。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 3.1/2
r3 = r2*(4 - √2)/4
x3 = -r2 + 2√2*r2 - 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
y3 = -r2 + 2√2*r2 + 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
a = r2*(2 + 3√2)/2
r1 = r2*(2√2 - 1)
x2 = r2*(3√2/2 + 3)
@printf("乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r2, 2r3)
@printf("r2 = %g; a = %g; r1 = %g; x2 = %g; r3 = %g; x3 = %g; y3 = %g\n", r2, a, r1, x2, r3, x3, y3)
plot([0, 2a, 2a, 0, 0], [0, 0, 2a, 2a, 0], color=:blue, lw=0.5)
circle(0, 0, 2a, beginangle=0, endangle=90)
circle(2a, 2a, 2a, beginangle=180, endangle=270)
circle(0, 2a, 2a, beginangle=270, endangle=360)
circle(a, a, r1, :green)
circle(x2, 2a - x2, r2, :brown)
circle(x3, y3, r3, :magenta)
circle(2a - y3, 2a - x3, r3, :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)
point(a, a, "甲円:r1,(a,a)", :green, :center, delta=-delta/2)
point(x2, 2a - x2, "乙円:r2,(x2,2a-x2)", :brown, :center, delta=-delta/2)
point(x3, y3, "丙円:r3\n(x3,y3)", :magenta, :center, delta=-delta/2)
point(2a, 2a, "(2a,2a)", :blue, :right, :bottom, delta=delta/2)
end
end;
