算額（その1449）
六十四 花泉町金沢 大門神社・大門観世音菩薩 明治26年(1893)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円1個，四分円3個，長方形
#Julia, #SymPy, #算額, #和算
問は「今有如図画勾股設等象限三個容円其円径若干問得黄積術如何」である。
補足も加え現代文風にすると，「図のように直角三角形の頂点（O, C, E）を中心とする 3 個の等四分円（注1）を描く。四分円は勾（OC）と股（OE）の辺上（A, B）で接する。3 つの四分円と弦で作られる中央部の隙間（ABDF）に『四分円と弦に接する円』(注2)を容れる。隙間の面積から円の面積を差し引いた面積（黄積と名付ける）を求める術をのべよ。」
注1：算額の図は「等四分円ではない」が，等四分円では以下の条件を満たす円は描けないのは明らかである。
注2：後に述べるが，この条件は強すぎる。
術は「位 = √5，黄積 = (sqrt(2*位 + 5)*位 - 円積率*7)*円径^2/2」
注3：円積率 = 0.79 である。当時は円の面積は直径（円形）を使って求めるのが普通であった。
円径 d の円の面積 = d^2 * 円積率 = d^2*0.79 = (r/2)^2*0.79 = r^2 *(4*0.79) = r^2 * 3.16
左下，右下，左上の四分円の半径を r1, r2, r3 とおく。
「等四分円」はありえないことは明らかであるが，試みに色々図を書いてみると，r3 = r1, r2 = r1 + r3 = 2r1 のときに図形が成立するようにみえる。
円の半径と中心座標を r, (x, y); y = r1
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1, r2, r3, r, x, y, θ1, S
r3 = r1
r2 = r1 + r3
y = r1
θ1 = atand(2r1/(r1 + r2))
eq1 = x^2 + y^2 - (r1 + r)^2
eq2 = (r1 + r2 - x)^2 + y^2 - (r2 + r)^2
# eq3 = x^2 + (r1 + r3 - y)^2 - (r3 + r)^2
# eq4 = dist2(0, r1 + r3, r1 + r2, 0, x, y, r)
res = solve([eq1, eq2], (r, x))[1] # 1 of 2
(-3*r1*(3/2 - sqrt(6)/4) + 3*r1, r1*(3/2 - sqrt(6)/4))
ans_r = res[1] |> factor
ans_r |> println
3*r1*(-2 + sqrt(6))/4
円の半径は 3*r1*(-2 + sqrt(6))/4 である。
ans_r(r1 => 1/2).evalf() |> println
0.168558653543692
黄積 S は ⊿OCE - 扇形CBD - 四分円OAB - 扇形EAF - 円 である。
∠OEC = θ1 として，以下のように計算される。
θ1 = atand(2r1/(r1 + r2))
r = 3*r1*(-2 + sqrt(6))/4
⊿OCE = (r1 + r2)*2r1/2
扇形CBD = π*r1^2*(90 - θ1)/360
四分円OAB = π*r1^2/4
扇形EAF = π*r2^2*(θ1/360)
円 = π*r^2
黄積 = ⊿OCE - 扇形CBD - 四分円OAB - 扇形EAF - 円
S = (r1 + r2)*r1 - π*r1^2*(90 - θ1)/360 - π*r1^2/4 - π*r2^2*θ1/360 - π*r^2 |> factor
S |> println
-(2*pi*r^2 - 6*r1^2 + 3*r1^2*atan(2/3) + pi*r1^2)/2
r1 = 1/2 のとき，黄積は 0.0475409496565828 である。
S(r => -3*r1*(3/2 - sqrt(6)/4) + 3*r1)(r1 => 1/2).evalf() |> println
0.0475409496565828
ここまでやってきて今更ではあるが，r1 = r3, r2 = r1 + r3 = 2r1 では，円は四分円とは接するが，弦とは接しない。差は僅かではあるが。
術は，
円積率 = 0.7854
位 = √5
黄積 = (sqrt(2*位 + 5)*位 - 円積率*7)*円径^2/2
であるが，これにより計算する場合はそもそも「円径」が必要であるが，その円径を求めるすべが書かれていない。円形が，上の近似値を取るとしても，黄径は 0.0786506985351767 になり，前述の値と一致しない。
術は明らかに不適切なものである。
@syms 円積率, 位, 円径
円積率 = 0.7854
位 = √5
黄積 = (sqrt(2*位 + 5)*位 - 円積率*7)*円径^2/2
黄積(円径 => 2*0.168558653543692) |> println
0.0786506985351767
function draw(r1, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r3 = r1
r2 = r1 + r3
y = r1
(r, x) = (-3*r1*(3/2 - sqrt(6)/4) + 3*r1, r1*(3/2 - sqrt(6)/4))
θ1 = atand(2r1/(r1 + r2))
S = (r1 + r2)*r1 - π*r1^2*(90 - θ1)/360 - π*r1^2/4 - π*r2^2*θ1/360 - π*r^2
@printf("S = %g, 0.6911*2r^2 = %g\n", S, 0.6911*2r^2)
@printf("r1 = %.15g; r2 = %.15g; r3 = %.15g; r = %.15g; x = %.15g; y = %.15g; θ1 = %.15g\n", r1, r2, r3, r, x, y, θ1)
plot([0, r1 + r2, r1 + r2, 0, 0], [0, 0, r1 + r3, r1 + r3, 0], color=:green, lw=0.5)
circle(0, 0, r1, beginangle=0, endangle=90)
circle(r1 + r2, 0, r2, :blue, beginangle=90, endangle=180)
circle(0, r1 + r3, r3, :magenta, beginangle=270, endangle=360)
circle(x, y, r, :orange)
segment(0, r1 + r3, r1 + r2, 0, :green)
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(r1 + r2, 0, "E:r1+r2", :green, :center, delta=-delta)
point(r1, 0, "A:r1", :green, :center, delta=-delta)
point(0, r1, "B:r1 ", :green, :right, :vcenter)
Dx = r1*cosd(θ1)
Dy = 2r1 - r1*sind(θ1)
point(Dx, Dy, "D", :green, :center, :bottom, delta=delta)
Fx = r1 + r2*(1 - cosd(θ1))
Fy = r2*sind(θ1)
point(Fx, Fy, "F", :green, :center, :bottom, delta=delta)
point(0, 2r1, "C:2r1 ", :green, :right, :vcenter)
point(x, y, "円:r1,(x,y)", :red, :center, delta=-delta)
point(0, 0, "O ", :red, :right, delta=-delta)
plot!(xlims=(-13delta, r1 + r2 + 5delta), ylims=(-6delta, 2r1 + 2delta))
end
end;
draw(1/2, true)
