算額（その914）

算額（その914）

一〇八 加須市騎西町 玉敷神社 大正4年(1915)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に甲円，乙円を 2 個ずつ入れ，甲円，乙円に外接する丙円を 2 個描く。甲円，乙円，外円の直径がそれぞれ 9 寸，3 寸，26 寸のとき，丙円の直径はいかほどか。

後述するが，上の図は外円，甲円，乙円の直径が 20, 6, 3 のときのもので，丙円の直径は 12.3435 である。
その他のパラメータは，R = 10;  r1 = 3;  r2 = 1.5;  y1 = -6.32456;  y2 = 8.3666;  r3 = 6.17176;  y3 = 2.28132

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (r3, y3)
とおき，以下の連立方程式を解く。

まとめて解くには SymPy の能力が足りない。

include("julia-source.txt");

using SymPy
@syms R::positive,
     r1::positive, y1::negative,
     r2::positive, y2::positive,
     r3::positive, y3::positive
eq1 = r1^2 + y1^2 - (R - r1)^2 |> expand
eq2 = r2^2 + y2^2 - (R - r2)^2 |> expand
eq3 = (r3 - r1)^2 + (y3 - y1)^2 - (r1 + r3)^2 |> expand
eq4 = (r3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand;

外円，甲円，乙円の直径が与えられるので，eq1, eq2 はそれぞれだけで未知数 y1, y2 を求めることができる。

(ans_y1, ans_y2) = solve([eq1, eq2], (y1, y2))[1]

   (-sqrt(R)*sqrt(R - 2*r1), sqrt(R)*sqrt(R - 2*r2))

y1, y2 が既知として，eq3, eq4 の連立方程式を解く。

(y1, y2) = (-sqrt(R)*sqrt(R - 2*r1), sqrt(R)*sqrt(R - 2*r2))
eq3 = (r3 - r1)^2 + (y3 - y1)^2 - (r1 + r3)^2 |> expand
eq4 = (r3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand;
(ans_r3, ans_y3) = solve([eq3, eq4], (r3, y3))[1]  # 1 of 2

   ((-R*r1 + R*r2 + (sqrt(R)*sqrt(R - 2*r1) + sqrt(R)*sqrt(R - 2*r2))*(-sqrt(2)*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2 + sqrt(R - 2*r1)*sqrt(R - 2*r2))/(r1 - r2) + sqrt(R)*(r1*sqrt(R - 2*r2) + r2*sqrt(R - 2*r1))/(r1 - r2)))/(2*(r1 - r2)), -sqrt(2)*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2 + sqrt(R - 2*r1)*sqrt(R - 2*r2))/(r1 - r2) + sqrt(R)*(r1*sqrt(R - 2*r2) + r2*sqrt(R - 2*r1))/(r1 - r2))

R, r1, r2 を与えれば，それぞれのパラメータの数値解が得られる。

甲円，乙円，外円の直径がそれぞれ 9 寸，3 寸，26 寸のとき，以下のようになる。

