算額（その985）
一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
楕円の中に2本の斜線を引き，区分された領域に天円 2 個，地円 1 個，人円 2 個を入れる。楕円の長径，短径がそれぞれ 69 寸，23 寸のとき，人円の直径はいかほどか。
楕円の長半径，短半径を a, b
天円の半径と中心座標を r1, (r2 + r1, 0)
地円の半径と中心座標を r2, (0, 0)
人円の半径と中心座標を r3, (r2 - r3, 0)
原点を通る，地円と人円の共通接線と楕円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。
1. 天円の半径
天円の半径は，「算法助術」の公式84により求めることができる。
楕円に内接する同じ大きさの 2 円の半径 r は，楕円の長半径 a，短半径 b，原点から円の中心までの距離 dである。
d^2 = (a^2 - b^2)*(b^2 - r^2)/b^2
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive,
x0::positive, y0::positive,
r1::positive, r2::positive, r3::positive,
x01::positive, y01::positive, d
r2 = b
d = r2 + r1
eq1 = d^2 - (a^2 - b^2)*(b^2 - r1^2)/b^2
r1 = solve(eq1, r1)[1]
r1 |> println
b - 2*b^3/a^2
r1(a => 69/2, b => 23/2).evalf() |> println
8.94444444444444
r1 は a と b の関数で，b - 2*b^3/a^2 と表すことができる。
a = 69/2, b = 23/2 を代入すると 8.94444444444444 と正しい答えになっている。
2. 人円の半径
次いで，r1 が確定したとして，人円 r3 を求める。これは，人円と天円の中心から斜線までの距離と，原点からそれぞれの中心までの距離の関係から導かれる 3元連立方程式を解くことで求めることができる。
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = r3/(r2 - r3) - y0/sqrt(x0^2 + y0^2)
eq3 = r1/(r2 + r1) - y0/sqrt(x0^2 + y0^2)
res = solve([eq1, eq2, eq3], (x0, y0, r3))[4]; # 4 of 4
println("x0: ", res[1])
println("y0: ", res[2])
println("r3: ", res[3])
x0: sqrt(a*sqrt(a + b) - (a^2 - 2*b^2)*sqrt(1/(a - b)))*sqrt(a*sqrt(a + b) + (a^2 - 2*b^2)*sqrt(1/(a - b)))/sqrt(a + b)
y0: b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a*sqrt(a + b))
r3: b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a^2*sqrt(1/(a - b)) - 2*b^2*sqrt(1/(a - b)) + 2*sqrt((a^4 - 2*a^2*b^2 + b^4)/(a^2 - b^2))*sqrt(a + b))
res[3](a => 69/2, b => 23/2).evalf() |> println
3.50000000000000
楕円の長径，短径が それぞれ 69 寸，23 寸のとき，人円の直径は 7 寸である。
3. 連立方程式を一挙に解くには
なお，SymPy の能力的には無理であるが，eq1〜eq6 を連立方程式として解けば，解は一挙に求まる。
x01, y01 は天円と楕円の接点座標である。
eq4 = x01^2/a^2 + y01^2/b^2 - 1
eq5 = (x01 - d)^2 + y01^2 - r1^2
eq6 = -b^2*x01/(a^2*y01) +(x01 - d)/y01;
--- 結果の検証のための描画プログラム
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, r2) = (69, 23) ./ 2
b = r2
r1 = b - 2*b^3/a^2
x0 = sqrt(a*sqrt(a + b) - (a^2 - 2*b^2)*sqrt(1/(a - b)))*sqrt(a*sqrt(a + b) + (a^2 - 2*b^2)*sqrt(1/(a - b)))/sqrt(a + b)
y0 = b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a*sqrt(a + b))
r3 = b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a^2*sqrt(1/(a - b)) - 2*b^2*sqrt(1/(a - b)) + 2*sqrt((a^4 - 2*a^2*b^2 + b^4)/(a^2 - b^2))*sqrt(a + b))
@printf("長径が %g 寸，短径が %g 寸のとき，人円の直径は %g 寸である。\n", 2a, 2b, 2r3)
plot()
ellipse(0, 0, a, b, color=:green)
circle2(r2 + r1, 0, r1)
circle2(r2 - r3, 0, r3, :magenta)
circle(0, 0, r2, :blue)
segment(-x0, -y0, x0, y0, :orange)
segment(-x0, y0, x0, -y0, :orange)
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(r2 + r1, 0, "天円:r1\n(r2+r1,r1)", :red, :center, delta=-2delta)
point(0, 3b/4, "地円:r2\n(0,0)", :blue, :center, delta=-2delta, mark=false)
point(r2 - r3, 0, "人円:r3\n(r2-r3,r3)", :black, :center, delta=-2delta, deltax=-4delta)
point(a, 0, " a", :green, :left, :vcenter)
point(0, b, "b = r2", :blue, :center, :bottom, delta=2delta)
end
end;
