算額（その722）
群馬県桐生市天神町 桐生天満宮 明治11年
https://gunmawasan.web.fc2.com/files/sangak-corner/tenman-g1.htm
正三角形内に 3 個の楕円が入っている。三角形の一辺の長さと楕円の短径が与えられたとき，楕円の長径を求めよ。
正三角形の重心を原点として，
正三角形の一辺の長さを 2c
楕円の長半径と短半径を 2a, 2b
下側の楕円と左上の楕円の接点座標を (x0, y0)
として，下側の楕円についての以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a, b, c, x0, y0
s3 = sqrt(Sym(3))
eq1 = x0^2/a^2 + (y0 - b + c/s3)^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - b + c/s3)) - 1/s3
eq3 = y0/x0 - 1/s3
res = solve([eq1, eq2, eq3], (a, x0, y0))
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(sqrt(c*(-6*sqrt(3)*b^2 + 9*b*c - sqrt(3)*c^2)/(3*b - sqrt(3)*c)), c*(-6*b + sqrt(3)*c)/(3*b - sqrt(3)*c), c*(-2*sqrt(3)*b + c)/(3*b - sqrt(3)*c))
(-sqrt(-c*(6*sqrt(3)*b^2 - 9*b*c + sqrt(3)*c^2)/(3*b - sqrt(3)*c)), c*(-6*b + sqrt(3)*c)/(3*b - sqrt(3)*c), c*(-2*sqrt(3)*b + c)/(3*b - sqrt(3)*c))
2 組の解が得られるが，最初のものが適解である。長半径 a の式は以下のようになる。
a = res[1][1]
a |> println
sqrt(c*(-6*sqrt(3)*b^2 + 9*b*c - sqrt(3)*c^2)/(3*b - sqrt(3)*c))
この式を簡約化するためにルートの中身について考える。
a2 = a^2
a2 |> println
c*(-6*sqrt(3)*b^2 + 9*b*c - sqrt(3)*c^2)/(3*b - sqrt(3)*c)
分数式であるが，割り切れる。結果の第1要素が商，第2要素が剰余であり，確かに剰余は0になっており，割り切れている。
temp = sympy.div(numerator(a2), denominator(a2))
(-2*sqrt(3)*b*c + c^2, 0)
楕円の長半径は商の平方根を取れば求まる。
楕円の長半径 = sqrt(temp[1])
楕円の長半径 |> println
sqrt(-2*sqrt(3)*b*c + c^2)
楕円の長半径は sqrt(-2*sqrt(3)*b*c + c^2) であるが，もう少しわかりやすく書くと，
sqrt(c*(c - √12*b))
となる。この式は「術」に述べられている「12の平方根に短径を掛けて，正三角形の一辺の長さから引き，短径を掛けて，開平する」というのと同じである。
なお，この解法の最初で「正三角形の一辺の長さを 2c，楕円の長半径と短半径を 2a, 2b」としているので，注意が必要である。
正三角形の一辺の長さが 1，楕円の短径が 0.2 のとき，c = 1/2, b = 0.2/1 を与えて長半径を得て，2倍して長径を得る。
楕円の長半径(c => 1/2, b => 0.2/2).evalf()*2 |> println
0.554237924516582
c = 1, b = 0.2 を代入すると，（2倍しなくても）長径が得られる。
楕円の長半径(c => 1, b => 0.2).evalf() |> println
0.554237924516582
パラメータをまとめると以下のようになる。
楕円の長径 = 0.554238; 長半径 = 0.277119
c = 0.5
b = 0.1
a = 0.2771189622582913
x0 = -0.2349942264126275
y0 = -0.13567397987733842
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
s3 = sqrt(3)
(c, b) = (1/2, 2/20)
(a, x0, y0) = (sqrt(-2*sqrt(3)*b*c + c^2), c*(-6*b + sqrt(3)*c)/(3*b - sqrt(3)*c), c*(-2*sqrt(3)*b + c)/(3*b - sqrt(3)*c))
@printf("楕円の長径 = %g; 長半径 = %g\n", 2a, a)
@show(c, b, a, x0, y0)
plot([c, 0, -c, c], [-c/s3, 2c/s3, -c/s3, -c/s3], color=:blue, lw=0.5)
ellipse(0, b - c/s3, a, b, color=:red)
d = c/s3 - b
ellipse(d*cosd(30), d*sind(30), a, b, color=:red, φ=120)
ellipse(-d*cosd(30), d*sind(30), a, b, color=:red, φ=240)
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, 2c/√3, " 2c/√3", :blue, :left, :vcenter)
point(c, -c/√3, " c", :blue, :left, :bottom, delta=delta/2)
point(x0, y0, "(x0,y0)", :red, :left, delta=-delta/2)
point(0, b - c/√3, " 楕円:a,b,(0,b-c/√3)", :red, :center, delta=-delta/2)
end
end;
