裏 RjpWiki

算額（その885）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

外円の中に，大円 2 個，小円 6 個を入れる。大円の直径が 2 寸のとき，小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1); r1 = R/2
小円の半径と中心座標を r2, (x2, 2r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive
R = 2r1
eq1 = x2^2 + 4r2^2 - (R - r2)^2
eq2 = x2^2 + (r1 - 2r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r2, x2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (2*r1/5, 4*sqrt(3)*r1/5)

小円の半径は大円の半径の 2/5 倍である。
大円の直径が 2 寸のとき，小円の直径は 4/5 寸 = 8分 である。

その他のパラメータは以下のとおりである。
   r1 = 1;  r2 = 0.4;  x2 = 1.38564;  R = 2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 2/2
   (r2, x2) = r1 .* (2/5, 4√3/5)
   R = 2r1
   @printf("大円の直径が %g のとき，小円の直径は %g である\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  R = %g\n", r1, r2, x2, R)
   plot()
   circle(0, 0, R, :blue)
   circle22(0, r1, r1)
   circle4(x2, 2r2, r2, :green)
   circle2(x2, 0, 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(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, "大円:r1\n(0,r1)", :red, :center, delta=-delta/2)
       point(0, -r1, "大円:r1\n(0,-r1)", :red, :center, delta=-delta/2)
       point(0, 0, "", :blue)
       point(x2, 2r2, "小円:r2\n(x2,2r2)", :green, :center, delta=-delta)
   end
end;

算額（その884）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

外円内に二等辺三角形と，大円を 2 個，小円を 2 個入れる。小円の直径が 1 寸のとき，大円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, -r1), (0, r1); r1 = R/2
小円の半径と中心座標を r2, (x2, y2); y2 < 0
二等辺三角形の底辺と外円の接点座標を (x0, y0); y0 < 0
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, x0::positive, y0::negative, r1::positive,
     r2::positive, x2::positive, y2::negative
R = 2r1
y0 = y2 + r2
L = sqrt(x0^2 + (2r1 - y0)^2)
eq1 = x0^2 + y0^2 - R^2
eq2 = r1/3r1 - x0/L
eq3 = x2^2 + (r1 + y2)^2 - (r1 + r2)^2
eq4 = x2^2 + y2^2 - (R - r2)^2;
res = solve([eq1, eq2, eq3, eq4], (r1, x2, y2, x0))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (9*r2, 8*r2, -15*r2, 8*sqrt(2)*r2)

大円の半径は小円の半径の 9 倍である。
小円の直径が 1 寸のとき，大円の直径は 9 寸である。

その他のパラメータは以下のとおりである。
r2 = 0.5;  r1 = 4.5;  x2 = 4;  y2 = -7.5;  R = 9;  x0 = 5.65685;  y0 = -7

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, x2, y2, x0) = r2 .* (9, 8, -15, 8√2)
   R = 2r1
   y0 = y2 + r2
   @printf("小円の直径が %g のとき，大円の直径は %g である\n", 2r2, 2r1)
   @printf("r2 = %g;  r1 = %g;  x2 = %g;  y2 = %g;  R = %g;  x0 = %g;  y0 = %g\n", r2, r1, x2, y2, R, x0, y0)
   plot()
   circle(0, 0, R, :blue)
   circle22(0, r1, r1)
   plot!([x0, 0, -x0, x0], [y0, R, y0, y0], color=:green, lw=0.5)
   circle2(x2, y2, r2, :magenta)
   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(0, r1, "大円:r1\n(0,r1)", :red, :center, delta=-delta/2)
       point(0, -r1, "大円:r1\n(0,-r1)", :red, :center, delta=-delta/2)
       point(0, 0, "", :blue)
       point(x0, y0, " (x0,y0)", :blue, :left, :vcenter)
       point(x2, y2, "  小円:r2,(x2,y2)", :magenta, :left, delta=-delta)
   end
end;

算額（その883）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

半円の中に二等辺三角形と等円 3 個を入れる。等円の直径が 1 寸のとき，半円の直径はいかほどか。

二等辺三角形の底辺の長さを 2a
半円の半径と中心座標を r1, (0, 0)
等円の半径と中心座標を r2, (0, r2), (x2, r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms d, a::positive, r1::positive, r2::positive,
     x2::positive
eq1 = numerator(apart(dist(a, 0, 0, r1, 0, r2) - r2^2, d))
eq2 = numerator(apart(dist(a, 0, 0, r1, x2, r2) - r2^2, d))
eq3 = x2^2 + r2^2 - (r1 - r2)^2
res = solve([eq1, eq2, eq3], (a, r1, x2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r2*sqrt(sqrt(2) + 2 + sqrt(2*sqrt(2) + 3))/sqrt(1 + sqrt(2)), r2*(1 + sqrt(2*sqrt(2) + 3)), r2*sqrt(2 + 2*sqrt(2)))

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

   r2*(sqrt(2) + 2)

半円の半径は等円の半径の 1 + sqrt(2√2 + 3) = √2 + 2 倍である。
等円の直径が 1 寸のとき，半円の直径は 3.414213562373095 寸である。

算額では「術曰置二個開平方加二個乗等円径得半円径合問」は，「等円の直径の √2 + 2 倍」である。

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

r2 = 0.5;  a = 0.776887;  r1 = 1.70711;  x2 = 1.09868

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, r1, x2) = (r2*sqrt(sqrt(2) + 2 + sqrt(2*sqrt(2) + 3))/sqrt(1 + sqrt(2)), r2*(1 + sqrt(2*sqrt(2) + 3)), r2*sqrt(2 + 2*sqrt(2)))
   @printf("等円の直径が %g のとき，半円の直径は %g である\n", 2r2, 2r1)
   @printf("r2 = %g;  a = %g;  r1 = %g;  x2 = %g\n", r2, a, r1, x2)
   plot(ylims=(-0.15, 1.85))
   circle(0, 0, r1, beginangle=0, endangle=180)
   circle(0, r2, r2, :blue)
   circle2(x2, r2, r2, :blue)
   segment(0, r1, a, 0, :green)
   segment(0, r1, -a, 0, :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(a, 0, "a", :green, :center, delta=-delta)
       point(0, r1, "r1", :red, :center, :bottom, delta=delta)
       point(r1, 0, "r1", :red, :center, delta=-delta)
       point(0, r2, "等円:r2,(0,r2)", :blue, :center, delta=-delta)
       point(x2, r2, "等円:r2,(x2,r2)", :blue, :center, delta=-delta)
       point(0, 0, "", :red)
   end
end;

算額（その882）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

長方形の中に対角線と斜線を引き，区画された領域に甲円，乙円 2 個ずつを入れる。乙円の直径が 1 寸のとき，長方形の長辺の長さはいかほどか。

長方形の長辺，短辺をそれぞれ a, b
斜線と長辺の交点座標を (c, 0)
甲円の半径と中心座標を r1, (a - r1, y1), (a/2, r1)
乙円の半径と中心座標を r2, (r2, b/2 + r1)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms d, a::positive, b::positive, c::positive,
     r1::positive, y1::positive, r2::positive
eq1 = numerator(apart(dist(a, 0, 0, b, a - r1, y1) - r1^2, d))
eq2 = numerator(apart(dist(a, 0, 0, b, a/2, r1) - r1^2, d))
eq3 = numerator(apart(dist(a, 0, 0, b, r2, b/2 + r2) - r2^2, d))
eq4 = numerator(apart(dist(a, b, c, 0, a - r1, y1) - r1^2, d))
eq5 = numerator(apart(dist(a, b, c, 0, a/2, r1) - r1^2, d));

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, b, c, r1, y1) = u
   return [
       -a^2*r1^2 + a^2*y1^2 - 2*a*b*r1*y1,  # eq1
       a^2*b^2 - 4*a^2*b*r1 - 4*b^2*r1^2,  # eq2
       a^2*b^2 - 4*a^2*b*r2 - 4*a*b^2*r2 + 8*a*b*r2^2,  # eq3
       a^2*b^2 - 2*a^2*b*y1 - a^2*r1^2 + a^2*y1^2 - 2*a*b^2*c - 2*a*b^2*r1 + 4*a*b*c*y1 + 2*a*b*r1*y1 + 2*a*c*r1^2 - 2*a*c*y1^2 + b^2*c^2 + 2*b^2*c*r1 - 2*b*c^2*y1 - 2*b*c*r1*y1 - c^2*r1^2 + c^2*y1^2,  # eq4
       a^2*b^2 - 4*a^2*b*r1 - 4*a*b^2*c + 12*a*b*c*r1 + 4*b^2*c^2 - 4*b^2*r1^2 - 8*b*c^2*r1,  # eq5
   ]
end;

r2 = 1/2
iniv = BigFloat[53, 28, 30, 6, 11] ./5
res = nls(H, ini=iniv)

   ([5.23606797749979, 2.618033988749895, 2.8944271909999157, 0.6180339887498949, 1.0], true)

長方形の長辺の長さは 5.23606797749979 である。

算額では「術曰置五個開平方加三個得直長合問」は √5 + 3 = 5.23606797749979 なので，数値解がこれに一致していることがわかる。

その他のパラメータは以下のとおりである。
r2 = 0.5;  a = 5.23607;  b = 2.61803;  c = 2.89443;  r1 = 0.618034;  y1 = 1

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, c, r1, y1) = res[1]
   @printf("乙円の直径が %g のとき，長方形の長辺は %g である\n", 2r2, a)
   @printf("r2 = %g;  a = %g;  b = %g;  c = %g;  r1 = %g;  y1 = %g\n", r2, a, b, c, r1, y1)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:magenta, lw=0.5)
   circle(a - r1, y1, r1)
   circle(a/2, r1, r1)
   circle(r2, b/2 + r2, r2, :green)
   circle(r2, b/2 - r2, r2, :green)
   segment(0, 0, a, b, :blue)
   segment(0, b, a, 0, :blue)
   segment(c, 0, a, b, :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(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
       point(c, 0, "  c", :magenta, :left, :bottom, delta=delta/2)
       point(0, b, " b", :magenta, :left, :bottom, delta=delta/2)
       point(a - r1, y1, "甲円:r1\n(a-r1,y1)", :red, :center, delta=-delta/2)
       point(a/2, r1, "甲円:r1\n(a/2,r1)", :red, :center, delta=-delta/2)
       point(r2, b/2 + r2, "乙円:r2\n(r2,b/2+r2)", :green, :center, :bottom, delta=delta/2)
   end
end;

算額（その881）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

直線の上に乙円 4 個が載り，その上に甲円，丙円，丁円が 2 個ずつ載っている。乙円の直径が 2 寸のとき，甲円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms d, r1::positive, r1::positive, y1::positive,
     r2::positive, r3::positive, y3::positive, r4::positive
eq1 = r2^2 + (y1 - r2)^2 - (r1 + r2)^2
eq2 = 4r2^2 + (y1 - y3)^2 - (r1 + r3)^2
eq3 = r2^2 + (y3 - r2)^2 - (r2 + r3)^2
eq4 = r4^2 + (y1 - y3)^2 - (r3 + r4)^2
eq5 = r1 + 2r4 - 2r2
solve([eq1, eq2, eq3, eq4, eq5], (r1, y1, r3, y3, r4))

   1-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (3*r2/2, r2*(2 + sqrt(21))/2, 7*r2/10, r2*(1 + 3*sqrt(21)/10), r2/4)

甲円の半径は乙円の半径の 3/2 倍である。
乙円の直径が 2 寸のとき，甲円の直径は 3 寸である。

その他のパラメータは以下のとおりである。
r1 = 1.5;  y1 = 3.29129;  r2 = 1;  r3 = 0.7;  y3 = 2.37477;  r4 = 0.25

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

算額（その880）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

二等辺三角形内に，甲円 1 個，乙円 4 個，丙円 2 個が入っている。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。

二等辺三角形の底辺と高さを 2a，b
甲円の半径と中心座標を r1, (0, y1)
乙円の半径と中心座標を r2, (r2, r2), (3r2, r2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms d, a::positive, b::positive,
     r1::positive, y1::positive,
     r2::positive, r3::positive, y3::positive
eq1 = numerator(apart(dist(a, 0, 0, b, 0, y1) - r1^2, d))
eq2 = numerator(apart(dist(a, 0, 0, b, 3r2, r2) - r2^2, d))
eq3 = numerator(apart(dist(a, 0, 0, b, 2r2, y3) - r3^2, d))
eq4 = r2^2 + (y1 - r2)^2 - (r1 + r2)^2 |> expand
eq5 = 4r2^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq6 = r2^2 + (y3 - r2)^2 - (r2 + r3)^2 |> expand;

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, b, r1, y1, r3, y3) = u
   return [
       a^2*b^2 - 2*a^2*b*y1 - a^2*r1^2 + a^2*y1^2 - b^2*r1^2,  # eq1
       a^2*b^2 - 2*a^2*b*r2 - 6*a*b^2*r2 + 6*a*b*r2^2 + 8*b^2*r2^2,  # eq2
       a^2*b^2 - 2*a^2*b*y3 - a^2*r3^2 + a^2*y3^2 - 4*a*b^2*r2 + 4*a*b*r2*y3 + 4*b^2*r2^2 - b^2*r3^2,  # eq3
       -r1^2 - 2*r1*r2 + r2^2 - 2*r2*y1 + y1^2,  # eq4
       -r1^2 - 2*r1*r3 + 4*r2^2 - r3^2 + y1^2 - 2*y1*y3 + y3^2,  # eq5
       r2^2 - 2*r2*r3 - 2*r2*y3 - r3^2 + y3^2,  # eq6
   ]
end;

r2 = 1/2
iniv = BigFloat[3, 3, 1, 1.5, 1, 2]
res = nls(H, ini=iniv)

   ([2.726702475495804, 2.665648262572091, 0.73841681234051, 1.6329943517456873, 0.3542486889354094, 1.1926332525571277], true)

乙円の直径が 1 のとき，甲円の直径は 1.47683362468102 である。

算額では「術曰置二十八個開平方加八個除以球個得甲円径合問」で，数式で表すと (sqrt(28) + 8)/9 = 1.47683362468102 となり，数値解が正しいことがわかる。

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

a = 2.7267;  b = 2.66565;  r1 = 0.738417;  y1 = 1.63299;  r3 = 0.354249;  y3 = 1.19263

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, r1, y1, r3, y3) = res[1]    
   @printf("乙円の直径が %g のとき，甲円の直径は %.15g である\n", 2r2, 2r1)
   @printf("a = %g;  b = %g;  r1 = %g;  y1 = %g;  r3 = %g;  y3 = %g\n", a, b, r1, y1, r3, y3)
   plot([a, 0, -a, a], [0, b, 0, 0], color=:green, lw=0.5)
   circle(0, y1, r1, :magenta)
   circle2(r2, r2, r2)
   circle2(3r2, r2, r2)
   circle2(2r2, 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, y1, "甲円:r1,(0, y1)", :magenta, :center, delta=-delta)
       point(r2, r2, "乙円:r2\n(r2,r2)", :red, :center, delta=-delta)
       point(3r2, r2, "乙円:r2\n(3r2,r2)", :red, :center, delta=-delta)
       point(2r2, y3, " 丙円:r3,(2r2,y3)", :blue, :left, :vcenter)
       point(a, 0, "a", :green, :left, :bottom, delta=delta)
       point(0, b, " b", :green, :left, :vcenter)
   end
end;

算額（その879）

六十三 羽生市須影 八幡神社 慶応元年(1865)

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

正方形の中に，四分円が 2 個，半円が 1 個，甲円が 3 個，乙円が 2 個入っている。乙円の直径が 3 寸のとき，甲円の直径はいかほどか。

四分円の半径と中心座標を 2r1, (0, 0)
半円の半径と中心座標を r1, (r1, 0)
甲円の半径と中心座標を r2, (r1, r1 + r2), (2r1 - r2, y2)
乙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive
eq1 = r1^2 + (r1 + r2)^2 - (2r1 - r2)^2
eq2 = (x3 - r1)^2 + y3^2 - (r1 + r3)^2
eq3 = (2r1 - r2)^2 + y2^2 - (2r1 + r2)^2
eq4 = x3^2 + y3^2 - (2r1 - r3)^2
eq5 = (x3 - r1)^2 + (r1 + r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, y2, x3, y3))

   1-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (11*r3/2, 11*r3/6, 11*sqrt(6)*r3/3, 8*r3, 6*r3)

