算額（その1103）

算額（その1103）

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

全円の中に水平な弦と斜線を引き，甲円 3 個，乙円 2 個，丙円 2 個を容れる。丙円の直径が 8 寸のとき，乙円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (0, R - 3r1), (0, r1 - R)
乙円の半径と中心座標を r2,(x2, y2); y2 = -r1
丙円の半径と中心座標を r3, (x3, R - 2r1 + r3)
斜線と全円の交点座標を (x0, -sqrt(R^2 - x0^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms R::positive, r1::positive, r2::positive,
     x2::positive, y2::positive, r3::positive,
     x3::positive, x0::positive
y2 = -r1  # (R - 2r1) - (2R - 2r1)/2
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x3^2 + (R - 2r1 + r3)^2 - (R - r3)^2
eq4 = dist2(sqrt(R^2 - (R - 2r1)^2), R - 2r1, -x0, -sqrt(R^2 - x0^2), 0, R - 3r1, r1) |> numerator
eq5 = dist2(sqrt(R^2 - (R - 2r1)^2), R - 2r1, -x0, -sqrt(R^2 - x0^2), 0, r1 - R, r1) |> numerator
eq6 = dist2(sqrt(R^2 - (R - 2r1)^2), R - 2r1, -x0, -sqrt(R^2 - x0^2), x2, y2, r2) |> numerator;

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)
   (R, r1, r2, x2, x3, x0) = u
   s = sqrt(R^2 - x0^2)
   t = sqrt(R - r1)
   u = sqrt(r1)
   v = r1^(3/2)
   w = r1^(5/2)
   return [
       x3^2 + (r1 - r3)^2 - (r1 + r3)^2,  # eq1
       r1^2 + x2^2 - (R - r2)^2,  # eq2
       x3^2 - (R - r3)^2 + (R - 2*r1 + r3)^2,  # eq3
       8R^5*r1 - 42R^4*r1^2 + 8R^4*r1*s + 4R^3*v*x0*t + 78R^3*r1^3 - 42R^3*r1^2*s - 6R^3*r1*x0^2 - 26R^2*w*x0*t + 4R^2*v*x0*t*s - 46R^2*r1^4 + 78R^2*r1^3*s + 47R^2*r1^2*x0^2/2 - 2R^2*r1*x0^2*s + 60R*r1^(7/2)*x0*t - 26R*w*x0*t*s - 2R*v*x0^3*t - 78R*r1^4*s - 36R*r1^3*x0^2 + 5R*r1^2*x0^2*s/2 - 36*r1^(9/2)*x0*t + 20*r1^(7/2)*x0*t*s + w*x0^3*t + 36*r1^5*s + 20*r1^4*x0^2 - r1^3*x0^2*s,  # eq4
       12R^4*r1^2 + 4R^4*x0^2 + 24R^3*v*x0*t - 20R^3*r1^3 - 20R^3*r1^2*s + 8R^3*r1*x0^2 + 4R^3*x0^2*s - 52R^2*w*x0*t - 40R^2*v*x0*t*s + 16R^2*u*x0^3*t + 4R^2*r1^4 + 44R^2*r1^3*s - 73R^2*r1^2*x0^2 - 40R^2*r1*x0^2*s + 24R*r1^(7/2)*x0*t + 76R*w*x0*t*s - 60R*v*x0^3*t - 28R*r1^4*s + 88R*r1^3*x0^2 + 85R*r1^2*x0^2*s - 8*r1^(9/2)*x0*t - 24*r1^(7/2)*x0*t*s + 50*w*x0^3*t + 8*r1^5*s - 24*r1^4*x0^2 - 50*r1^3*x0^2*s,  # eq5
       4R^5*r1 - 4R^4*u*x0*t - 8R^4*u*x2*t - 12R^4*r1^2 + 4R^4*r1*s - 4R^4*r2^2 + R^4*x0^2 + 4R^4*x0*x2 + 4R^4*x2^2 + 24R^3*v*x0*t + 24R^3*v*x2*t - 4R^3*u*x0*t*s - 8R^3*u*x2*t*s + 28R^3*r1^3 - 20R^3*r1^2*s + 8R^3*r1*r2^2 - 6R^3*r1*x0^2 - 20R^3*r1*x0*x2 - 12R^3*r1*x2^2 - 4R^3*r2^2*s + R^3*x0^2*s + 4R^3*x0*x2*s + 4R^3*x2^2*s - 20R^2*w*x0*t - 32R^2*w*x2*t + 8R^2*v*x0*t*s + 24R^2*v*x2*t*s - 8R^2*u*r2^2*x0*t + 6R^2*u*x0^3*t + 12R^2*u*x0^2*x2*t + 4R^2*u*x0*x2^2*t - 20R^2*r1^4 + 20R^2*r1^3*s - 8R^2*r1^2*r2^2 + 21R^2*r1^2*x0^2 + 24R^2*r1^2*x0*x2 + 12R^2*r1^2*x2^2 + 8R^2*r1*r2^2*s - 16R^2*r1*x0^2*s - 20R^2*r1*x0*x2*s - 12R^2*r1*x2^2*s + 2R^2*r2^2*x0^2 - 2R^2*x0^3*x2 - 3R^2*x0^2*x2^2 + 8R*r1^(7/2)*x0*t - 28R*w*x0*t*s - 16R*w*x2*t*s - 24R*v*x0^3*t - 32R*v*x0^2*x2*t - 8R*v*x0*x2^2*t - 8R*u*r2^2*x0*t*s + 8R*u*x0^2*x2*t*s + 4R*u*x0*x2^2*t*s - 12R*r1^4*s - 40R*r1^3*x0^2 - 24R*r1^3*x0*x2 + 33R*r1^2*x0^2*s + 48R*r1^2*x0*x2*s + 12R*r1^2*x2^2*s - 16R*r1*r2^2*x0^2 + 14R*r1*x0^3*x2 + 8R*r1*x0^2*x2^2 - R*x0^2*x2^2*s - 8*r1^(9/2)*x0*t + 24*r1^(7/2)*x0*t*s + 16*r1^(7/2)*x2*t*s + 18*w*x0^3*t + 32*w*x0^2*x2*t + 8*w*x0*x2^2*t + 16*v*r2^2*x0*t*s - 12*v*x0^2*x2*t*s - 8*v*x0*x2^2*t*s - 2*u*x0^3*x2^2*t + 8*r1^5*s + 24*r1^4*x0^2 + 16*r1^4*x0*x2 - 18*r1^3*x0^2*s - 32*r1^3*x0*x2*s - 8*r1^3*x2^2*s + 16*r1^2*r2^2*x0^2 - 12*r1^2*x0^3*x2 - 8*r1^2*x0^2*x2^2 + 2*r1*x0^2*x2^2*s,  # eq6
   ]
end;

r3 = 8/2
iniv = BigFloat[18.8, 5.8, 6.4, 10.8, 9.6, 11.8]
res = nls(H, ini=iniv)

   ([18.77777777777778, 5.777777777777778, 6.5, 10.833333333333334, 9.614803401237305, 11.786666666666667], true)

甲円の直径は乙円の直径の 13/8 倍である。
乙円の直径が 8 寸のとき，甲円の直径は 13 寸である。

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

   R = 169/9, r1 = 52/9, r2 = 13/2, x2 = 65/6, x3 = 9.614803401237305, x0 = 884/75

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 8/2
   (R, r1, r2, x2, x3, x0) = [169/9, 52/9, 13/2, 65/6, 9.614803401237305, 884/75]
   (R, r1, r2, x2, x3, x0) = [18.77777777777778, 5.777777777777778, 6.5, 10.833333333333334, 9.614803401237305, 11.786666666666667]
   y2 = (R - 2r1) - (2R - 2r1)/2
   plot()
   circle(0, 0, R, :blue)
   circle(0, R - r1, r1, :magenta)
   circle(0, R - 3r1, r1, :magenta)
   circle(0, r1 -R, r1, :magenta)
   circle2(x2, y2, r2)
   circle2(x3, R - 2r1 + r3, r3, :orange)
   segment(-sqrt(R^2 - (R - 2r1)^2), R - 2r1, sqrt(R^2 - (R - 2r1)^2), R - 2r1)
   segment(x0, -sqrt(R^2 - x0^2), -sqrt(R^2 - (R - 2r1)^2), R - 2r1)
   segment(-x0, -sqrt(R^2 - x0^2), sqrt(R^2 - (R - 2r1)^2), R - 2r1)
   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, " R", :blue, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(sqrt(R^2 - (R - 2r1)^2), R - 2r1, "(√(R^2-(R-2r1)^2),R-2r1)", :black, :right, delta=-delta/2)
       point(-x0, -sqrt(R^2 - x0^2), "(-x0,-√(R^2-x0^2))", :black, :center, delta=-delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :magenta, :center, delta=-delta/2)
       point(0, R - 3r1, "甲円:r1,(0,R-3r1)", :magenta, :center, delta=-delta/2)
       point(0, r1 - R, "甲円:r1,(0,r1-R)", :magenta, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :red, :center, delta=-delta/2)
       point(x3, R - 2r1 + r3, "丙円:r3\n(x3,R-2r1+r3)", :black, :center, delta=-delta/2)
   end
end;

算額（その1102）

算額（その1102）

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

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

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (R - r2, 0)
丙円の半径と中心座標を r3, (x2, y3)
丁円の半径と中心座標を r4, (r5 + r4, 0)
戊円の半径と中心座標を r5, (0, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, r5::positive
R = 2r1 + r5
r4 = r1 - r2
eq1 = (R - r2)^2 + (R - r1)^2 - (r1 + r2)^2 |> expand
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2 |> expand
eq3 = (r4 + r5)^2 + (R - r1)^2 - (r1 + r4)^2 |> expand
eq4 = (x3 - R + r2)^2 + y3^2 - (r2 + r3)^2 |> expand
eq5 = x3^2 + y3^2 - (R - r3)^2 |> expand;

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   4*r1^2 - 6*r1*r2 + 6*r1*r5 - 2*r2*r5 + 2*r5^2,  # eq1
   -2*r1*r3 + 2*r1*r5 - 2*r1*y3 - r3^2 + r5^2 - 2*r5*y3 + x3^2 + y3^2,  # eq2
   -2*r1^2 + 2*r1*r2 + 4*r1*r5 - 2*r2*r5 + 2*r5^2,  # eq3
   4*r1^2 - 4*r1*r2 + 4*r1*r5 - 4*r1*x3 - 2*r2*r3 - 2*r2*r5 + 2*r2*x3 - r3^2 + r5^2 - 2*r5*x3 + x3^2 + y3^2,  # eq4
   -4*r1^2 + 4*r1*r3 - 4*r1*r5 - r3^2 + 2*r3*r5 - r5^2 + x3^2 + y3^2,  # eq5

res = solve([eq1, eq3], (r1, r2))[1];

ans_r1 = res[1] |> simplify
ans_r1 |> println

   r5*(3 + 2*sqrt(3))

ans_r2 = res[2]|> simplify
ans_r2 |> println

   r5*(5 + 3*sqrt(3))/2

eq12 = eq2(r1 => ans_r1, r2 => ans_r2) |> simplify
eq14 = eq4(r1 => ans_r1, r2 => ans_r2) |> simplify
eq15 = eq5(r1 => ans_r1, r2 => ans_r2) |> simplify
res2 = solve([eq12, eq14, eq15], (r3, x3, y3))[1];

#= r3 =# res2[1] |> sympy.sqrtdenest |> simplify |> println

   r5*(sqrt(3) + 3)/2

#= x3 =# res2[2] |> sympy.sqrtdenest |> simplify |> println

   r5*(5*sqrt(3) + 9)/2

#= y3 =# res2[3] |> sympy.sqrtdenest |> simplify |> println

   2*sqrt(3)*r5 + 4*r5

丙円の半径 r3 は 戊円の半径 r5 の (√3 + 3)/2 倍である。
戊円の直径が 1 寸のとき，丙円の直径は 2.3660254037844384 寸である。

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

  r5 = 0.5;  r1 = 3.23205;  r2 = 2.54904;  r3 = 1.18301;  x3 = 4.41506;  y3 = 3.73205

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