算額（その1074）

算額（その1074）

九十八 江刺市 雨宝堂 現中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円2個，二等辺三角形，正方形，界斜

正方形内に二等辺三角形と甲円，乙円，界斜を容れる。乙円の直径が 1 寸のとき，界斜の長さはいかほどか。

注：算額（山村）の図では，界斜は甲円に接していないように描かれているがそんなことはない。界斜は乙円と甲円の両方に接している。

正方形の一辺の長さを 2a
界斜と正方形の一辺との交点座標を (2a, b)
界斜を延長して x 軸との交点座標を (c, 0)
甲円の半径と中心座標を r1, (a, r1)
乙円の半径と中心座標を r2, (x2, 2a - r2)
とおき，以下の連立方程式を解く。

円の中心から直線までの距離を求める際に，同じ直線でも定義に使う点の座標によって連立方程式が複雑になることがある。この算額でいえば界斜を (0, 2a) - (2a, b) で定義するのと (0, 2a) - (c, 0) で定義する場合に当たる。連立方程式が少しでも複雑になると SymPy で解けなくなってしまうので，工夫が必要になる。

図において ⊿ABX と ⊿OCX は相似で，界斜の長さは BC*2a/C になる。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, c::positive,
     r1::positive, r2::positive, x2::positive
#eq1 = dist2(0, 2a, 2a, b, a, r1, r1)
#eq2 = dist2(0, 2a, 2a, b, x2, 2a - r2, r2)
eq1 = dist2(0, 2a, c, 0, a, r1, r1)/4a
eq2 = dist2(0, 2a, c, 0, x2, 2a - r2, r2)/4a
eq3 = dist2(0, 0, a, 2a, a, r1, r1)
eq4 = dist2(0, 0, a, 2a, x2, 2a - r2, r2)
eq5 = r1*x2 - r2*(c - a)
res = solve([eq2, eq3, eq4, eq5], (a, r1, x2, c))[2];

@syms d
res[1] |> factor |> println
apart(res[2]) |> factor |>  println
apart(res[3]) |> factor |> println
res[4] |> factor |> println

   r2*(-sqrt(10*sqrt(5) + 23) + 5 + 3*sqrt(5) + sqrt(5)*sqrt(10*sqrt(5) + 23))/4
   -r2*(-3*sqrt(10*sqrt(5) + 23) - 5 - sqrt(5) + sqrt(5)*sqrt(10*sqrt(5) + 23))/4
   r2*(-sqrt(10*sqrt(5) + 23) + sqrt(5) + 3 + sqrt(5)*sqrt(10*sqrt(5) + 23))/4
   r2*(7 + 4*sqrt(5) + sqrt(5)*sqrt(10*sqrt(5) + 23))/2

以上をまとめて，a, r1, x2, c は以下のようにして求めることができる。

r2 = 1/2
t = sqrt(10√5 + 23)
a = r2*(5 - t + √5(3 + t))/4
r1 = r2*(5 + 3t + √5(1 - t))/4
x2 = r2*(3 - t + √5(1 + t))/4
c = r2*(7 + √5(4 + t))/2
(r2, a, r1, x2, c)

   r2 = 0.5;  a = 2.50415;  r1 = 1.54765;  x2 = 1.69513;  c = 7.75107

乙円の直径が 1 寸のとき，界斜の長さは sqrt(4a^2 +c^2)*(2a/c) = 5.962810866213448 寸である。

sqrt(4a^2 +c^2)*(2a/c)

   5.962810866213448

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 0.5
   t = sqrt(10√5 + 23)
   a = r2*(5 - t + √5(3 + t))/4
   r1 = r2*(5 + 3t + √5(1 - t))/4
   x2 = r2*(3 - t + √5(1 + t))/4
   c = r2*(7 + √5(4 + t))/2
   len = sqrt(4a^2 + c^2)*2a/c
   @printf("乙円の直径が %g のとき，界斜は %g である。\n", 2r2, len)
   @printf("r2 = %g;  a = %g;  r1 = %g;  x2 = %g;  c = %g\n", r2, a, r1, x2, c)
   plot([0, 2a, 2a, 0, 0], [0, 0, 2a, 2a, 0], color=:green, lw=0.5)
   plot!([0, a, 2a], [0, 2a, 0], color=:blue, lw=0.5)
   circle(a, r1, r1)
   circle(x2, 2a - r2, r2, :magenta)
   segment(0, 2a, c, 0)
   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(2a, 0, " 2a", :black, :left, :bottom, delta=delta/2, deltax=-0.5delta)
       point(2a, 2a*(c - 2a)/c, " (2a,b)", :black, :left, :vcenter)
       point(c, 0, " C", :red, :left, :bottom, delta=delta/2)
       point(a, 2a, " A", :red, :left, :bottom, delta=delta/2)
       point(0, 2a, "B ", :red, :right, :bottom, delta=delta/2)
       point(0, 0, "O ", :red, :right, :bottom, delta=delta/2)
       (x0, y0) = intersection(0, 2a, c, 0, 0, 0, a, 2a)
       point(float(x0), float(y0), "  X", :red, :left, :vcenter)
       plot!(xlims=(-5delta, c + delta))
   end
end;