算額（その681）
山口正義: やまぶき4，第58号，2018/12/06.
千葉県君津市鹿野山 鹿野山神野寺 万延二年(1861)
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf
山口正義：やまぶき，和算と歴史随想 からリンク
https://yamabukiwasan.sakura.ne.jp/page3.html
短径が 3 寸の楕円の中に，直径 2寸4分の等円が 2 個入っている。楕円の周長はいかほどか。ただし，円積率として 7分8厘5毛4糸を使うものとする。
楕円の長半径，短半径，中心座標を a, b, (0, 0)
等円の半径と，中心座標を r, (r, 0)
楕円と等円の接点の座標を (x1, y1)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, r::positive, x1::positive, y1::positive
eq1 = (x1 - r)^2 + y1^2 - r^2
eq2 = -b^2*x1/(a^2*y1) - (r - x1)/y1
eq3 = x1^2/a^2 + y1^2/b^2 - 1
res = solve([eq1, eq2, eq3], (x1, y1, a))
2-element Vector{Tuple{Sym, Sym, Sym}}:
(b^2/r, b*sqrt(-b^2 + 2*r^2)/r, -b^2*sqrt(1/(b^2 - r^2)))
(b^2/r, b*sqrt(-b^2 + 2*r^2)/r, b^2*sqrt(1/(b^2 - r^2)))
2 組の解が得られるが，2 番目のものが適解である（最初のものは x 軸に関して線対称な解である）。
x1, y2, a はそれぞれ短半径 b の b/r 倍，sqrt(2r^2 - b^2) 倍，b/sqrt(b^2 - r^2) 倍である。
短径（差渡し径）が 3 寸（短半径が 1.5 寸），等円の直径が 2.4 寸のとき，長半径は 2.5 寸である（差渡し径 5 寸）。
b = 1.5; r = 1.2; x1 = 1.875; y1 = 0.992157; a = 2.5
楕円の周長は第二種完全楕円積分による正確な値として 12.763499431699064 が得られる。
問において，「円周率は 0.7854 を使う」とあるのは，円積率のことで円周率になおすと 3.1416 である。当時，円積率としては 0.79（円周率は 3.16）が使われることが多かったであろうにしては，正確すぎる。
しかし，この円積率を使ったとしても，「答」にある 12.7634994317113 有奇を導いた近似式などが不明（術を読み解く力がない）。有効数字 11 桁が一致している。驚異的。
長径と短径の比がきれいな値なので，事前に計算された数表などがあったのかもしれない。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(b, r) = (3//2, 24/20)
(x1, y1, a) = b .* (b/r, sqrt(2r^2 - b^2)/r, b/sqrt(b^2 - r^2))
@printf("b = %g; r = %g; x1 = %g; y1 = %g; a = %g\n", b, r, x1, y1, a)
plot()
ellipse(0, 0, a, b, color=:blue)
circle(r, 0, r)
circle(-r, 0, r)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
point(x1, y1, " (x1,y1)", :black, :left, :bottom)
point(r, 0, "r", :red, :center, :bottom, delta=delta)
point(0, b, " b", :blue, :left, :bottom, delta=delta)
point(a, 0, " a", :blue, :left, :bottom, delta=delta)
end
end;
