算額（その1175）
三八 諏訪大神社 熊谷市新堀 弘化4年(1847)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，外円，斜線4本
外円の中に原点の中心を原点として，x 軸と y 軸を線対称とする 4 本の斜線を引き，区画された領域に等円 4 個を容れる。外円の直径を 1.6 寸としたとき，等円が最も大きくなるときの等円の直径を求めよ。
外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
とおき，以下の連立方程式を立てる。
include("julia-source.txt");
using SymPy
@syms R::positive, r::positive, x::positive, y::positive
eq1 = x^2 + y^2 - (R - r)^2
eq2 = r^2*(R^2 + y^2) - x^2*y^2;
eq1 は x^2 = (R - r)^2 - y^2 なので，eq2 の x^2 に代入して整理する。
eq2 = eq2(x^2 => (R - r)^2 - y^2) |> expand
eq2 |> println
R^2*r^2 - R^2*y^2 + 2*R*r*y^2 + y^4
eq2 は，R, y で r が決まるという式なので，r について解けば R, y の関数になる。
f = solve(eq2, r)[1]
f |> println
y*(R - y)/R
これは単純に R が定数のときは放物線の式であり，y がどのような値を取ると f(y) が最大になるかは計算しなくてもわかる。
たとえば R = 1.6/2 のとき y = 0.4 のとき r = 0.2 が最大値になる。
pyplot(size=(300, 200), grid=true, aspectratio=:none, label="")
plot(f(R => 1.6/2), xlims=(0, 0.8), xlabel="y", ylabel="r")
r が最大になるときの y，そして y がその値を取るときの r は以下のようにして求める。
g = diff(f, y)
y0 = solve(g, y)[1]
y0 |> println
R/2
r0 = f(y => y0)
r0 |> println
R/4
R = 1.6/2 寸 のとき，y0 = 0.4 寸, r0 = 0.2 寸ゆえ，等円の直径は 0.4 寸である。
@syms x0, y0, x, y, r, R
y0 = sqrt(R^2 - x0^2)
x = sqrt((R - r)^2 - y^2)
eq = dist2(-R, 0, x0, y0, x, y, r)
res = solve(eq, x0)[1] |> println
(R*(8*R^4 - 16*R^3*r + 8*R^3*sqrt(R^2 - 2*R*r + r^2 - y^2) + 8*R^2*r^2 - 8*R^2*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 8*R^2*y^2 + 4*R*r*y^2 - 4*R*y^2*sqrt(R^2 - 2*R*r + r^2 - y^2) - r^4 + 2*r^2*y^2) - 4*sqrt(2)*r*y*sqrt(R^3*(4*R^3 - 8*R^2*r + 4*R^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + 5*R*r^2 - 4*R*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 3*R*y^2 - r^3 + r^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + r*y^2 - y^2*sqrt(R^2 - 2*R*r + r^2 - y^2))))/(8*R^4 - 16*R^3*r + 8*R^3*sqrt(R^2 - 2*R*r + r^2 - y^2) + 12*R^2*r^2 - 8*R^2*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 4*R^2*y^2 - 4*R*r^3 + 4*R*r^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + r^4)
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(y, r) = (R/2, R/4)
x = sqrt((R - r)^2 - y^2)
t = sqrt(R^2 - 2*R*r + r^2 - y^2)
u1 = R*(8R^4 - 16R^3*r + 8R^3*t + 8R^2*r^2 - 8R^2*r*t - 8R^2*y^2 + 4R*r*y^2 - 4R*y^2*t - r^4 + 2r^2*y^2)
u2 = 4√2r*y*sqrt(R^3*(4R^3 - 8R^2*r + 4R^2*t + 5R*r^2 - 4R*r*t - 3R*y^2 - r^3 + r^2*t + r*y^2 - y^2*t))
v = 8R^4 - 16R^3*r + 8R^3*t + 12R^2*r^2 - 8R^2*r*t - 4R^2*y^2 - 4R*r^3 + 4R*r^2*t + r^4
x0 = (u1 - u2)/v
y0 = sqrt(R^2 - x0^2)
plot()
circle(0, 0, R, :red)
circle4(x, y, r, :blue)
segment(-R, 0, x0, y0, :green)
segment(-R, 0, x0, -y0, :green)
segment(R, 0, -x0, y0, :green)
segment(R, 0, -x0, -y0, :green)
#=
abline(0, y, -y/R, -R, R, :green)
abline(0, y, y/R, -R, R, :green)
abline(0, -y, -y/R, -R, R, :green)
abline(0, -y, y/R, -R, R, :green)
=#
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(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
point(0, y, "y", :blue, :center, delta=-delta/2)
point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
end
end;
draw(1.6/2, true)
