算額（その1644）
三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
全円の中に合同な直角三角形を 4 個，大円 4 個，中円 2 個，小円 4 個を容れる。小円の直径が与えられたとき，中円の直径を得る術を述べよ。
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, 0), (0, √3r1)
中円の半径と中心座標を r2, (R - r2, 0)
とおき，以下の方程式を順に解いていく。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, r1::positive, r2::positive,
r3::positive, x3::positive, y3::positive,
l::positive;
# r1
eq1 = dist2(0, R, R/√Sym(3), 0, r1, 0, r1)
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
R*(2 - sqrt(3))
# r2
eq2 = r2/(R - r2 - R/√Sym(3)) - 1//2
eq2 = 2r2 - (R - r2 - R/√Sym(3))
ans_r2 = solve(eq2, r2)[1]
ans_r2 |> println
R*(3 - sqrt(3))/9
# (x0, y0)
eq3 = (R/√Sym(3) + l*cosd(Sym(30)))^2 + (l*sind(Sym(30)))^2 - R^2
ans_l = solve(eq3, l)[1]
ans_l |> println
x0 = R/√Sym(3) + ans_l*cosd(Sym(30)) |> simplify
y0 = ans_l*sind(Sym(30))
(x0, y0) |> println
R*(-3 + sqrt(33))/6
(sqrt(3)*R*(1 + sqrt(33))/12, R*(-3 + sqrt(33))/12)
# r3, x3, y3
eq4 = dist2(0, R, x0, y0, x3, y3, r3)
eq5 = y3/x3*(y0 - R)/x0 + 1
eq5 = y3/x3 + x0/(y0 - R)
eq6 = x3^2 + y3^2 - (R - r3)^2
(ans_r3, ans_x3, ans_y3) = solve([eq4, eq5, eq6], (r3, x3, y3))[1]
# r3
ans_r3 |> println
R*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144)/12288
中円の直径は，小円の直径の (R*(3 - sqrt(3))/9) / (R*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144)/12288) = 1.30332288183929 倍である。
p = ans_r2/ans_r3
p |> println
4096*(3 - sqrt(3))/(3*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144))
p.evalf() |> println
1.30332288183929
山村の術の解説は以下のように，「中円径 = 1.1882963312409864 * 小円径」となっているが，比は小さすぎる。
また，無理数は √3，√5 しか出てこない。
@syms 天, A, 小円径, 中円径
小円径 = 1
天 = √3
A = sqrt((√5天 - 3)*2)
# 中円径 = (天 - 1)*2/((天 - A)*3)*小円径
# 中円径/小円径
(天 - 1)*2/((天 - A)*3) |> println
1.1882963312409864
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = R*(2 - sqrt(3))
plot([0, R/√3, 0, -R/√3, 0], [R, 0, -R, 0, R], color=:blue, lw=0.5)
circle(0, 0, R, :green)
circle2(r1, 0, r1)
circle22(0, √3r1, r1)
l = R*(-3 + sqrt(33))/6
(x0, y0) = (R/√3 + l*cosd(30), l*sind(30))
plot!([0, x0, R/√3], [R, y0, 0], color=:blue, lw=0.5)
plot!([0, x0, R/√3], -[R, y0, 0], color=:blue, lw=0.5)
plot!(-[0, x0, R/√3], [R, y0, 0], color=:blue, lw=0.5)
plot!(-[0, x0, R/√3], -[R, y0, 0], color=:blue, lw=0.5)
r2 = R*(3 - sqrt(3))/9
circle2(R - r2, 0, r2, :magenta)
(r3, x3, y3) = (R*(-43*sqrt(14814 - 1794*sqrt(33)) - 15*sqrt(54318 - 6578*sqrt(33)) + 6144)/12288, R*(128*sqrt(3) + 9*sqrt(4938 - 598*sqrt(33)) + 384*sqrt(11) + 7*sqrt(162954 - 19734*sqrt(33)))/6144, sqrt(33)*R/48 + 5*sqrt(33)*R*sqrt(14814 - 1794*sqrt(33))/12288 + 3*R/16 + 43*R*sqrt(14814 - 1794*sqrt(33))/12288)
circle4(x3, y3, r3, :orange)
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, 0, "大円:r1\n(r1,0)", :red, :center, delta=-delta/2)
point(0, √3r1, "大円:r1\n(0,√3r1)", :red, :center, delta=-delta/2)
point(R - r2, 0, "中円:r2,(R-r2,0)", :black, :center, delta=-delta/2)
point(x3, y3, "小円:r3,(x3,y3)", :black, :center, delta=-delta/2)
point(x0, y0, "(x0,y0)", :blue, :right, :vcenter, deltax=-2delta)
point(R/√3, 0, "R/√3", :blue, :left, :bottom, delta=2.5delta, deltax=-1.5delta)
point(0, R, "R", :green, :center, :bottom, delta=delta/2)
end
end;
draw(1, true)
