算額（その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;

算額（その1073）

算額（その1073）

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

外円の中に，甲円 2 個，乙円 4 個，丙円 2 個，丁円 4 個を容れる。丁円の直径が 1 寸のとき，丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0); R - 2r1
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, 0)
丁円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。

一度に解けない。しかし，丙円は最終的には丁円との関連を問われるが，甲円，乙円，丁円は丙円とは独立に決めることができるので，まず eq1, eq2, eq3, eq5, eq7 の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive,
     r4::positive, x4::positive, y4::positive
R = 2r1
eq1 = x2^2 + r2^2 - (R -r2)^2 |> expand
eq2 = x4^2 + y4^2 - (R - r4)^2 |> expand
eq3 = x2^2 + r2^2 - (r1 + r2)^2 |> expand
eq4 = x3^2 + (R - r1)^2 - (r1 + r3)^2 |> expand
eq5 = x4^2 + (y4 - R + r1)^2 - (r1 + r4)^2 |> expand
eq6 = (x2 - x3)^2 + r2^2 - (r2 + r3)^2 |> expand
eq7 = (x2 - x4)^2 + (y4 - r2)^2 - (r2 + r4)^2 |> expand;
res = solve([eq1, eq2, eq3, eq5, eq7], (r1, r2, x2, x4, y4))[1]

   (sqrt(2)*r4 + 5*r4/2, sqrt(2)*r4/2 + 5*r4/4, 2*r4 + 5*sqrt(2)*r4/2, 2*r4*sqrt(2*sqrt(2) + 3), 2*r4*(1 + sqrt(2)))

x4 は二重根号を外すことができる。

res[1] |> simplify |> println
res[2] |> simplify |> println
res[3] |> simplify |> println
res[4] |> sympy.sqrtdenest |> println
res[5] |> simplify |> println

   r4*(2*sqrt(2) + 5)/2
   r4*(2*sqrt(2) + 5)/4
   r4*(4 + 5*sqrt(2))/2
   2*r4*(1 + sqrt(2))
   2*r4*(1 + sqrt(2))

r1, r2, x2, x4, y4 が決まったので，残りの r3, x3 を求める。

r1 = √Sym(2)r4 + 5r4/2
r2 = √Sym(2)r4/2 + 5r4/4
x2 = 2r4 + 5√Sym(2)r4/2
x4 = 2r4*(1 + √Sym(2))
y4 = 2r4 + 2√Sym(2)r4
R = 2r1
eq4 = x3^2 + (R - r1)^2 - (r1 + r3)^2 |> expand
eq6 = (x2 - x3)^2 + r2^2 - (r2 + r3)^2 |> expand
res2 = solve([eq4, eq6], (r3, x3))[1]

   (-5*r4 - 2*sqrt(2)*r4 + 2*sqrt(2)*(8*r4/7 + 10*sqrt(2)*r4/7), 8*r4/7 + 10*sqrt(2)*r4/7)

r3 は更に簡約化できる。

res2[1] |> simplify |> println
res2[2] |> simplify |> println

   r4*(2*sqrt(2) + 5)/7
   2*r4*(4 + 5*sqrt(2))/7

丙円の半径 r3 は，丁円の半径 r4 の (2√2 + 5)/7 倍である。
丁円の直径が 1 寸のとき，丙円の直径は (2√2 + 5)/7 = 1.1183467321065985 寸である。