甲円の半径は乙円の半径の 11/6 倍である。
乙円の直径が 3 寸のとき，甲円の直径は 11/2 = 5.5 寸である。

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

r1 = 8.25;  r2 = 2.75;  y2 = 13.4722;  x3 = 12;  y3 = 9

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 3/2
   (r1, r2, y2, x3, y3) = r3 .* [11/2, 11/6, 11√6/3, 8, 6]
   @printf("乙円の直径が %g のとき，甲円の直径は %g である\n", 2r3, 2r2)
   @printf("r1 = %g;  r2 = %g;  y2 = %g;  x3 = %g;  y3 = %g\n", r1, r2, y2, x3, y3)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:green, lw=0.5)
   circle(0, 0, 2r1, :magenta, beginangle=0, endangle=90)
   circle(2r1, 0, 2r1, :magenta, beginangle=90, endangle=180)
   circle(r1, 0, r1, :green, beginangle=0, endangle=180)
   circle(r1, r1 + r2, r2)
   circle(2r1 - r2, y2, r2)
   circle(r2, y2, r2)
   circle(x3, y3, r3, :blue)
   circle(2r1 - x3, 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(r1, r1 + r2, "甲円:r2,(r1,r1+r2)", :red, :center, delta=-delta/2)
       point(2r1 - r2, y2, "甲円:r2,(2r1-r2,y2)", :red, :center, delta=-delta/2)
       point(x3, y3, "乙円:r3\n(x3,y3)", :red, :center, delta=-delta/2)
       point(r1, 0, "r1", :green, :center, :bottom, delta=delta/2)
       point(2r1, 0, "2r1 ", :green, :right, :bottom, delta=delta/2)
   end
end;

算額（その878）

六十三 羽生市須影 八幡神社 慶応元年(1865)

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

正方形内に半円 2 個，等円 2 個を入れる。等円の直径が 1 寸のとき，正方形の一辺の長さはいかほどか。

正方形の一辺の長さを 2r1
半円の半径と中心座標を r1, (r1, 0)
等円の半径と中心座標を r2, (r2, y2), (2r1 - y2, 2r1 - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive
eq1 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq2 = (2r1 - y2 - r2)^2 + (2r1 - r2 - y2)^2 - 4r2^2;
res = solve([eq1, eq2], (r1, y2))

   3-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r2/2, sqrt(2)*r2)
    (r2*(2 - sqrt(2))^2/4, r2*(2 - sqrt(2)))
    (r2*(sqrt(2) + 2)^2/4, r2*(sqrt(2) + 2))

3 組の解が得られるが，3 番目のものが適解である。

res[3][1]/r2 |> expand |> println

   sqrt(2) + 3/2

半円の半径は等円の半径の (√2 + 3/2) 倍である。
等円の直径が 1 寸のとき，半円の直径は 2.914213562373095 寸で，それは正方形の一辺の長さでもある。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, y2) = (r2*(sqrt(2) + 2)^2/4, r2*(sqrt(2) + 2))
   @printf("等円の直径 = %g;  正方形の一辺の長さ = %g\n", 2r2, 2r1)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:green, lw=0.5)
   circle(r1, 0, r1, beginangle=0, endangle=180)
   circle(2r1, r1, r1, beginangle=90, endangle=270)
   circle(r2, y2, r2, :blue)
   circle(2r1 - y2, 2r1 - r2, r2, :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(r1, 0, "r1", :red, :center, :bottom, delta=delta/2)
       point(2r1, 0, " 2r1", :red, :left, :bottom, delta=delta/2)
       point(0, 2r1, " 2r1", :red, :left, :bottom, delta=delta/2)
       point(2r1, r1, "", :red)
       point(r2, y2, "等円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(2r1 - y2, 2r1 - r2, "等円:r2\n(2r1-y2,2r1-r2)", :blue, :center, delta=-delta/2)
   end
end;

算額（その877）

七 川越市石田本郷折戸 地蔵堂 文化元甲子歳

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

外円の中に甲方（正方形）と乙円が入っている。正方形の一辺の長さは乙円の直径より 25 寸短く，正方形の下辺が作る円弧の矢が 5 寸である。このとき，外円の直径はいかほどか。

正方形の一辺の長さを 2a
乙円の半径と中心座標を r, (0, R - r)
乙円と正方形の一辺の長さとの差をそのまま「差」
矢をそのまま「矢」
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, 差::positive, 矢::positive, r::positive,
     a::positive
r  = a + 差//2
eq1 = (矢 - R)^2 + a^2 - R^2
eq2 = 2r + 2a + 矢 - 2R;
res = solve([eq1, eq2], (R, a));
res[2]  # 2 of 2

   (差/2 + 2*sqrt(矢)*sqrt(差 + 4*矢) + 9*矢/2, sqrt(矢)*sqrt(差 + 4*矢) + 2*矢)

2 組の解が得られるが，2 番目のものが適解である。

与えられた条件のもとでは，外円の直径は 130 寸である。

ちなみに，乙円の直径は 75 寸，正方形の一辺の長さは 50 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   差 = 25
   矢 = 5
   (R, a) = (差/2 + 2sqrt(矢*差 + 4矢^2) + 9矢/2, sqrt(矢*差 + 4矢^2) + 2矢)
   r  = a + 差/2
   @printf("乙円の直径 = %g;  矢 = %g;  外円の直径 = %g;  正方形の一辺の長さ = %g\n", 2r, 矢, 2R, 2a)
   plot()
   circle(0, 0, R)
   plot!([a, a, -a, -a, a], (矢 - R) .+ [0, 2a, 2a, 0, 0], color=:blue, lw=0.5)
   circle(0, R - r, r, :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, "R", :red, :center, :bottom, delta=delta/2)
       point(0, R - r, "乙円:r,(0,R-r)", :green, :center, delta=-delta/2)
       point(0, R - 2r, "R-2r=2a+矢-R", :blue, :center, delta=-delta/2)
       point(0, 矢-R, "矢-R", :blue, :center, :bottom, delta=delta/2)
   end
end;

算額（その876）

新潟県長岡市 蒼柴神社 亨和元年(1801)3月
http://www.wasan.jp/niigata/aoshi.html

涌田和芳，外川一仁: 長岡蒼柴神社の算額，長岡工業高等専門学校研究紀要，第42巻，第2号(2006)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_41-45/vol_42_2/42_2_1wakuta.pdf

外円内に，弦と 2 本の斜線を引き，区分された領域に甲円と乙円を 2 個ずつ入れる。大円の直径が 521 寸で，弦と矢の長さの差を最大にするとき，乙円の直径を求めよ。
注：矢（し）とは，弓形の孤の中点から弦におろした垂線
外円の半径と中心座標を R, (0, 0)
弦，矢の長さおよびその差を（そのまま）弦，矢，弦矢差
甲円の半径と中心座標を r1, (0, R - r1), (0, r - 矢 + r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき以下の連立方程式を解く。

まず，弦と矢の長さの差が最大になるときの矢と弦を決定する。

include("julia-source.txt");

using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
     r1::positive, r2::positive, x2::positive, y2::positive, d
R = 521//2
弦 = 2sqrt(R^2 - (R - 矢)^2)
弦矢差 = 弦 - 矢;

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(弦矢差, xlims=(9, R), xlabel="矢", ylabel="弦 - 矢")

矢が 150 前後のときに弦矢差が最大になることがわかる。

最大値を取るときの矢の値を求めるには，導関数を求め，導関数が 0 になるときの値を求めればよい。

g = diff(弦矢差, 矢);
g |> println

   2*(521/2 - 矢)/sqrt(271441/4 - (521/2 - 矢)^2) - 1

mx_a = solve(g, 矢)[1]
mx_a |> factor |> println
mx_a.evalf() |> println

   -521*(-5 + sqrt(5))/10
   144.000858372261

弦矢差(矢 => mx_a.evalf()) |> println

   321.995708138695

弦(矢 => mx_a.evalf()) |> println

   465.996566510956

矢が 521/2 - 521*sqrt(5)/10 = 144.000858372261 のとき，弦と矢の差が最大値 321.995708138695 になる。ちなみにそのときの弦は 465.996566510956 である。

465.996566510956 - 144.000858372261, 321.995708138695

   (321.99570813869497, 321.995708138695)

弦と矢の差が最大となるときの矢，弦が決まったので，その状況における甲円と乙円の半径を求める。

using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
     r1::positive, r2::positive, x2::positive, y2::positive, d
矢 = (5 - sqrt(Sym(5)))R/5
弦 = 2sqrt(R^2 - (R - 矢)^2)
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - r1) - r1^2
eq3 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - 矢 + r1) - r1^2
eq4 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, x2, y2) - r2^2
eq5 = dist(-x, sqrt(R^2 - x^2), 弦/2, R - 矢, x2, y2) - r2^2;

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

