算額（その1648）
三十二 一関市舞川 観福寺内地蔵堂 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
今有如図 03086
https://w.atwiki.jp/sangaku/pages/158.html
キーワード：円3個，半円，斜線2本
#Julia, #SymPy, #算額, #和算
水平線上に合同な直角三角形（緑）を縦・横に配し，甲乙の 2 本の斜線（紫）を引いてそれぞれの直角三角形の面積を等分割する。分割された区画に大円 1 個，小円 2 個を容れる。小円の直径が与えられたとき，大円の直径を求める術を述べよ。
山村の図は細部が不適切なので（算額へのリスペクトがまったくない），今回も「今有如図」が示す図に基づいて解を求める。その結果，術と同じ解が得られた。山村の術の解説には誤りがあることもわかった。
直角三角形の直角を挟む二辺のうち，短い方を「鈎」，長い方を「股」とする。
大円の半径と中心座標を r1, (x1, y1)
小円の半径と中心座標を r2, (股/2, r2), (股 + 鈎 - r2, 股/2)
とおき，以下の連立法手式を解く。
鈎（または股）が既知として，順次パラメータを確定する手順を示し，最後に一挙に連立方程式を解く手順を示す。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms 鈎::positive, 股::positive, l::positive,
r1::positive, x1::positive, y1::positive,
r2::positive;
まず最初に，図のように甲斜，乙斜（2 本の斜線）がそれぞれの直角三角形の面積を等分割するということは，甲斜は 2 つの合同な直角三角形の弦（斜辺）の中点 A, B を結ぶ直線が左側の直角三角形の直角の頂点 O を通らなければならない。すなわち，鈎，股（直角を挟む二辺）の間に特定の関係が成り立っている必要がある。
ここでは仮に，鈎が事前に決められている場合について考える。
eq1 = 鈎/股 - (股/2)/(股 + 鈎/2);
# 股
ans_股 = solve(eq1, 股)[1]
ans_股 |> println
鈎*(1 + sqrt(2))
股は鈎の (1 + √2) 倍でなければならない。
# 鈎
ans_鈎 = solve(eq1, 鈎)[1]
ans_鈎 |> println
股*(-1 + sqrt(2))
逆に，股が事前に決められている場合には 鈎は股の (√2 - 1)倍でなければならない。
2変数の間の関係式（方程式）は同じで，どちらを未知数として解くかが違うだけである。
次に，小円が内接している 2 つの二等辺三角形は合同であり，小円の半径は弦の長さによって決まる（言うまでもなく，弦は鈎，股によって決まる）。
# r2
弦 = sqrt(鈎^2 + 股^2)
eq2 = 鈎*股/4 - (股 + 弦/2 + 弦/2)*r2/2
ans_r2 = solve(eq2, r2)[1]
ans_r2 |> println
股*鈎/(2*(股 + sqrt(股^2 + 鈎^2)))
また，大円が内接している三角形は，直角二等辺三角形であり，内接円の半径 r1 は，等辺 = 弦/2，斜辺 = √2等辺 により決定される。
2r1 = 等辺 + 等辺 - √2等辺
# r1
等辺 = 弦/2
eq3 = 2r1 - (等辺 + 等辺 - √Sym(2)等辺)
ans_r1 = solve(eq3, r1)[1]
ans_r1 |> println
(2 - sqrt(2))*sqrt(股^2 + 鈎^2)/4
変数は 鈎，股，r1, r2 の 4 変数で，方程式は eq1, eq2, eq3 の 3 本である。
変数のうちの 1 つが既知で，残りの 3 つが未知として連立方程式を立てて一挙に解を求めることができる。
上の例では，鈎が既知で，股，r1, r2 を求める場合は以下のようになる。
res1 = solve([eq1, eq2, eq3], (股, r1, r2))[1]
(鈎*(1 + sqrt(2)), 鈎*(-1 + sqrt(2))*sqrt(sqrt(2) + 2)/2, 鈎*(1 + sqrt(2))/(2*(1 + sqrt(2) + sqrt(2)*sqrt(sqrt(2) + 2))))
同じ連立方程式で，r2 が 既知で，r1，鈎，股を求める場合は以下のようになる。
算額では，鈎，股を求めることは要求していないが，既知の r2 から未知の r1 を求める術になっている。
res2 = solve([eq1, eq2, eq3], (r1, 鈎, 股))[1]
(r2*(-2 - sqrt(sqrt(2) + 2) + sqrt(2*sqrt(2) + 4) + 2*sqrt(2)), 2*r2*(-sqrt(2*sqrt(2) + 4) + 1 + 2*sqrt(sqrt(2) + 2)), 2*r2*(1 + sqrt(2) + sqrt(2*sqrt(2) + 4)))
# r1
res2[1] |> factor
大円の半径 r1 は，小円の半径 r2 の (-1 + sqrt(2))*(sqrt(sqrt(2) + 2) + 2) = 1.59379398947637 倍である。
小円の直径が 1 のとき，大円の直径は 1.59379398947637 である。
# r1
res2[1] |> factor |> println
r2*(-1 + sqrt(2))*(sqrt(sqrt(2) + 2) + 2)
res2[1](r2 => 1).evalf() |> println
1.59379398947637
「術」は以下のように「大円径は小円径の 1.5937939894763695 倍」としている。上に述べた結果と同じである。
なお，山村は術の解説で，「方斜率 = √3」としている。方斜率は √2 なので，とんでもない間違いである。山村の本は，いろんな間違いがたくさんある，稀代の悪書である。
# 術
小円径 = 1
方斜率 = √2
天 = 2 - √2
大円径 = (√天 + 天*方斜率)*小円径
1.5937939894763695
図を描くためだけに大円の中心座標 (x1, y1) の解析解を求めるのは，労多くして益少しなので数値解を求める。
eq1 = 鈎/股 - (股/2)/(股 + 鈎/2);
弦 = sqrt(鈎^2 + 股^2)
eq2 = 鈎*股/4 - (股 + 弦/2 + 弦/2)*r2/2
等辺 = 弦/2
eq3 = 2r1 - (等辺 + 等辺 - √Sym(2)等辺)
eq4 = dist(0, 鈎, 股, 0, x1, y1) - r1^2
eq5 = dist(股, 0, 股 + 鈎, 股, x1, y1) - r1^2;
function H(u)
(r1, 鈎, 股, x1, y1) = u
return [
-股/(2*(股 + 鈎/2)) + 鈎/股, # eq1
-r2*(股 + sqrt(股^2 + 鈎^2))/2 + 股*鈎/4, # eq2
2*r1 - sqrt(股^2 + 鈎^2) + sqrt(2)*sqrt(股^2 + 鈎^2)/2, # eq3
-r1^2 + (x1 - 股*(x1*股 - 鈎*(y1 - 鈎))/(股^2 + 鈎^2))^2 + (y1 - 鈎 + 鈎*(x1*股 - 鈎*(y1 - 鈎))/(股^2 + 鈎^2))^2, # eq4
-r1^2 + (y1 - 股*(y1*股 + 鈎*(x1 - 股))/(股^2 + 鈎^2))^2 + (x1 - 股 - 鈎*(y1*股 + 鈎*(x1 - 股))/(股^2 + 鈎^2))^2, # eq5
]
end;
r2 = 1 # 0.11769
iniv = BigFloat[1.59, 4.16, 10.05, 9.19, 2.08]
res = nls(H, ini=iniv)
([1.5937939894763693, 4.164784400584788, 10.054678984251696, 9.192123892710637, 2.082392200292394], true)
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, 鈎, 股, x1, y1) = r2 .* [1.5937939894763693, 4.164784400584788, 10.054678984251696, 9.192123892710637, 2.082392200292394]
plot([0, 股 + 鈎, 股 + 鈎, 股, 0, 0],
[0, 0, 股, 0, 鈎, 0], color=:green, lw = 0.5)
plot!([0, 股 + 鈎/2, 股 + 鈎],
[0, 股/2, 0], color=:magenta, lw = 0.5)
circle(x1, y1, r1)
circle(股/2, r2, r2, :blue)
circle(股 + 鈎 - r2, 股/2, r2, :blue)
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(股/2, 鈎/2, "A", :green, :center, :bottom, delta=delta)
point(股 + 鈎/2, 股/2, " B", :green, :left, :vcenter)
point(0, 0, "O", :green, :center, delta=-2delta)
point(0, 鈎, "鈎 ", :green, :right, :vcenter)
point(股, 0, "股", :green, :center, delta=-2delta)
point(股 + 鈎, 0, "股+鈎", :green, :center, delta=-2delta)
point(股 + 鈎, 股, "(股+鈎,股) ", :green, :right, :vcenter)
point(x1, y1, "大円:r1\n(x1,y1)", :red, :center, delta=-delta)
point(股/2, r2, "小円:r2,(股/2,r2)", :blue, :center, delta=-delta, deltax=-10delta)
point(股 + 鈎 - r2, 股/2, "小円:r2,(股+鈎-r2,股2)", :blue, :right, :bottom, delta=delta, deltax=-8delta)
plot!(xlims=(-8delta, 股 + 鈎 + 5delta), ylims=(-8delta, 股 + 5delta))
end
end;
draw(1/2, true)
