算額（その1140）

算額（その1140）

二五 武州妻沼 聖天宮 文政11年(1828)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，半円，斜線2本

半円の中に二等辺三角形になるように 2 本の斜線を引き，大円，中円，小円の 3 個を容れる。大円の直径が 5 寸のとき，中円の直径はいかほどか。

半円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
とおき以下の連立方程式を解く。

SymPy の性能上，一挙に連立方程式を解けないので，まず中円のパラメータを求める。

ただし，複数の条件式を使える場合，どの条件を用いれば適切な解が得られるか適切に選択する必要がある。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative, r3::positive, x3::positive
tand45by2 = √Sym(Sym(2)) - 1  # tand(Sym(45)//2)
# tand45by2 = √2 - 1  # tand(Sym(45)//2)
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = x1^2 + r1^2 - (R - r1)^2 |> expand
# eq3 = r1/(x1 + R) - tand45by2
eq3 = dist2(-R, 0, 0, R, x1, r1, r1)
# eq4 = r2 - tand45by2*(x2 + R)
eq4 = r2 - (r1/(x1 + R))*(x2 + R)
# eq4 = dist2(-R, 0, 0, R, x2, r2, r2)
##
eq5 = (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq6 = r3 - (R - x3)*tand45by2;
# eq6 = dist2(R, 0, 0, R, x3, r3, r3)

res1 = solve([eq1, eq2, eq3, eq4], (R, x1, r2, x2))[1]

   (r1*(4 + 3*sqrt(2))/4, sqrt(2)*r1/4, r1*(-4*sqrt(2) - 4*sqrt(2 - sqrt(2)) + 2*sqrt(4 - 2*sqrt(2)) + 7), r1*(-8*sqrt(4 - 2*sqrt(2)) - 8 + 9*sqrt(2))/4)

res1[1](r1 => 5/2).evalf() |> println
res1[2](r1 => 5/2).evalf() |> println
res1[3](r1 => 5/2).evalf() |> println
res1[4](r1 => 5/2).evalf() |> println

   5.15165042944955
   0.883883476483184
   1.11615673042922
   -2.45700971311331

中円の半径 r2 は，大円の半径 r1 の，(-4*sqrt(2) - 4*sqrt(2 - sqrt(2)) + 2*sqrt(4 - 2*sqrt(2)) + 7) = 0.4464626921716892 倍である。
大円の直径が 5 寸のとき，中円の直径は 5*0.4464626921716892 = 2.232313460858446 寸である。

「術」は，この倍数 0.4464626921716892 を以下のように求めている。

天 = 7 - √32
天 - sqrt(天^2 - 1)

   0.44646269217168977

この式を展開すると以下のようになる。

(7 - √32) - sqrt((7 - √32)^2 - 1)

   0.44646269217168977

求めた式と違う項は

- sqrt((7 - √Sym(32))^2 - 1) |> sympy.sqrtdenest |> simplify |> println

   -2*sqrt(20 - 14*sqrt(2))

と

2*sqrt(4 - 2*sqrt(Sym(2))) - 4*sqrt(2 - sqrt(Sym(2))) |> println

   -4*sqrt(2 - sqrt(2)) + 2*sqrt(4 - 2*sqrt(2))

であるが，後者を simplify すると

2*sqrt(4 - 2*sqrt(Sym(2))) - 4*sqrt(2 - sqrt(Sym(2))) |> simplify |> println

   -2*(2 - sqrt(2))^(3/2)

さらに (2 - √2)^3 を expand すると 20 - 14√2 になり，一致する。

(2 - √Sym(2))^3 |> expand |> println

   20 - 14*sqrt(2)

算額の「問」に答えるためにはここまででよいが，図を描くために，小円のパラメータを求める。

eq5 = eq5(R => res1[1], x1 => res1[2])
eq6 = eq6(R => res1[1], x1 => res1[2]);

res2 = solve([eq5, eq6], (r3, x3))[1]

   (sqrt(2)*r1/4 + r1/2 + (1 - sqrt(2))*(r1*(8 - 7*sqrt(2))/4 + sqrt(6)*r1*sqrt(2 - sqrt(2))), r1*(8 - 7*sqrt(2))/4 + sqrt(6)*r1*sqrt(2 - sqrt(2)))

res2[1](r1 => 5/2).evalf() |> println
res2[2](r1 => 5/2).evalf() |> println

   0.684255560080320
   3.49971137617445

r3 は若干簡約化できる。

res2[1] |> sympy.sqrtdenest |> simplify |> println

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

 全てのパラメータは以下の通りである。

   R = 5.15165; r1 = 2.5; x1 = 0.883883;  r2 = 1.11616;  x2 = -2.45701;  r3 = 0.684256;  x3 = 3.49971

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x1, r2, x2) = (r1*(4 + 3√2)/4, √2r1/4, r1*(-4√2 - 4*sqrt(2 - √2) + 2sqrt(4 - 2√2) + 7), r1*(-8sqrt(4 - 2√2) - 8 + 9√2)/4)
   (r3, x3) = (r1*(-7√2 - 4sqrt(6 - 3√2) + 2sqrt(12 - 6√2) + 12)/2, r1*(8 - 7√2)/4 + √6r1*sqrt(2 - √2))
   @printf("大円の直径が %g のとき，中円の直径は %g である。\n", 2r1, 2r2)
   @printf("R = %g; r1 = %g; x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", R, r1, x1, r2, x2, r3, x3)
   plot([-R, 0, R], [0, R, 0], color=:green, lw=0.5)
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, r1, r1, :blue)
   circle(x2, r2, r2, :magenta)
   circle(x3, r3, 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(x1, r1, "大円:r1,(x1,r1)", :blue, :center, delta=-delta/2)
       point(x2, r2, "中円:r2(x2,r2)", :magenta, :center, delta=-delta/2)
       point(x3, r3, "小円:r3\n(x3,r3)", :orange, :center, delta=-delta/2)
       point(0, R, "R", :red, :center, :bottom, delta=delta)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
   end
end;