R = 521//2
iniv = BigFloat[35, 37, 122, 188, 127]
res = nls(H, ini=iniv)

   ([35.179523802100135, 36.00021459306524, 121.93467698217249, 188.49957081386952, 126.65410684822777], true)

外円の直径が 521 寸のとき，乙円の直径は 72 寸有奇である。

その他のパラメータは以下のとおりである。
r1 = 35.1795;  r2 = 36.0002;  x2 = 121.935;  y2 = 188.5;  x = 126.654

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 521//2
   矢 = (5 - √5)R/5
   弦 = 2sqrt(R^2 - (R - 矢)^2)
   (r1, r2, x2, y2, x) = res[1]
   @printf("弦 = %g;  矢 = %g;  弦と矢の差 = %g\n", 弦, 矢, 弦 - 矢)
   @printf("乙円の直径 = %.15g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x = %g\n", r1, r2, x2, y2, x)
   plot()
   circle(0, 0, R)
   segment(-弦/2, R - 矢, 弦/2, R - 矢, :blue)
   circle2(x2, y2, r2, :green)
   circle(0, R - r1, r1, :magenta)
   circle(0, R - 矢 + r1, r1, :magenta)
   y = sqrt(R^2 - x^2)
   segment(x, y, -弦/2, R - 矢)
   segment(-x, y, 弦/2, R - 矢)
   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, :center, delta=-delta/2)
       point(x, y, "(x,y)", :green, :left, :bottom, delta=delta/2)
       point(弦/2, R - 矢, "弦/2 ", :blue, :right, delta=-delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :black, :center, :bottom, delta=delta/2)
       point(0, R - 矢/2, "R-矢/2   ", :black, :right, :vcenter)
       point(0, R - 矢 + r1, "甲円:r1,(0,R-矢+r1)", :black, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :black, :center, delta=-delta/2)
       point(0, R, "R", :red, :center, :bottom, delta=delta/2)
 end
