算額（その1523）
四十 岩手県一関市牧沢 牧沢八幡神社 明治5年(1872)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
今有如図 03060
https://w.atwiki.jp/sangaku/pages/308.html
キーワード：円5個，斜線2本
#Julia, #SymPy, #算額, #和算
菱形の中に合同な 3 個の楕円を容れる。楕円の長径が 2 寸，短径が 1 寸 のとき，菱形の長い方の対角線はいかほどか。
山村の図は正確ではない上に，問題文で「長径三寸」と書いているのは「長径二寸」の誤りである（長径三寸では，答があわない）。「今有如図」の正しい図（3 つの楕円が合同に見えないのは錯覚。下図と同じ）に基づいて解を求める。
菱形の対角線の長い方を「長」，短い方を「平」
楕円の長半径を p，短半径を q
とおき，以下の連立方程式を解く。
菱形の辺を延長し，x = -q および x = q の直線による 2 個の二等辺三角形を考える。
大きい二等辺三角形には緑の楕円が内接し，小さい二等辺三角形には赤の楕円が内接する。
これらに，三角形に内接する楕円と三角形の辺と高さに関する「算法助術の公式97」を適用する。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms A::positive, B::positive, C::positive, H::positive,
P::positive, Q::positive
eq97 = -(B^2 - C^2)^2*H*Q^2 + (B^2 - C^2)^2*Q^3 + A^4*H*(2H - Q)^2 - A^4*Q*(2H - Q)^2 - A^2*H*P^2*(2H - Q)^2;
@syms a::positive, b::positive, 長::positive, 平::positive,
p::positive, q::positive
b = 平/2
a = 長/2
factor1 = (a - q/2)/a
bc1 = sqrt((平/2*factor1)^2 + (長/2*factor1)^2)
eq1 = eq97(A => 平*factor1, H => 長/2*factor1, B => bc1, C => bc1, P => q, Q => p);
factor2 = (a + q/2)/a
bc2 = sqrt((平/2*factor2)^2 + (長/2*factor2)^2)
eq2 = eq97(A => 平*factor2, H => 長/2*factor2, B => bc2, C => bc2, P => p, Q => q);
res = solve([eq1, eq2], (平, 長))[2] # 平と長を逆にすると有限の時間内に求まらない
(sqrt((p - q)/(p^2 - p*q + q^2))*(-2*p^2 + p*q - q^2)/(2*(p - q)*sqrt(1/(p - q))), (2*p^2 - p*q + q^2)/(p - q))
p，q を長半径，短半径とすると，平，長として菱長，菱平の 1/2 が求まる。
菱長/2 は (2*p^2 - p*q + q^2)/(p - q)
res[2] |> println
res[2](p => 2/2, q => 1/2) |> println
(2*p^2 - p*q + q^2)/(p - q)
3.50000000000000
術は，「長径の二乗の二倍を，長径と短径の差で割り，短径□□を引く」で，長径が 2，短径が 1 のとき，菱長は 7 になり，答えと一致する。
山村は「短径□□」を「短径二段（二倍）」と読み取ったようだ。
「今有如図」は「額文は『和算 岩手の現存算額のすべて』によった」として，単に「短径」としている。
また，山村は「長径三寸」と読み取ったようだが，「今有如図」では「長径二寸」としている。
結果的には，山村は二重に間違っている。
2長径^2/(長径 - 短径) - 短径 = (2長径^2 - 長径*短径 - 短径^2)/(長径 - 短径) なので，上で得られた (2*p^2 - p*q + q^2)/(p - q) と同じである。
長径 = 2
短径 = 1
菱長1 = 2長径^2/(長径 - 短径) - 短径 # = (2長径^2 - 長径*短径 - 短径^2)/(長径 - 短径)
7.0
村山は問で，「長径三寸」と書いているので，上の術では答えと合わないので「術を修正」している。
長径 = 3
短径 = 1
菱長2 = 2長径^2/(長径 - 短径) - 2短径 |> simplify
7.0
どちらの術が正しいのか，長径 5，短径 2 のときに計算してみる。
2res[2](p => 5/2, q => 2/2) |> println # res[2] は 菱長/2 なので，2倍する
14.6666666666667
長径 = 5
短径 = 2
菱長1 = 2長径^2/(長径 - 短径) - 短径 |> println
菱長2 = 2長径^2/(長径 - 短径) - 2短径 |> println
14.666666666666668
12.666666666666668
算額の術が正しい。
村山はいつもの悪い癖で，問の解釈を誤った図，あるいは問の読み取りの誤りによる解が「答・術」に一致しない場合に，自分の解を正当化するために「答・術」の方を修正するという暴挙を繰り返している。
function draw(p, q, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
# p, q は，長半径，短半径，得られるのは 平/2， 長/2
(平2, 長2) = (sqrt((p - q)/(p^2 - p*q + q^2))*(2*p^2 - p*q + q^2)/(2*(p - q)*sqrt(1/(p - q))), (2*p^2 - p*q + q^2)/(p - q))
@printf("楕円の長径と短径は %g と %g; 菱形の対角線の長さは %g と %g\n", 2p, 2q, 2平2, 2長2)
plot([長2, 0, -長2, 0, 長2], [0, 平2, 0, -平2, 0], color=:blue, lw=0.5)
ellipse(0, 0, q, p, color=:green)
ellipse((p + q), 0, p, q, color=:red)
ellipse(-(p + q), 0, p, q, color=:red)
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, 平2, "平/2", :blue, :left, :bottom, delta=delta)
point(長2, 0, "長/2", :blue, :left, :bottom, delta=delta)
point(q + p, 0, "q+p", :red, :center, delta=-delta)
point(q, 0, "q ", :red, :right, delta=-delta)
point(-q, 0, " -q", :red, :left, delta=-delta)
xlims!(-長2 - 5delta, 長2 + 15delta)
end
end;
draw(2/2, 1/2, true)