(2√2 + 5)/7

   1.1183467321065985

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

   r4 = 0.5;  r1 = 1.95711;  r2 = 0.978553;  x2 = 2.76777;  x4 = 2.41421;  y4 = 2.41421;  r3 = 0.559173;  x3 = 1.58158;  R = 3.91421

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r4 = 1/2
   r1 = √2r4 + 5r4/2
   r2 = √2r4/2 + 5r4/4
   x2 = r4*(2 + 5√2/2)
   x4 = 2r4*(1 + √2)
   y4 = 2r4*(1 + √2)
   r3 = r4*(2√2 + 5)/7
   x3 = 2r4*(4 + 5√2)/7
   R = 2r1
   @printf("丁円の直径が %g のとき，丙円の直径は %g である。\n", 2r4, 2r3)
   @printf("r4 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  x4 = %g;  y4 = %g;  r3 = %g;  x3 = %g;  R = %g\n",
       r4, r1, r2, x2, x4, y4, r3, x3, R)
   plot()
   circle(0, 0, R, :orange)
   circle22(0, R - r1, r1)
   circle4(x2, r2, r2, :blue)
   circle2(x3, 0, r3, :magenta)
   circle4(x4, y4, r4, :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(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(x2, r2, "乙円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
       point(x3, 0, "丙円:r3,(x3,0)", :magenta, :left, delta=-delta/2, deltax=-4delta)
       point(x4, y4, "丁円:r4,(x4,y4)", :green, :left, delta=-delta/2, deltax=-4delta)
       point(R, 0, " R", :black, :left, :bottom, delta=delta/2, deltax=-0.5delta)
   end
end;

算額（その1072）

算額（その1072）

九十六 大船渡市立根 気仙安養寺稲荷堂 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正方形の中に楕円，元円，利円，貞円，享円を容れる。貞円の直径が 1 寸のとき，元円の直径はいかほどか。

注1：利円，享円の条件については何も述べられていない。実際，制限はあるもののどんな値も取りうる。
注2：算額の図（山村の図）では，楕円が楕円らしくないが，それはよくあることかもしれないが，結果に影響を与えるかもしれない。

正方形の一辺の長さを 2a
楕円の長半径，短半径，中心座標を a, b, (0, b); b = a - r1
元円の半径と中心座標を r1, (0, 2a - r1), (r1, a), (r1, a - 2r1)
貞円の半径と中心座標を r2, (a - r3, r3)
楕円と元円の接点座標を (x1, y1)
楕円と貞円の接点座標を (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive
b = a - r1
eq1 = x2^2/a^2 + (y2 - a + r1)^2/b^2 - 1
eq2 = (x2 - a + r2)^2 + (y2 - r2)^2 - r2^2
eq3 = -b^2*x2/(a^2*(y2 - a + r1)) - (a - r2 - x2)/(y2 - r2)
eq4 = x1^2/a^2 + (y1 - a + r1)^2/b^2 - 1
eq5 = (x1 - r1)^2 + (y1 - a)^2 - r1^2
eq6 = -b^2*x1/(a^2*(y1 - a + r1)) - (r1 - x1)/(y1 - a);

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, r1, x1, y1, x2, y2) = u
   return [
       -1 + (-a + r1 + y2)^2/(a - r1)^2 + x2^2/a^2,  # eq1
       -r2^2 + (-r2 + y2)^2 + (-a + r2 + x2)^2,  # eq2
       -(a - r2 - x2)/(-r2 + y2) - x2*(a - r1)^2/(a^2*(-a + r1 + y2)),  # eq3
       -1 + (-a + r1 + y1)^2/(a - r1)^2 + x1^2/a^2,  # eq4
       -r1^2 + (-a + y1)^2 + (-r1 + x1)^2,  # eq5
       -(r1 - x1)/(-a + y1) - x1*(a - r1)^2/(a^2*(-a + r1 + y1)),  # eq6
   ]
end;

r2 = 1/2
iniv = BigFloat[2.1, 0.63, 0.83, 2.6, 1.5, 0.52]
res = nls(H, ini=iniv)

   ([3.5742972321311055, 1.1362300699229417, 1.470601088978576, 4.6602138004791485, 2.755899335723288, 0.8855162507120506], true)

元円の直径は 2*1.1362300699229417 = 2.2724601398458835 である。「答」では 2.243 寸有奇となっており，かなり異なった結論になった。

ちなみに，楕円の長半径，短半径は 3.5743, b = 2.43807 とこれまた，あまりきれいな数ではない。

2*1.1362300699229417

   2.2724601398458835

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

   r2 = 0.5;  a = 3.5743;  b = 2.43807;  r1 = 1.13623;  r2 = 0.5;  x1 = 1.4706;  y1 = 4.66021;  x2 = 2.7559;  y2 = 0.885516

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, r1, x1, y1, x2, y2) = res[1]
   b = a - r1
   @printf("貞円の直径が %g のとき，元円の直径は %g である。\n", 2r2, 2r1)
   @printf("r2 = %g;  a = %g;  b = %g;  r1 = %g;  r2 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n",
       r2, a, b, r1, r2, x1, y1, x2, y2)
   plot([a, a, -a, -a, a], [0, 2a, 2a, 0, 0], color=:green, lw=0.5)
   ellipse(0, a - r1, a, b, color=:blue)
   circle(0, 2a - r1, r1)
   circle2(r1, a, r1)
   circle2(r1, a - 2r1, r1)
   circle2(a - r2, r2, r2, :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(x2,  y2, "(x2,y2) ", :green, :right, :bottom, delta=2delta, deltax=4delta)
       point(x1,  y1, "(x1,y1)", :green, :left, :bottom, delta=delta/2)
       point(r1, a, "元円:r1,(r1, a)", :red, :center, delta=-delta/2)
       point(r1, a - 2r1, "元円:r1,(r1, a-2r1)", :red, :center, :bottom, delta=delta/2)
       point(0, 2a - r1, "元円:r1,(0,2a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, r2, "貞円:r3,(a-r3,r3)", :green, :center, delta=-delta/2)
       point(0, 2a - 2r1, "2a-2r1", :blue, :center, delta=-delta/2)
       point(0, a - r1, "a-r1", :blue, :center, delta=-delta/2)
   end
end;