ans_y1(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_y2(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_r3(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_y3(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println

   -7.21110255092798
   11.4017542509914
   7.73565831477412
   4.58897582448691

「術」では，「1385.4 を 89.57 で割れば丙円の直径が得られる」とあるが，「甲円，乙円，外円の直径がそれぞれ 9 寸，3 寸，26 寸のとき」の解であるし 1385.4/89.57 = 15.467232332254104 で，「答」の「丙円径十五寸四分二厘」とも小数点以下第2位で不一致である。

丙円の直径の正確な値は上に示すように 2*7.73565831477412 = 15.47131662954824 である。

任意の 甲円，乙円，外円の直径に対しては，上述の数式解に定数を代入すればよい。

y1, y2, r3, y3 の式は長く，SymPy では簡約化できないが，以下のようにすればプログラム的には短くなる。

   s = √(R - 2r1)
   t = √(R - 2r2)
   u = -√(2R*r1*r2*(R - r1 - r2 + s*t))
   v = √R*(r1*t + r2*s)
   w = r1 - r2
   (y1, y2, r3, y3) = (
       -s√R,
       t√R,
       (s + t)*(u + v)√R/(2w^2) - R/2,
       (u + v)/w
   )

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #(R, r1, r2) = (26, 9, 3) ./ 2
   #(y1, y2, r3, y3) = (-7.21110255092798,
   #                    11.4017542509914,
   #                    7.73565831477412,
   #                    4.58897582448691)
   (R, r1, r2) = (20, 6, 3) ./ 2
   s = √(R - 2r1)
   t = √(R - 2r2)
   u = -√(2R*r1*r2*(R - r1 - r2 + s*t))
   v = √R*(r1*t + r2*s)
   w = r1 - r2
   (y1, y2, r3, y3) = (
       -s√R,
       t√R,
       (s + t)*(u + v)√R/(2w^2) - R/2,
       (u + v)/w
   )

   @printf("外円，甲円，乙円の直径が %g, %g, %g のとき，丙円の直径は %g\n", 2R, 2r1, 2r2, 2r3)
   @printf("R = %g;  r1 = %g;  r2 = %g;  y1 = %g;  y2 = %g;  r3 = %g;  y3 = %g\n", R, r1, r2, y1, y2, r3, y3)
   plot()
   circle(0, 0, R, :brown)
   circle2(r1, y1, r1, :green)
   circle2(r2, y2, r2)
   circle2(r3, y3, r3, :blue)
   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", :green, :left, :bottom, delta=delta/2)
       point(r1, y1, "甲円:r1,(r1,y1)", :green, :center, delta=-delta/2)
       point(r2, y2, "乙円:r2,(r2,y2)", :black, :center, :bottom, delta=delta/2)
       point(r3, y3, "丙円:r3,(r3,y3)", :blue, :center, delta=-delta/2)
   end
end;

算額（その913）

算額（その913）

一〇八 加須市騎西町 玉敷神社 大正4年(1915)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に正方形 3 個，甲円 3 個，乙円 6 個を入れる。甲円の直径が 1 寸のとき，乙円の直径はいかほどか。

正方形の一辺の長さを a
円周上にある正方形の頂点の座標を (x0, y0)
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive, r1::positive,
     r2::positive, x2::positive, y2::positive,
     x0::positive, y0::positive
x0 = √Sym(3)a
y0 = (2√Sym(3)/3 + 1)a

eq1 = dist2(0, 2a√Sym(3)/3, x0, y0, 0, R - r1, r1)
eq2 = dist2(0, 2a√Sym(3)/3, x0, y0, x2, y2, r2)
eq3 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq4 = x2^2 + (y2 - R + r1)^2 - (r1 + r2)^2 |> expand
eq5 = x0^2 + y0^2 - R^2 |> simplify;

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, R, r2, x2, y2) = u
   return [
       3*R^2/4 - sqrt(3)*R*a - 3*R*r1/2 + a^2 + sqrt(3)*a*r1 - r1^2/4,  # eq1
       a^2 + a*x2 - sqrt(3)*a*y2 - r2^2 + x2^2/4 - sqrt(3)*x2*y2/2 + 3*y2^2/4,  # eq2
       x2^2 + y2^2 - (R - r2)^2,  # eq3
       x2^2 - (r1 + r2)^2 + (-R + r1 + y2)^2,  # eq4
       -R^2 + 3*a^2 + a^2*(1 + 2*sqrt(3)/3)^2,  # eq5
   ]
end;

r1 = 0.5
iniv = BigFloat[0.7, 1.9, 0.2, 0.7, 1.5]
res = nls(H, ini=iniv)

   ([0.6692244041071277, 1.8501040489086598, 0.2283119839009551, 0.7210516428693637, 1.4526851105581506], true)

甲円の直径が 1 寸のとき，乙円の直径は 2*0.2283119839009551 = 0.4566239678019102 寸である。

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

   r1 = 0.5;  a = 0.669224;  R = 1.8501;  r2 = 0.228312;  x2 = 0.721052;  y2 = 1.45269

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 0.5
   (a, R, r2, x2, y2) = res[1]
   @printf("甲円の直径が %g のとき，乙円の直径は %g\n", 2r1, 2r2)
   @printf("r1 = %g;  a = %g;  R = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", r1, a, R, r2, x2, y2)
   x0 = √3a
   y0 = (2√3/3 + 1)a
   plot()
   circle(0, 0, R, :brown)
   rotate(0, R - r1, r1, :green)
   rotate(x2, y2, r2)
   rotate(-x2, y2, r2)
   x = a.*[0, √3, 1 + √3, 1, 0]
   y = a.*[2/√3, (3 + 2√3)/3, (3 - √3)/3, -√3/3, 2/√3]
   for i = 1:3
       plot!(x, y, color=:blue, lw=0.5)
       (x, y) = transform(x, y, deg=120)
   end
   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", :green, :left, :bottom, delta=delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :green, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :black, :center, :bottom, delta=delta/2)
       point(0, 2a/√3, " (0,2a/√3)", :blue, :left, :vcenter, delta=-delta/2)
       point(x0, y0, " (x0,y0)", :blue, :left, :vcenter)
   end
end;

算額（その912）

算額（その912）

一〇八 加須市騎西町 玉敷神社 大正4年(1915)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に 2 本の水平な弦を引き，甲円，乙円，丙円をそれぞれ 2 個，丁円を 1 個入れる。丁円の直径が 1 寸のとき，甲円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive, r4::positive
eq1 = r4 + 2r1 - R
eq2 = x2^2 + (R - 2r4 - 2r1 + r2)^2 - (R - r2)^2
eq3 = x3^2 + (R - 2r4 - r3)^2 - (R - r3)^2
eq4 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq6 = (x2 - x3)^2 + (2r1 - r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, r1, r2, x2, r3, x3))[1]

   (9*r4, 4*r4, 5*r4/2, 2*sqrt(10)*r4, 8*r4/5, 8*sqrt(10)*r4/5)

甲円の半径は丁円の半径の 4 倍である。
丁円の直径が 1 寸ならば，甲円の直径は 4 寸である。

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

R = 4.5;  r1 = 2;  r2 = 1.25;  x2 = 3.16228;  r3 = 0.8;  x3 = 2.52982

function draw(more=false)
   pyplot(size=(600, 600), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r4 = 1/2
   (R, r1, r2, x2, r3, x3) = r4 .* (9, 4, 5/2, 2√10, 8/5, 8√10/5)
   @printf("丁円の直径が %g 寸のとき，甲円の直径は %g である。\n", 2r4, 2r1)
   @printf("r4 = %g;  R = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r4, R, r1, r2, x2, r3, x3)
   y01 = R - 2r4
   x01 = sqrt(R^2 - y01^2)
   y02 = y01 - 2r1
   x02 = sqrt(R^2 - y02^2)
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r4, r4, :magenta)
   circle(0, y01 - r1, r1)
   circle(0, y02 - r1, r1)
   circle2(x2, y02 + r2, r2, :orange)
   circle2(x3, y01 - r3, r3, :brown)
   segment(-x01, y01, x01, y01)
   segment(-x02, y02, x02, y02)
   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(R, 0, " R", :green, :left, :bottom, delta=delta/2)
       point(0, R, " R", :green, :left, :bottom, delta=delta/2)
       point(0, R - 2r4 - r1, "甲円:r1,(0,R-2r4-r1)", :red, :center, delta=-delta/2)
       point(0, R - 2r4 - 3r1, "甲円:r1,(0,R-2r4-3r1)", :red, :center, delta=-delta/2)
       point(x2, R - 2r4 - 2r1 + r2, "乙円:r2\n(x2,R-2r4-2r1+r2)", :orange, :center, :bottom, delta=delta/2)
       point(x3, R - 2r4 - r3, "乙円:r2\n(x3,R-2r4-r3)", :brown, :center, :bottom, delta=delta/2)
       point(0, R - r4, " 丙円:r4,(0,R-r4)", :black, :left, :vcenter)
       point(0, R - 2r4, "R-2r4", :green, :center, delta=-delta/2)
       point(0, R - 2r4 - 2r1, "R-2r4-2r1", :green, :center, delta=-delta/2)
   end
end;

算額（その911）

算額（その911）

一一三 熊谷市熊谷 大神宮 昭和辛丑（昭和36年，1961年）

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

楕円内に正方形 2 個，等円 2 個を入れる。等円の直径が与えられたとき，正方形の一辺の長さを求めよ。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
等円の半径と中心座標を r, (0, b - r)
とおき，以下の連立方程式を解く。

なお，正方形の一辺の長さは楕円の長半径の 1/√2 倍である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r::positive
eq1 = (b - r)/√Sym(2) - r
eq2 = (a/2)^2/a^2 + (a/2)^2/b^2 - 1
res = solve([eq1, eq2], (a, b))[1]

   (r*sqrt(6*sqrt(2) + 9), r*(1 + sqrt(2)))

res[1] |> sympy.sqrtdenest |> factor |> println

   r*(sqrt(3) + sqrt(6))

楕円の長半径は等円の半径の √3 + √6 倍，楕円の短半径は等円の半径の 1 + √2 倍である。
等円の直径が 1 のとき，楕円の長径，短径はそれぞれ 4.18154055035206，2.41421356237309 である。

正方形の一辺の長さは楕円の長半径の 1/√2 倍であるから，等円の半径の √2(√3 + √6)/2 倍である。
等円の直径が 1 のとき，正方形の一辺の長さは 2.9567956789604666 である。

√2(√3 + √6)/2

   2.9567956789604666

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (sq, a, b) = r .* (√2(√3 + √6), √3(1 + √2) , 1 + √2)
   @printf("正方形の一辺の長さ = %g;  楕円の長径 = %g，楕円の短径 = %g\n", sq, 2a, 2b)
   plot(a .* [1, 1/2, 0, -1/2, -1, -1/2, 0, 1/2, 1],
        a/2 .* [0, 1, 0, -1, 0, 1, 0, -1, 0], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   circle22(0, b - r, r, :green)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a, 0, " a", :red, :left, :bottom, delta=delta/2)
       point(0, b, " b", :red, :left, :bottom, delta=delta/2)
       point(0, b - r, "等円:r\n(0,b-r)", :green, :center, delta=-2delta)
       point(a/2, a/2, "(a/2,a/2))", :blue, :left, :bottom, delta=delta)
   end
end;