end;

算額（その875）

新潟県長岡市 蒼柴神社 亨和元年(1801)3月
http://www.wasan.jp/niigata/aoshi.html

涌田和芳，外川一仁: 長岡蒼柴神社の算額，長岡工業高等専門学校研究紀要，第42巻，第2号(2006)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_41-45/vol_42_2/42_2_1wakuta.pdf

角錐台の中に大球 1 個，小球 4 個が入っている。上底，下底の正方形の一辺の長さがそれぞれ 2 寸，6 寸のとき，大球の直径はいかほどか。

角錐と角錐台の高さを h, h2
上底，下底の正方形の一辺の長さを 2a, 2b
大球の半径と中心座標を r1, (0, 0, h2 - r1)
小球の半径と中心座標を r2, (r2, r2, r2)
と置く。

eq1: x 軸 に対して 45° 方向の負の方向から投影すると，大円と小円は互いに外接している。

eq2, eq3, eq4: x 軸の正の方向から y-z 平面に投影された図形は下図のように大円，小円は角錐台の面に内接している。

これらの連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive, h::positive, h2::positive, d
(a, b) = (1, 3)
eq1 = 2r2^2 + (h2 - r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = (b + h - 2r2)^2 - (h^2 + b^2)  # 3 + h - sqrt(h^2 + 9) - 2r2
eq3 = b/h - a/(h - h2)
eq4 = numerator(apart(dist(0, h, b, 0, 0, h2 - r1) - r1^2, d))
res = solve([eq1, eq2, eq3, eq4], (r1, r2, h, h2));
res[1]  # 1 of 3

   (3/2, 6/5, 36/5, 24/5)

3 組の解が得られるが，最初のものが適解である。

上底，下底の正方形の一辺の長さがそれぞれ 2 寸，6 寸のとき，大球の直径は 3 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, h, h2) = (3/2, 6/5, 36/5, 24/5)
   plot([b, 0, -b, b], [0, h, 0, 0], color=:blue, lw=0.5)
   circle(0, h2-r1,r1, :green)
   circle(r2, r2, r2)
   circle(-r2, r2, r2, :lightpink)
   circle(0, r2, r2, :lightpink)
   segment(-a, h2, a, h2, :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, h2-r1, "大球:r1,(0,0,h2-r1)", :green, :center, delta=-delta/2)
       point(r2, r2, "小球:r2,(r1,r1,r1)", :red, :center, delta=-delta/2)
       point(0, h2, "h2 ", :blue, :right, :bottom, delta=delta/2)
       point(a, h2, "(a,h2) ", :blue, :right, :bottom, delta=delta/2)
       point(b, 0, "b ", :blue, :right, :bottom, delta=delta/2)
       point(0, h, "h ", :blue, :right, :bottom, delta=delta/2)
   end
