https://sangaku0418.hatenablog.com/entry/2025/03/17/120416

宮城県白石市小原字馬頭山三瀧神社奉納算額（萬蔵稲荷神社所蔵）明治8年(1875)

徳竹亜紀子，谷垣美保：宮城県白石市小原地区の算額調査，仙台高等専門学校名取キャンパス研究紀要，第57号，2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf
キーワード：円4個，楕円，斜線
#Julia #SymPy #算額 #和算 #数学

楕円の中に斜線を 2 本引き，甲円 2 個，乙円 2 個を容れる。長径二□□寸、短径一十三寸，乙円径五寸のとき，甲円径はいかほどか。

楕円の長半径，短半径，中心座標を \(a, b, (0, 0)\)
甲円の半径と中心座標を \(r_1, (x_1, 0)\)
乙円の半径と中心座標を \(r_2, (0, b - r_2)\)
楕円と斜線の交点座標（第一象限）を \( (x_0, y_0)\)
とおき，以下の連立方程式を解く。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms a::positive, b::positive,
x0::positive, y0::positive,
r1::positive, x1::positive,
r2::positive
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = dist2(0, 0, x0, y0, x1, 0, r1)
eq3 = dist2(0, 0, x0, y0, 0, b - r2, r2)
eq4 = x1^2 - (a^2 - b^2)*(b^2 - r1^2)/b^2
res = solve([eq1, eq2, eq3, eq4], (r1, x1, x0, y0))[14] # 14 of 16

(b*sqrt(-(a - b)*(-b + 2*r2)/(a^2*b - 2*a^2*r2 + b*r2^2))*sqrt(a + b), -sqrt(b)*sqrt( (a - b)/(a^2*b - 2*a^2*r2 + b*r2^2))*sqrt(a + b)*(-b + r2), a*sqrt(b)*r2*sqrt(1/(a^2*b - 2*a^2*r2 + b*r2^2)), a*b*sqrt(-(-b + 2*r2)/(a^2*b - 2*a^2*r2 + b*r2^2)))

# r1
res[1] |> simplify

res[1](a => 26/2, b => 13/2, r2 => 5/2).evalf()*2 |> println

10.4519421124809

ルートの中身の分子は \( (a^2 - b^2)(b - 2r_2)\) と因数分解できるので，
甲円の半径は \(b\sqrt{ (a^2 - b^2)(b - 2r_2)/(a^2 b - 2a^2 r_2 + b r2^2)}\) と簡約化できる。

ans_r1 = b*sqrt( (a^2 - b^2)*(b - 2r2)/(a^2*b - 2*a^2*r2 + b*r2^2))

長径 = 26，短径 = 13，乙円径 = 5 のとき，甲円の直径は 10.4519421124809 である。

ans_r1(a => 26/2, b => 13/2, r2 => 5/2).evalf()*2 |> println

10.4519421124809

描画関数プログラムのソースを見る

function draw(a, b, r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, x0, y0) = (b*sqrt(-(a - b)*(-b + 2*r2)/(a^2*b - 2*a^2*r2 + b*r2^2))*sqrt(a + b), -sqrt(b)*sqrt( (a - b)/(a^2*b - 2*a^2*r2 + b*r2^2))*sqrt(a + b)*(-b + r2), a*sqrt(b)*r2*sqrt(1/(a^2*b - 2*a^2*r2 + b*r2^2)), a*b*sqrt(-(-b + 2*r2)/(a^2*b - 2*a^2*r2 + b*r2^2)))
@printf("長径 = %g，短径 = %g，乙円径 = %g のとき，甲円の直径は %g である。\n", 2a, 2b, 2r2, 2r1)
@printf("a = %g; b = %g; r2 = %g; r1 = %g; x1 = %g; x0 = %g; y0 = %g\n", a, b, r2, r1, x1, x0, y0)
plot()
ellipse(0, 0, a, b, color=:red)
circle2(x1, 0, r1, :blue)
circle22(0, b - r2, r2, :green)
segment(x0, y0, -x0, -y0)
segment(-x0, y0, x0, -y0)
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, 0, "a", :red, :left, :bottom, delta=delta, deltax=delta/2)
point(0, b, "b", :red, :center, :bottom, delta=delta)
point(x0, y0, "(x0,y0)", :red, :center, :bottom, delta=1.5delta)
point(x1, 0, "甲円:r1,(x1,0)", :blue, :center, delta=-delta)
point(0, b-r2, "乙円:r2\n(0,b-r2)", :green, :center, delta=-delta)
end
end;