算額（その1060）

算額（その1060）

九十四 大船渡市立根町 五葉神社 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

正三角形の中に正方形と大円，小円を容れる。大円の直径が 1 寸のとき，小円の直径はいかほどか。

正三角形の一辺の長さを 2a
正方形の対角線の長さを b
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (0, b + r2)
とおき，まず eq1, eq2, eq3 の連立方程式を解き a, b, x1 を求める（4元連立方程式を一度に解けないので）。

include("julia-source.txt");

using SymPy
@syms a::positive, b,
     r1::positive, x1::positive, r2::positive
eq1 = dist2(0, 0, b/2, b/2, x1, r1, r1)
eq2 = dist2(a, 0, 0, √Sym(3)a, x1, r1, r1)
eq3 = √Sym(3)a - b/2 - √Sym(3)b/2
res = solve([eq1, eq2, eq3], (a, b, x1))[2]  # 2 of 2

   (-sqrt(6)*r1*sqrt(sqrt(3) + 2)/3 + sqrt(3)*r1/3 + r1 + sqrt(2)*r1 + sqrt(2)*r1*sqrt(sqrt(3) + 2), -sqrt(6)*r1 + 2*sqrt(2)*r1/sqrt(sqrt(3) + 2) + 2*r1 + 3*sqrt(2)*r1, r1*(1 + sqrt(2)))

得られた a, b を eq4 の a, b に代入し，方程式を解き r2 を求める。

eq4 = √Sym(3)a - b - 3r2
eq = eq4(a => res[1], b => res[2]);
ans_r2 = solve(eq, r2)[1]
ans_r2 |> println

   r1*(-2*sqrt(2) + sqrt(sqrt(3) + 2)*(-3*sqrt(2) - sqrt(2*sqrt(3) + 4) - 1 + sqrt(3) + sqrt(6*sqrt(3) + 12) + 2*sqrt(6)))/(3*sqrt(sqrt(3) + 2))

@syms d
final = apart(ans_r2, d) |> simplify |> sympy.sqrtdenest |> simplify |> factor
final |> println

   -r1*(-2*sqrt(6) - 3 + sqrt(3) + 3*sqrt(2))/3

小円の半径 r2 は，大円の半径の (2√6 + 3 - √3 - 3√2)/3 倍である。
大円の直径が 1 寸のとき，小円の直径は 0.6414293302927311 寸である。

2final(r1 => 1/2).evalf() |> println

   0.641429330292731

その他のパラメータは以下のとおりである。

   r1 = 0.5;  r2 = 0.320715;  a = 2.07313;  b = 2.62863;  x1 = 1.20711

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   t = sqrt(√3 + 2)
   (a, b, x1) = r1 .* (
       -√6t/3 + √3/3 + 1 + √2 + √2t,
       -√6 + 2√2/t + 2 + 3√2,
       1 + √2)
   r2 = r1*(2√6 + 3 - √3 - 3√2)/3
   @printf("大円の直径が %g のとき，小円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  a = %g;  b = %g;  x1 = %g\n", r1, r2, a, b, x1)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   plot!([0, b/2, 0, -b/2, 0], [0, b/2, b, b/2, 0], color=:magenta, lw=0.5)
   circle2(x1, r1, r1, :green)
   circle(0, b + r2, r2)
   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(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
       point(0, b, "b", :magenta, :center, :bottom, delta=delta/2)
       point(b/2, b/2, " (b/2,b/2)", :magenta, :left, :vcenter)
       point(0, b + r2, "小円:r2\n(0,b+r2)", :red, :center, :bottom, delta=delta/2)
       point(x1, r1, "大円:r1\n(x1,r1)", :green, :center, :bottom, delta=delta/2)
   end
end;