end;

function draw2(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, h, h2) = (3/2, 6/5, 36/5, 24/5)
   plot(√2 .* [b, a, -a, -b, b], [0, h2, h2, 0, 0], color=:blue, lw=0.5)
   circle(√2r2, r2, r2)
   circle(-√2r2, r2, r2, :lightpink)
   circle(0, r2, r2, :lightpink)
   circle(0, h2 - r1, r1, :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, h2-r1, "大球:r1,(0,0,h2-r1)", :green, :center, delta=-delta/2)
       point(√2r2, r2, "小球:r2,(√2r1,r1,r1)", :red, :center, delta=-delta/2)
       point(0, h2, "h2 ", :blue, :right, :bottom, delta=delta/2)
       point(√2a, h2, "(√2a,h2) ", :blue, :right, :bottom, delta=delta/2)
       point(√2b, 0, "b ", :blue, :right, :bottom, delta=delta/2)
       #point(0, h, "h ", :blue, :right, :bottom, delta=delta/2)
   end
end;

算額（その874）

秋田県大仙市角間川町 角間川前寺

矢野充志：いろいろな算額，平成30年度 微分積分Ⅰ (2S・2I・2C)，2019/02/25.
https://www.libe.nara-k.ac.jp/~yano/biseki1_2018/meisatsu2.pdf

正方形の内部に斜線，四分円，半円，大円，小円を入れる。小円の直径が与えられたとき，大円の直径を求めよ。

正方形の一辺の長さを 2r1
四分円の半径と中心座標を 2r1, (0, 0)
半円の半径と中心座標を r1, (2r1, r1)
大円の半径と中心座標を r2, (r2, r2)
小円の半径と中心座標を r3, (2r1 - r3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, y3::positive, x::positive
eq1 = (2r1 - r3)^2 + y3^2 - (2r1 + r3)^2
eq2 = r3^2 + (y3 - r1)^2 - (r1 - r3)^2
eq3 = dist(0, 2r1, x, 0, r2, r2) - r2^2
eq4 = dist(0, 2r1, x, 0, 2r1, r1) - r1^2
res = solve([eq1, eq2, eq3, eq4], (r1, r2, y3, x))

   2-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (25*r3/8, 25*r3/16, 5*r3, 75*r3/16)
    (25*r3/8, 75*r3/8, 5*r3, 75*r3/16)

2 組の解が得られるが， 2 番目のものが適解である。

大円の半径は，小円の半径の 25/16 倍である。
小円の直径が 16 のとき，大円の直径は 25 である

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

   r3 = 8;  r1 = 25;  r2 = 12.5;  y3 = 40;  x= 37.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 8
   (r1, r2, y3, x) = r3 .* (25/8, 25/16, 5, 75/16)
   @printf("小円の直径が %g のとき，大円の直径は %g である\n", 2r3, 2r2)
   @printf("r3 = %g;  r1 = %g;  r2 = %g;  y3 = %g;  x= %g\n", r3, r1, r2, y3, x)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:blue, lw=0.5)
   circle(0, 0, 2r1, beginangle=0, endangle=90)
   circle(2r1, r1, r1, :green, beginangle=90, endangle=270)
   circle(r2, r2, r2, :magenta)
   circle(2r1 - r3, y3, r3, :orange)
   segment(0, 2r1, x, 0, :brown)
   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(2r1, 0, "2r1 ", :blue, :right, :bottom, delta=delta/2)
       point(x, 0, " x", :blue, :left, :bottom, delta=delta/2)
       point(0, 2r1, " 2r1", :blue, :left, :bottom, delta=delta/2)
       point(r2, r2, "大円:r2,(r2,r2)", :magenta, :center, delta=-delta/2)
       point(2r1 - r3, y3, "小円:r3,(2r1-r3,y3)", :orange, :center, delta=-delta/2)
       point(2r1, r1, "半円:r1\n(2r1,r1)", :green, :right, :vcenter)
   end
end;

算額（その873）

会田安明：算法天生法指南，文化7年(1810)
京都大学貴重資料デジタルアーカイブ
https://rmda.kulib.kyoto-u.ac.jp/item/rb00028519#?c=0&m=0&s=0&cv=95&r=0&xywh=-487%2C-149%2C4481%2C2977

横長の紙をおみくじ上に結ぶ。紙の幅が一寸のとき正五角形の一辺の長さはいかほどか。

正五角形の一辺の長さ(ab)を x，紙の幅(ac)を y とする。
図において，∠bac = 18° である。
x*cosd(18) = y より，y が与えられたとき x = 2√2y/sqrt(√5 + 5)を求める。

include("julia-source.txt");

using SymPy
@syms x, y
eq = x*cosd(Sym(18)) - y
x = solve(eq, x)[1]
x |> println

   2*sqrt(2)*y/sqrt(sqrt(5) + 5)

紙の幅が 1 寸のとき，正五角形の一辺の長さは 1.0514622242382672 寸である。

x(y => 1).evalf() |> N

   1.0514622242382672

算額（その872）

大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html

大球の中に，小球 3 個と，中球 1 個が入っている。大球と小球の直径がそれぞれ 1 尺 2 寸と 4 寸 8 分のとき，中球の直径を求めよ。

eq1, eq2: y 軸の負の無限大方向から x-z 平面を見ると，小球は大球に内接し，中級と外接している。

eq3: 3 次元空間，x-y-z 軸で考える。z 軸の正の無限大方向から x-y 平面を見ると，小球 3 個は互いに接しており，その中心を結ぶと正三角形になる。

大球の半径と中心座標を r1, (0, 0, 0)
中球の半径と中心座標を r2, (0, 0, r1 - r2)
小球の半径と中心座標を r3, (x3, 0, z3); z3 < 0
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive, z3::negative
(r1, r3) = (120, 48) .// 20
eq1 = x3^2 + z3^2 - (r1 - r3)^2
eq2 = x3^2 + (r1 - r2 - z3)^2 - (r2 + r3)^2
eq3 = x3√Sym(3)/2 - r3
(r2, x3, z3) = solve([eq1, eq2, eq3], (r2, x3, z3))[1]

   (3*sqrt(33)/17 + 39/17, 8*sqrt(3)/5, -2*sqrt(33)/5)

中球の直径は 2(3√33 + 39)/17  = 6.615727992895775 寸である。

2(3√33 + 39)/17

   6.615727992895775

算額の「答」は「中球径七寸二分」となっているが，これは中球と小球の 1 個が大球の直径上に並ぶ解（12 尺 = 7.2 + 4.8）で，残りの 2 個の小球が入る余地がない。「術」の記述はない。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (120, 48) .// 20
   (r2, x3, z3) = (3*sqrt(33)/17 + 39/17, 8*sqrt(3)/5, -2*sqrt(33)/5)
   plot()
   circle(0, 0, r1, :blue)
   circle(0, r1 - r2, r2, :orange)
   circle(x3, z3, r3, :brown)
   if more        
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r1, " r1", :blue, :center, :bottom, delta=delta/2)
       point(0, r1 - r2, "中球:r2,(0,0,r1-r2)", :orange, :center, delta=-delta/2)
       point(x3, z3, "小球:r3,(x3,0,z3)", :brown, :center, delta=-delta/2)
   end
end;

function draw2(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (120, 48) .// 20
   (r2, x3, z3) = (3*sqrt(33)/17 + 39/17, 8*sqrt(3)/5, -2*sqrt(33)/5)
   plot()
   circle(0, 0, r1, :blue)
   rotate(x3, 0, r3, :brown)
   (x, y) = x3*cosd(120), x3*sind(120)
   segment(x, y, x3, 0, :gray90)
   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, r1, " r1", :blue, :center, :bottom, delta=delta/2)
       point(x3, 0, "小球:r3,(x3,0,z3)", :brown, :center, delta=-delta/2)
       point(x, y, "", :brown)
   end
end;

算額（その871）

奈良県大和郡山市 庚申堂 明治13年(1880)
http://www.wasan.jp/nara/kosindo.html

牧下英世：数学教育を通して取り組んだ総合的な学習とその実証的な研究-算額を―算額を用いた課題学習とそのフィールドワークの実践から―，2002筑波大学附属駒場論集第42集
https://tsukuba.repo.nii.ac.jp/record/6466/files/13.pdf

大円と小円が交差してできる隙間に甲円から始まる累円を入れる。
大円，小円，甲円の直径をそれぞれ 7 尺 2 寸，6 尺 1 寸，1 尺 7 寸としたとき，乙円，丙円，丁円，戊円の直径を求めよ。

大円の半径と中心座標を R1, (0, 0)
小円の半径と中心座標を R2, (0, 2r1 - R1 + R2)
甲円の半径と中心座標を r1, (0, r1 - R1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
:

include("julia-source.txt");

using SymPy
@syms R1, R2, r1, r2, x2, y2

(R1, R2, r1) = (72, 61, 17) .// 2
eq1 = x2^2 + y2^2 - (R1 -r2)^2
eq2 = x2^2 + (2r1 - R1 + R2 - y2)^2 - (R2 + r2)^2
eq3 = x2^2 + (y2 - r1 + R1)^2 - (r1 + r2)^2
#(r2, x2, y2) = 
solve([eq1, eq2, eq3], (r2, x2, y2))[2]

   (36465/4681, 204*sqrt(130845)/4681, -109509/4681)

漸化式的に累円のパラメータを数式で求めることが SymPy では困難なので，関数内で solve() を使う，半自動的なやり方を取る。
一つの累円のパラメータが得られたら，それを関数の引数に指定して，次の累円のパラメータを得る。
注意点として（Julia の SymPy 特有であるが），xx/yy の分数は xx//yy にすることと，sqrt(zz) は sqrt(Sym(zz)) にすること。そのようにしなくても良いが，処理時間がかかる。

function nextcircle(R1, R2, r1, r2, x2, y2)
   @syms r3, x3, y3
   eq4 = x3^2 + y3^2 - (R1 -r3)^2
   eq5 = x3^2 + (2r1 - R1 + R2 - y3)^2 - (R2 + r3)^2
   eq6 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2;
   solve([eq4, eq5, eq6], (r3, x3, y3))
end;

# r3
nextcircle(R1, R2, r1, 36465//4681, 204*sqrt(Sym(130845))/4681, -109509//4681)

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (8690825/1404066, 50932*sqrt(130845)/702033, -19854559/1404066)
    (17/2, 0, -55/2)

# r4
nextcircle(R1, R2, r1, 8690825//1404066, 50932*sqrt(Sym(130845))/702033, -19854559//1404066)

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (18641819625/4108026529, 353486916*sqrt(130845)/4108026529, -18850643421/4108026529)
    (36465/4681, 204*sqrt(130845)/4681, -109509/4681)

# r5
nextcircle(R1, R2, r1, 18641819625//4108026529, 353486916*sqrt(Sym(130845))/4108026529, -18850643421//4108026529)

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (4442967010625/1379049831714, 62216188328*sqrt(130845)/689524915857, 4167487120889/1379049831714)
    (8690825/1404066, 50932*sqrt(130845)/702033, -19854559/1404066)

# r6
nextcircle(R1, R2, r1, 4442967010625//1379049831714, 62216188328*sqrt(Sym(130845))/689524915857, 4167487120889//1379049831714)

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (9530164237790625/4195069355668441, 378487752059964*sqrt(130845)/4195069355668441, 35723160673770891/4195069355668441)
    (18641819625/4108026529, 353486916*sqrt(130845)/4108026529, -18850643421/4108026529)

# r7
nextcircle(R1, R2, r1, 9530164237790625//4195069355668441, 378487752059964*sqrt(Sym(130845))/4195069355668441, 35723160673770891//4195069355668441)

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (2271355810006765625/1413029171738162306, 62689304552989212*sqrt(130845)/706514585869081153, 17460791512813260881/1413029171738162306)
    (4442967010625/1379049831714, 62216188328*sqrt(130845)/689524915857, 4167487120889/1379049831714)

# r8
nextcircle(R1, R2, r1, 2271355810006765625//1413029171738162306, 62689304552989212*sqrt(Sym(130845))/706514585869081153, 17460791512813260881//1413029171738162306)

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (4872058212464512265625/4255498594259851496689, 369959068275411941556*sqrt(130845)/4255498594259851496689, 63967589464505474522739/4255498594259851496689)
    (9530164237790625/4195069355668441, 378487752059964*sqrt(130845)/4195069355668441, 35723160673770891/4195069355668441)

累円の直径と中心座標は以下のようになる。算額の答えとはずいぶん異なる。

乙円: 直径 = 15.58000, (15.7641, -23.3944)
丙円: 直径 = 12.37951, (26.2429, -14.1408)
丁円: 直径 =  9.07580, (31.1257, -4.58873)
戊円: 直径 =  6.44352, (32.6386,  3.02200)
己円: 直径 =  4.54351, (32.6356,  8.51551)
庚円: 直径 =  3.21487, (32.0960, 12.35700)
辛円: 直径 =  2.28977, (31.4472, 15.03170)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   names = Char["甲乙丙丁戊己庚辛壬癸"...]
   (R1, R2, r1) = (72, 61, 17) .// 2
   plot(showaxis=false)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   circlef(0, 0, R1, :lightgreen)
   circlef(0, 2r1 - R1 + R2, R2, :cadetblue1)
   circle(0, 2r1 - R1 + R2, R2, :lightslateblue)
   circle(0, 0, R1, :mediumseagreen)
   circlef(0, r1 - R1, r1, :gold)
   number = 1
   for (r, x, y) in [
       (7.790002136295663, 15.764133064156171, -23.394360179448835)
       (6.189755324892134, 26.242896581838256, -14.140759052637128)
       (4.537901470061319, 31.12566720233952, -4.588734587745892)
       (3.221759582902747, 32.63863611252851, 3.021998933649333)
       (2.2717536779012537, 32.63557363626892, 8.515511340831878)
       (1.607437309459633, 32.095998094235, 12.356992949646468)
       (1.1448854005113114, 31.447186145793232, 15.031749640521548)]
       number += 1
       circlef(x, y, r, number)
       circlef(-x, y, r, number)
       point(x, y, names[number], :white, :center, :vcenter, mark=false)
       @printf("%s円: 直径 = %8.5f, (%g, %g)\n", names[number], 2r, x, y)
   end
   if more        
       #hline!([0], color=:gray80, lw=0.5)
       #vline!([0], color=:gray80, lw=0.5)
       point(0, 0, "大円:R1,(0,0)", :red, :center, delta=-1)
       point(0, 2r1 - R1 + R2, "小円:R2,(0,2r1-R1+R2)", :green, :center, delta=-1)
       point(0, r1 - R1, " 甲円:r1\n(0,r1-R1)", :white, :center, delta=-1)
   end
end;