算額（その829）

算額（その829）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に水平な弦，弧 2 個，大円 1 個，中円 1 個，小円 3 個を入れる。弧は全円と同じ半径で，全周の 1/3 である。小円の直径が 13 寸のとき，大円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
弦の両端の座標を (x0, y0)
円弧の半径と中心座標を R, (x0, -y0), (-x0, -y0); x0 = R*cosd(30), y0 = R*sind(30)
大円の半径と中心座標を r1, (x1, y1)
中円の半径と中心座標を r2, (0, R - r2)
小円の半径と中心座標を r3, (x3,  R - 2r2 + r3), (x3,  R - 2r2 - r3)

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, x1::positive, y1::negative,
     r2::positive, r3::positive, x3::positive,
     x0::positive, y0::positive
x0 = R*cosd(Sym(30))
y0 = R*sind(Sym(30))
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = (x1 - x0)^2 + (y1 + y0)^2 - (R - r1)^2
eq3 = x0^2 + (R - 2r2 - r3 + y0)^2 - (R +r3)^2
eq4 = x3^2 + (R - 2r2 + r3)^2 - (R - r3)^2
eq5 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq6 = (x1 + x0)^2 + (y1 + y0)^2 - (R + r1)^2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, r1, x1, y1, r2, x3))

   1-element Vector{NTuple{6, Sym{PyCall.PyObject}}}:
    (16*r3/3, 32*r3/13, 64*sqrt(3)*r3/39, -16*r3/39, 4*r3/3, 4*sqrt(3)*r3/3)

大円の半径 r1 は，小円の半径 r3 の 32/13 倍である。
小円の直径が 13 寸のとき，大円の直径は 32 寸である。

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

R =104/3;  r1 = 16;  x1 =  32√3/3;  y1 = -8/3;  r2 = 26/3;  x3 = 26√3/3

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 13/2
   (R, r1, x1, y1, r2, x3) = [104/3, 16, 32√3/3, -8/3, 26/3, 26√3/3]
   @printf("大円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  x1 = %g;  y1 = %g;  r2 = %g;  x3 = %g\n", R, r1, x1, y1, r2, x3)
   x0 = R*cosd(30)
   y0 = R*sind(30)
   plot()
   circle(0, 0, R, :black)
   circle2(x1, y1, r1, :green)
   circle(0, R - r2, r2)
   circle2(x3, R - 2r2 + r3, r3, :blue)
   circle(0, R - 2r2 - r3, r3, :blue)
   circle(x0, -y0, R, :magenta, beginangle=90, endangle=210)
   circle(-x0, -y0, R, :magenta, beginangle=-30, endangle=90)
   segment(-x0, y0, x0, y0, :magenta)
   if more == true
       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(x0, y0, " (x0,y0)", :magenta, :left, :vcenter)
       point(x0, -y0, " (x0,-y0)", :black, :left, :vcenter)
       point(-x0, -y0, " (-x0,-y0)", :black, :left, :vcenter)
       point(x1, y1, "大円:r1,(x1,y1)", :green, :center, delta=-delta/2)
       point(0, R - r2, "中円:r2\n(0,R-r2)", :red, :center, delta=-delta/2)
       point(x3, R - 2r2 + r3, " 小円:r3,(x3,R-2r2+r3)", :blue, :left, :vcenter)
       point(0, R - 2r2 - r3, " 小円:r3,(0,R-2r2-r3)", :black, :left, :vcenter)
       point(R, 0, " R", :black, :left, :vcenter)
       point(0, R, " R", :black, :center, :bottom, delta=delta/2)
   end
end;

算額（その828）

算額（その828）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に弧を 2 個，甲円 1 個，乙円 3 個を入れる。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (r2 - R), (x2, y2)
全円の円周上にある弧の中心座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, r2::positive, x2::positive, y2::negative, x0::positive, y0::negative
eq1 = x0^2 + (r2 - R - y0)^2 - (R - r2)^2
eq2 = (x0 - x2)^2 + (y2 - y0)^2 - (R - r2)^2
eq3 = x0^2 +(R - r1 - y0)^2 - (r1 + R)^2
eq4 = (x2 + x0)^2 + (y2 - y0)^2 - (R + r2)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = x0^2 + y0^2 - R^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)
   (R, r1, x2, y2, x0, y0) = u
   return [
x0^2 - (R - r2)^2 + (-R + r2 - y0)^2,  # eq1
-(R - r2)^2 + (x0 - x2)^2 + (-y0 + y2)^2,  # eq2
x0^2 - (R + r1)^2 + (R - r1 - y0)^2,  # eq3
-(R + r2)^2 + (x0 + x2)^2 + (-y0 + y2)^2,  # eq4
x2^2 + y2^2 - (R - r2)^2,  # eq5
-R^2 + x0^2 + y0^2,  # eq6
   ]
end;
r2 = 1/2
iniv = BigFloat[2, 1.248, 2.5, -0.7, 0.8, -1.8]
res = nls(H, ini=iniv)

   ([1.3512071919596575, 0.6257139728088581, 0.821951090990108, -0.22124666701222107, 0.821951090990108, -1.0724538589718786], true)

乙円の半径が 1/2 のとき，甲円の半径は 0.6257139728088581 である。
乙円の直径が 1 のとき，甲円の直径は 1.2514279456177162 である。

「術曰」の出だし，「置二分五厘，開立方，名天」は「0.25 の 3 乗根を天と名付ける」でよいだろう。

天 = cbrt(0.25)

   0.6299605249474366

続いて，「四之，」は「これを 4 倍する」
「加一個」は「1を加える」

4天 + 1

   3.5198420997897464

「乗天与一個差」は「1 と天の差を掛ける」（演算数，被演算数が逆になるので注意）

(4天 + 1)*(1 - 天)

   1.3024805228741103

「以除天与一ケ和」は「天と1の和を割る（演算数，被演算数が逆になるので注意）

(天 + 1)/((4天 + 1)*(1 - 天))

   1.2514279456177162

最後に，一般化のために「乗乙円径」で「乙円径を掛ける」

乙円径 = 1
(天 + 1)/((4天 + 1)*(1 - 天))*乙円径

   1.2514279456177162

「答曰」で，「甲円径一寸二分五厘有奇」とあるので，数値解とも一致する。

円弧の描き始めと描き終わりの角度の計算。

@syms R::positive, r1::positive, r2::positive, x2::positive, y2::negative, x0::positive, y0::negative
@syms xa::positive, ya::positive, xb::negative, yb::negative
(R, r1, x2, y2, x0, y0) = [1.3512071919596575, 0.6257139728088581, 0.821951090990108, -0.22124666701222107, 0.821951090990108, -1.0724538589718786]
eq1 = xa^2 + ya^2 - R^2
eq2 = (xa - x0)^2 + (ya - y0)^2 - R^2
res1 = solve([eq1, eq2], (xa, ya))
θ1 = atand(res1[1][2] - y0, res1[1][1] - x0).evalf()
eq1 = xb^2 + yb^2 - R^2
eq2 = (xb - x0)^2 + (yb - y0)^2 - R^2
res2 = solve([eq1, eq2], (xb, yb))
θ2 = atand(res2[1][2] - y0, res2[1][1] - x0).evalf()
(θ1, θ2) |> println

   (67.4673112109364, -172.532688789064)

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (R, r1, x2, y2, x0, y0) = res[1]
   @printf("甲円の直径 = %g\n", 2r1)
   (θ1, θ2) = (67.4673112109364, -172.532688789064)
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r1, r1, :green)
   circle2(x2, y2, r2)
   circle(0, r2 - R, r2)
   circle(x0, y0, R, :blue, beginangle=θ1, endangle=360 + θ2, n=500)
   circle(-x0, y0, R, :blue, beginangle=-180 - θ2, endangle=180 - θ1, n=500)
   if more == true
       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)", :green, :center, delta=-delta)
       point(0, r2 - R, "乙円:r2,(0,r2-R)", :red, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :red, :center, delta=-delta)
       point(x0, y0, " (x0,y0)", :blue, :left, :vcenter)
   end
end;

算額（その827）

算額（その827）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に水平な弦を引き，その上に正三角形を置き，その内部に内接する円を入れる。弦の下に三角形の内部の円と同じ大きさの円（等円）が弦に接し互いにも接し，さらに両端の円は全円とも接している。
等円の直径が 6 寸のとき，全円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, a + r)
弦と y 軸の交点座標を (0, a)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r::positive, a::positive
eq1 = R - a - 3r  # 正三角形の高さは内接する円の半径の 3 倍である
eq2 = (2r)^2 + (a - r)^2 - (R -r)^2  # 両端の等円が全円に内接する
res = solve([eq1, eq2], (R, a))[1]

   (19*r/6, r/6)

全円の半径 R は，等円の半径 r の 19/6 倍である。
等円の直径が 6 寸であれば，全円の直径は 19 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 6/2
   (R, a) = (19*r/6, r/6)
   @printf("全円の直径 = %g\n", 2R)
   plot(√3r .* [1, 0, -1, 1], a .+ [0, 3r, 0, 0], color=:blue, lw=0.5)
   circle(0, 0, R, :black)
   circle(0, a + r, r)
   circle(0, a - r, r, :green)
   circle2(2r, a - r, r, :green)
   x = sqrt(R^2 - a^2)
   segment(-x, a, x, a, :magenta)
   if more == true
       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, a, " a", :magenta, :center, :bottom, delta=delta)
       point(0, a + r, " 等円:r,(0,a+r)", :red, :center, :bottom, delta=delta)
       point(0, a - r, " 等円:r,(0,a-r)", :red, :center, :bottom, delta=delta)
       point(2r, a - r, " 等円:r,(2r,a-r)", :red, :center, :bottom, delta=delta)
       point(0, R, " R", :black, :left, :bottom, delta=delta/2)
   end
end;

算額（その826）

算額（その826）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

正方形の中に四分円と正方形と小円が入っている。小円の直径が 5 寸のとき，内部の正方形の一辺の長さはいかほどか。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r1::positive, r2::positive, a::positive
eq1 = (r1/2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1/2 + a/2)^2 + a^2 - r1^2;
res = solve([eq1, eq2], (r1, a))[1]

   (16*r2, 48*r2/5)

正方形の一辺の長さ a は小円の半径 r2 の 48/5 倍である。
小円の直径が 5 寸のとき，正方形の一辺の長さは 24 寸である。
ちなみに，外の正方形の一辺の長さは 40 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 5/2
   (r1, a) = (16*r2, 48*r2/5)
   @printf("正方形の一辺の長さ = %g\n", a)
   plot(r1 .* [0, 1, 1, 0, 0], r1 .* [0,0,1,1,0], color=:green, lw=0.5)
   circle(r1/2, r1 - r2, r2)
   plot!((r1 - a)/2 .+ [0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5)
   circle(0, 0, r1, :blue, beginangle=0, endangle=90)
   circle(r1, 0, r1, :blue, beginangle=90, endangle=180)
   if more == true
       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 - a)/2, 0, " (r1-a)/2", :magenta, :left, :bottom, delta=delta/2)
       point((r1 + a)/2, 0, " (r1+a)/2", :magenta, :left, :bottom, delta=delta/2)
       point(r1/2, 0, "r1/2", :magenta, :center, :bottom, delta=delta/2)
       point(r1/2, a, "(r1/2,a)", :magenta, :center, :bottom, delta=delta/2)
       point(r1, 0, " r1", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :blue, :left, :bottom, delta=delta/2)
       point(r1/2, r1 - r2, " 小円:r1,(r1/2,r1-r2", :black, :center, :bottom, delta=delta)
   end
end;

算額（その825）

算額（その825）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

直線上に甲円と乙円が載っており互いに接している。乙円の上に正方形が載り，正方形の二辺は甲円に接している。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (0, r2)
正方形の大きさはこの問題では無関係であるが，右端の頂点座標を例えば (x1, 2r2 + x1)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r1::positive, x1::positive, r2::positive
eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = dist(0, 2r2, x1, 2r2 + x1, x1, r1) - r1^2;
res = solve([eq1, eq2], (r1, x1))[1]

   (2*r2, 2*sqrt(2)*r2)

甲円の半径 r1 は乙円の半径 r2 の 2 倍である。
乙円の直径が 1 寸のとき，甲円の直径は 2 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, x1) = (2*r2, 2*sqrt(2)*r2)
   @printf("甲円の直径 = %g\n", 2r1)
   plot()
   circle(0, r2, r2, :green)
   circle2(x1, r1, r1)
   plot!([0, x1, 0, -x1, 0], 2r2 .+ [0, x1, 2x1, x1, 0], color=:blue, lw=0.5)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0, r1], color=:gray80, lw=0.5)
       vline!([0, x1], color=:gray80, lw=0.5)
       point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(0, r2, "乙円:r2\n(0,r2)", :green, :center, delta=-delta/2)
       point(x1, 2r2 + x1, " (x1,2r2+x1)", :blue, :left, :vcenter)
   end
end;

算額（その824）

算額（その824）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に 2 本の斜線と，圭（二等辺三角形），甲円 3 個，乙円 2 個を入れる。乙円の直径が 3 寸のとき，甲円の直径はいかほどか。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative, x0, y0
eq1 = x2^2 + (R - 2r1 +r2)^2 - (R - r2)^2
eq2 = dist(x0, y0, 0, -R, x1, R - 3r1) - r1^2
eq3 = dist(x0, y0, 0, -R, 0, R - r1) - r1^2
eq4 = dist(x0, y0, 0, -R, x2, R - 2r1 + r2) - r2^2
eq5 = dist(sqrt(R^2 - (R - 2r1)^2), R - 2r1, 0, -R, x1, R - 3r1) - r1^2
eq6 = x0^2 + y0^2 - R^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)
   (R, r1, x1, x2, x0, y0) = u
   return [
       x2^2 - (R - r2)^2 + (R - 2*r1 + r2)^2,  # eq1
       -r1^2 + (-x0 + x0*(-x0*(-x0 + x1) + (-R - y0)*(R - 3*r1 - y0))/(x0^2 + (-R - y0)^2) + x1)^2 + (R - 3*r1 - y0 - (-R - y0)*(-x0*(-x0 + x1) + (-R - y0)*(R - 3*r1 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq2
       -r1^2 + (x0*(x0^2 + (-R - y0)*(R - r1 - y0))/(x0^2 + (-R - y0)^2) - x0)^2 + (R - r1 - y0 - (-R - y0)*(x0^2 + (-R - y0)*(R - r1 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq3
       -r2^2 + (-x0 + x0*(-x0*(-x0 + x2) + (-R - y0)*(R - 2*r1 + r2 - y0))/(x0^2 + (-R - y0)^2) + x2)^2 + (R - 2*r1 + r2 - y0 - (-R - y0)*(-x0*(-x0 + x2) + (-R - y0)*(R - 2*r1 + r2 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq4
       -r1^2 + (-r1 - (-2*R + 2*r1)*(-r1*(-2*R + 2*r1) - sqrt(R^2 - (R - 2*r1)^2)*(x1 - sqrt(R^2 - (R - 2*r1)^2)))/(R^2 + (-2*R + 2*r1)^2 - (R - 2*r1)^2))^2 + (x1 + sqrt(R^2 - (R - 2*r1)^2)*(-r1*(-2*R + 2*r1) - sqrt(R^2 - (R - 2*r1)^2)*(x1 - sqrt(R^2 - (R - 2*r1)^2)))/(R^2 + (-2*R + 2*r1)^2 - (R - 2*r1)^2) - sqrt(R^2 - (R - 2*r1)^2))^2,  # eq5
       -R^2 + x0^2 + y0^2,  # eq6
   ]
end;

r2 = 3/2
iniv = BigFloat[7.9, 2.1, 3.5, 3.5, 2.3, 7.7]
res = nls(H, ini=iniv)

   ([8.0, 2.0, 3.4641016151377544, 3.4641016151377544, 2.262270442538942, 7.673469387755102], true)

乙円の直径が 3 寸のとき，甲円の直径は 4 寸である。

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

R = 8;  r1 = 2;  x1 = 3.4641;  x2 = 3.4641;  x0 = 2.26227;  y0 = 7.67347

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 3//2
   (R, r1, x1, x2, x0, y0) = res[1]
   @printf("甲円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  x1 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", R, r1, x1, x2, x0, y0)
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r1, r1)
   circle2(x2, R - 2r1 + r2, r2, :blue)
   circle2(x1, R - 3r1, r1)
   plot!([-x0, 0, x0], [y0, -R, y0], color=:magenta, lw=0.5)
   y00 = R - 2r1
   x00 = sqrt(R^2 - y00^2)
   plot!([-x00, 0, x00, -x00], [y00, -R, y00, y00], color=:magenta, lw=0.5)
   if more == true
       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\n(0,R-r1)", :red, :center, :bottom, delta=delta/2)
       point(x1, R - 3r1, " 甲円:r1\n(x1,R-3r1)", :red, :center, delta=-delta/2)
       point(x2, R - 2r1 + r2, " 乙円:r2\n(x2,R-2r1+r2)", :blue, :center, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :magenta, :left, :bottom, delta=delta/2)
       point(sqrt(R^2 - (R - 2r1)^2), R - 2r1, "sqrt(R^2-(R-2r1)^2),R-2r1 ", :black, :right, delta=-delta/2)
       point(0, R, " R", :green, :left, :bottom, delta=delta/2)
   end
end;

算額（その823）

算額（その823）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に 2 本の斜線と，大円 1 個，小円 3 個を入れる。　大円と小円は，水平な弦に接している。小円の直径が 1 寸のとき，大円の直径はいかほどか。

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

「答」，「術」は水平な弦が全円の中心を通る（つまり直径）と考えているようであるが，実際は水平な弦と y 軸の交点座標は (0, 2r1 - R) である。そして，r1 は R/2 ではない。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, r2::positive,
     x2::positive, x0::positive, y0::negative
sinθ = x0/sqrt(x0^2 + (R - y0)^2)
eq1 = (2R - r1)*sinθ - r1
eq2 = (2R - 2r1 - r2)*sinθ - r2
eq3 = dist(0, R, x0, y0, x2, 2r1 - R - r2) - r2^2
eq4 = x2^2 + (2r1 - R - r2)^2 - (R - r2)^2
eq5 = x0^2 + y0^2 - R^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)
   (R, r1, x2, x0, y0) = u
   return [
       -r1 + x0*(2*R - r1)/sqrt(x0^2 + (R - y0)^2),  # eq1
       -r2 + x0*(2*R - 2*r1 - r2)/sqrt(x0^2 + (R - y0)^2),  # eq2
       -r2^2 + (-x0*(x0*x2 + (-R + y0)*(-2*R + 2*r1 - r2))/(x0^2 + (-R + y0)^2) + x2)^2 + (-2*R + 2*r1 - r2 - (-R + y0)*(x0*x2 + (-R + y0)*(-2*R + 2*r1 - r2))/(x0^2 + (-R + y0)^2))^2,  # eq3
       x2^2 - (R - r2)^2 + (-R + 2*r1 - r2)^2,  # eq4
       -R^2 + x0^2 + y0^2,  # eq5
   ]
end;

r2 = 1/2
iniv = BigFloat[2, 1, 1.4, 1.27, -1.6]
res = nls(H, ini=iniv)

   ([2.00016219313451, 1.0090868249613552, 1.4206243873461804, 1.2703923742568435, -1.5449116525791091], true)

大円の直径は 2.0181736499227103 である。2 ではない。

問には「乃小円者請至多数」とある。小円は3個で，直径が 1 寸と言っているので，この意味がわからなかった。
なお，小円，大円の直径がそれぞれ 1 寸，2 寸のときには，図のイメージはかなり異なるものになる。

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

R = 2.00016;  r1 = 1.00909;  r2 = 0.5;  x2 = 1.42062;  x0 = 1.27039;  y0 = -1.54491

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (R, r1, x2, x0, y0) = res[1]
   @printf("大円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", R, r1, r2, x2, x0, y0)
   plot()
   circle(0, 0, R, :green)
   circle(0, r1 - R, r1)
   circle(0, r2 + 2r1 - R, r2, :blue)
   circle2(x2, 2r1 - R - r2, r2, :blue)
   plot!([-x0, 0, x0], [y0, R, y0], color=:black, lw=0.5)
   y = 2r1 - R
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :green)
   if more == true
       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, 2r1 - R, " 2r1-R", :green, :center, delta=-delta/2)
       point(x2, 2r1 - R - r2, "小円:r2\n(x2,2r1-R-r2)", :blue, :center, delta=-delta/2)
       point(0, r2 + 2r1 - R, "小円:r2\n(0,r2+2r1-R)", :blue, :center, delta=-delta/2)
       point(0, r1 - R, "大円:r1,(0,r1-R)", :red, :center, delta=-delta/2)
       point(0, -R, " -R", :red, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
   end
end;

算額（その822）

算額（その822）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

正方形の中に，四分円，半円，大円，小円を入れる。小円の直径が 16 寸のとき，大円の直径はいかほどか。

四分円の半径と中心座標を r1, (r1, 0)
半円の半径と中心座標を r1/2, (r1, r1/2), (r1/2, r1)
大円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

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

   2-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (33*r3/4, 33*r3/16, -33*sqrt(2)*r3/8 + 33*r3/4, 13*r3/4, 21*r3/4)
    (33*r3/4, 33*r3/16, 33*sqrt(2)*r3/8 + 33*r3/4, 13*r3/4, 21*r3/4)

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

大円の半径 r2 は小円の半径 r3 の 33/16 倍である。
小円の直径が 16 寸のとき，大円の直径は 33 寸である。

2res[1][2](r3 => 16/2).evalf() |> println

   33.0000000000000

なお，大円は 2 つある半円の片方とだけ接する。

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

r1 = 66;  r2 = 16.5;  x2 = 19.331;  x3 = 26;  y3 = 42

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 16//2
   (r1, r2, x2, x3, y3) = r3 .* (33/4, 33/16, 33(2 - √2)/8, 13/4, 21/4)
   @printf("大円の直径 = %g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  x3 = %g;  y3 = %g\n", r1, r2, x2, x3, y3)
   plot(r1 .* [0, 1, 1, 0, 0], r1 .* [0, 0, 1, 1, 0], color=:black, lw=0.5)
   circle(r1, 0, r1, beginangle=90, endangle=180)
   circle(r1, r1/2, r1/2, :magenta, beginangle=90, endangle=270)
   circle(r1/2, r1, r1/2, :magenta, beginangle=180, endangle=360)
   circle(x2, r2, r2, :green)
   circle(x3, y3, r3, :orange)
   if more == true
       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/2, r1, "半円:r1/2,(r1/2,r1)", :magenta, :center, :bottom, delta=delta/2)
       point(r1, r1/2, "半円:r1/2 \n(r1,r1/2) ", :magenta, :right, :vcenter)
       point(x2, r2, "大円:r2,(x2,r2)", :green, :center, delta=-delta/2)
       point(x3, y3, "小円:r3,(x3,y3)", :orange, :center, delta=-delta/2)
   end
end;

算額（その821）

算額（その821）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

外円内に円弧，甲円，乙円を入れる。円弧は外円と同じ半径で周長は外円の 1/3 である。乙円の直径が 3 寸のとき，甲円の直径はいかほどか。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, y1::positive,
     r2::positive, x2::positive, y2::negative
#r2 = 3
eq1 = r1^2 + (y1 + R)^2 - (r1 + R)^2
eq2 = x2^2 + (R + y2)^2 - (r2 + R)^2
eq3 = (x2 - r1)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq4 = r1^2 + y1^2 - (R - r1)^2
eq5 = x2^2 + y2^2 - (R - r2)^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (R, r1, y1, x2, y2));
res[2]

   (4*r2*(109282*sqrt(2107560*sqrt(3) + 3650401) + 63094*sqrt(6322680*sqrt(3) + 10951203) + 319330453 + 184365523*sqrt(3))/(3*(14764*sqrt(3)*sqrt(2107560*sqrt(3) + 3650401) + 25572*sqrt(2107560*sqrt(3) + 3650401) + 47909161*sqrt(3) + 82981101)), r2*(258564*sqrt(2107560*sqrt(3) + 3650401) + 149282*sqrt(6322680*sqrt(3) + 10951203) + 755543964 + 436213511*sqrt(3))/(3*(34932*sqrt(2107560*sqrt(3) + 3650401) + 20168*sqrt(3)*sqrt(2107560*sqrt(3) + 3650401) + 113354292 + 65445131*sqrt(3))), 4*r2*(23094*sqrt(6322680*sqrt(3) + 10951203) + 40000*sqrt(2107560*sqrt(3) + 3650401) + 67482465*sqrt(3) + 116883058)/(3*(14764*sqrt(3)*sqrt(2107560*sqrt(3) + 3650401) + 25572*sqrt(2107560*sqrt(3) + 3650401) + 47909161*sqrt(3) + 82981101)), -260*sqrt(3)*r2*sqrt(62400*sqrt(3)*sqrt(2107560*sqrt(3) + 3650401) + 108080*sqrt(2107560*sqrt(3) + 3650401) + 187572840*sqrt(3) + 324885689) + 1351*r2*sqrt(62400*sqrt(3)*sqrt(2107560*sqrt(3) + 3650401) + 108080*sqrt(2107560*sqrt(3) + 3650401) + 187572840*sqrt(3) + 324885689)/3, -5404*r2*sqrt(6322680*sqrt(3) + 10951203)/9 - 10*sqrt(3)*r2/9 + 2*r2 + 1040*sqrt(3)*r2*sqrt(6322680*sqrt(3) + 10951203)/3)

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

桁数の多い数値を含む複雑な式が得られるが，簡約化すると非常に単純な解になる。
R, r1, y1, x2, y2 の結果は以下のようになる。

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

   28*sqrt(3)*r2/9
   7*r2/3
   -14*sqrt(3)*r2/9 + 14*r2/3
   13*r2/3
   2*r2*(9 - 7*sqrt(3))/9

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

2res[2][2](r2 => 3/2).evalf() |> println

   7.00000000000000

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

R = 8.0829;  r1 = 3.5;  y1 = 2.95855;  x2 = 6.5;  y2 = -1.04145

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 3//2
   (R, r1, y1, x2, y2) = r2 .* (
       28√3/9,
       7/3,
       (42 - 14√3)/9,
       13/3,
       (18 - 14√3)/9)
   @printf("甲円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n", R, r1, y1, x2, y2)
   plot()
   circle(0, 0, R)
   circle(0, -R, R, beginangle=30, endangle=150)
   circle2(r1, y1, r1, :green)
   circle2(x2, y2, r2, :blue)
   if more == true
       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, y1, "甲円:r1,(r1,y1)", :green, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :center, delta=-delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
   end
end;

 

算額（その820）

算額（その820）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

正方形内に四分円，大円，小円を入れる。小円の直径が 4 寸のとき，大円の直径はいかほどか。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r1::positive, r2::positive, r3::positive,
     y3::positive

eq1 = (r1/2)^2 + r2^2 - (r1 - r2)^2
eq2 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = r3^2 + y3^2 - (r1 - r3)^2
res = solve([eq1, eq2, eq3], (r1, r2, y3))[1]

   (6*r3, 9*r3/4, 2*sqrt(6)*r3)

大円の半径 r2 は小円の半径 r3 の 9/4 倍である。
小円の直径が 4 寸なら，大円の直径は 9 寸である。

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

   r1 = 12;  r2 = 4.5;  y3 = 9.79796

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 4//2
   (r1, r2, y3) = (6*r3, 9*r3/4, 2*sqrt(6)*r3)
   @printf("大円の直径 = %g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  y3 = %g\n", r1, r2, y3)
   plot([0, r1, r1, 0, 0], [0, 0, r1, r1, 0], color=:black, lw=0.5)
   circle(0, 0, r1,  beginangle=0, endangle=90)
   circle(r1, 0, r1,  beginangle=90, endangle=180)
   circle(r1 - r3, y3, r3, :green)
   circle(r3, y3, r3, :green)
   circle(r1/2, r2, r2, :blue)
   if more == true
       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/2, r2, "大円:r2,(r1/2,r2)", :blue, :center, delta=-delta)
       point(r1-r3, y3, "小円:r3,(r1-r3,y3)", :green, :center, delta=-delta)
       point(r1, 0, " r1", :red, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :red, :left, :bottom, delta=delta/2)
   end
end;

算額（その819）

算額（その819）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

正三角形内に円弧（弧背），甲円，乙円を入れる。乙円の直径が 2 寸のとき，円弧が最大（最長）になるときの甲円の直径はいかほどか。

甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x2, r2)
円弧がその一部である円の半径と中心座標を r0, (2r1 - r0, 0)
二等辺三角形の一辺の長さと高さを a, √3a
とおき，以下の連立方程式を解く。

なお，円弧が最大というのは，円弧が斜辺に接する状態である。図では ∠BAH が直角である。AB = r0

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a::positive, r0::positive, r1::positive,
     r2::positive, x2::positive

r0 = √Sym(3)a - 2r1
eq1 = 2r1 - (sqrt(Sym(3))a - 3r1)
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = x2^2 + (r2 - 2r1 + r0)^2 - (r0 - r2)^2;

res = solve([eq1, eq2, eq3], (r1, a, x2))

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

甲円の直径は乙円の直径の 3/2 = 1.5 倍である。
乙円の直径が 2 寸なら，甲円の直径は 3 寸である。

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

r1 = 1.5;  a = 4.33013;  x2 = 2.44949;  r0 = 4.5

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2//2
   (r1, a, x2) = (3/2, 5√3/2, √6)
   r0 = √3a - 2r1
   @printf("甲円の直径 = %g\n", 2r1)
   @printf("r1 = %g;  a = %g;  x2 = %g;  r0 = %g\n", r1, a, x2, r0)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, 3r1, r1)
   circle2(x2, r2, r2, :blue)
   circle(0, 2r1 - r0, r0, :green, beginangle=30, endangle=150)
   (x0, y0) = (r0*cosd(30), (2r1 - r0) + r0*sind(30))
   segment(0, 2r1 - r0, x0, y0, :orange)
   if more == true
       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(x0, y0, " A", :green, :left, :vcenter)
       point(-x0, y0, "C ", :green, :right, :vcenter)
       point(a, 0, " a", :black, :left, :vcenter)
       point(0, √3a, " H:√3a", :black, :left, :vcenter)
       point(0, 2r1 - r0, " B:2r1-r0", :orange, :left, :vcenter)
       point(0, r1, "甲円:r1,(0,r1)", :red, :center, delta=-delta)
       point(0, 3r1, "甲円:r1,(0,3r1)", :red, :center, delta=-delta)
       point(x2, r2, "乙円:x2,(0,r2)", :blue, :center, delta=-delta)
       point(0.6x2, 2r1, "弧背", :green, mark=false)
   end
end;

 

算額（その818）

算額（その818）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

長方形と，長方形の長辺を底辺として共有する二等辺三角形の中に，甲円，乙円，丙円を入れる。
丙円の直径が 1 寸のとき，乙円の直径はいかほどか。

甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (0, 2r1 + r2)
丙円の半径と中心座標を r3, (0, r3), (x3, y3)
長方形の長辺と短辺を 2r1, r1
二等辺三角形の高さを h
とおき，以下の連立方程式を解く。

乙円についての方程式 eq5 は，他の方程式とは独立なので，まずは eq1〜eq4 の連立方程式を解いて
h, r1, x3, y3 を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms h::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive, d
r3 = 1//2
eq1 = r1^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = (x3 - r1)^2 + (y3 - r1)^2 - (r1 - r3)^2
eq3 = (x3 - r1)*2r1 - h*(y3 - r1)
eq4 = dist(0, h, 2r1, 0, x3, y3) - r3^2
eq4 = numerator(apart(eq4, d));
(h, r1, x3, y3) = solve([eq1, eq2, eq3, eq4], (h, r1, x3, y3))[2]

   (2*sqrt(2*sqrt(7) + 8)/3 + 2*sqrt(7)*sqrt(2*sqrt(7) + 8)/3, 2, 3*sqrt(2*sqrt(7) + 8)/8 + 2, 3*sqrt(7)/8 + 13/8)

h, r1, x3, y3 が求められた前提で，次の方程式を解く。

eq5 = r2/(h - 2r1 - r2) - 2r1/sqrt(h^2 + 4r1^2)
res2 = solve(eq5, r2)[1]
res2 |> println

   -144/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2)) + 24*sqrt(2*sqrt(7) + 8)/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2)) + 24*sqrt(7)*sqrt(2*sqrt(7) + 8)/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2))

長い式は簡約化できる。

res2 |> factor |> println

   4*(-6 + sqrt(2)*sqrt(sqrt(7) + 4) + sqrt(14)*sqrt(sqrt(7) + 4))/(6 + sqrt(32*sqrt(7) + 128))

乙円の直径は 2.83398951148328 である。

2res2.evalf() |> println

   2.83398951148328

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

   h = 8.861;  r1 = 2;  r2 = 1.41699;  x3 = 3.36716;  y3 = 2.61716

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (h, r1, x3, y3) = (2*sqrt(2*sqrt(7) + 8)/3 + 2*sqrt(7)*sqrt(2*sqrt(7) + 8)/3, 2, 3*sqrt(2*sqrt(7) + 8)/8 + 2, 3*sqrt(7)/8 + 13/8)
   r2 = 4*(-6 + sqrt(2)*sqrt(sqrt(7) + 4) + sqrt(14)*sqrt(sqrt(7) + 4))/(6 + sqrt(32*sqrt(7) + 128))
   @printf("乙円の直径 = %g\n", 2r2)
   @printf("h = %g;  r1 = %g;  r2 = %g;  x3 = %g;  y3 = %g\n", h, r1, r2, x3, y3)
   plot([2r1, 2r1, -2r1, -2r1, 2r1], [0, 2r1, 2r1, 0, 0], color=:black, lw=0.5)
   circle2(r1, r1, r1)
   circle(0, 2r1 + r2, r2, :blue)
   circle(0, r3, r3, :green)
   circle2(x3, y3, r3, :green)
   plot!([-2r1, 0, 2r1], [0, h, 0], color=:orange, lw=0.5)
   if more == true
       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(x3, y3, "丙円:r3,(x3,y3) ", :black, :right, :vcenter)
       point(0, r3, " 丙円:r3,(0,r3)", :black, :left, :vcenter)
       point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(0, 2r1 + r2, "乙円:r2,(0,2r1+r2)", :blue, :center, delta=-delta/2)
       point(0, h, " h", :orange, :left, :vcenter)
       point(2r1, 2r1, "(2r1,2r1)", :black, :right, :bottom, delta=delta/2)
   end
end;

算額（その817）

算額（その817）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

楕円内に等円 2 個が入っている。楕円の短径が 1 のとき，原点を点対称とする円と楕円の 2 接点間の距離を求めよ。

条件が「楕円の短径」のみなので，解けるかどうか不安であったが，「楕円の長径」は無関係であることがわかった。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a::positive, b::positive, x0::positive, y0::positive, r::positive
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/y0
eq3 = (x0 - r)^2 + y0^2 - r^2
res = solve([eq1, eq2, eq3], (x0, y0, r))

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

2 組の解が得られるが，符号の違いのみなのでどちらでもよい。

接点間の距離 diag は以下のようになり，点対称な接点間の距離は，楕円の短径のみで決まる。常に，楕円の短径の √2 倍である。
すなわち，問のように短径が 1 寸ならば，接点間の距離は 1.41421356... である。

@syms d
diag = 2sqrt(res[2][1]^2 + res[2][2]^2)
diag |> println

   2*sqrt(a^2*b^2/(a^2 - b^2) + b^2*((a^2 - 2*b^2)/(a - b))/(a + b))

長い式は，以下のように簡約化される。

apart(diag, d) |> simplify |> println

   2*sqrt(2)*b

b に特定の値を代入して数値解を求める。

apart(diag, d)(b => 1/2) |> simplify |> println

   1.41421356237310

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

   b = 1.234; x0 = 1.35383;  y0 = 1.1012;  r = 1.12477

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (3, 1.234)
   (x0, y0, r) = (a*b/sqrt(a^2 - b^2), b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)
   println("点対称な接点間の距離 = $(2sqrt(x0^2 + y0^2)) = $(2√2b)")
   @printf("b = %g; x0 = %g;  y0 = %g;  r = %g\n", b, x0, y0, r)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle2(r, 0, r, :blue)
   segment(x0, y0, -x0, -y0, :green)
   if more == true
       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, b, " b", :red, :left, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
       point(-x0, -y0, "(-x0,-y0)", :red, :left, :bottom, delta=delta)
       #plot!(xlims=(-0.1,0.1), ylims=(-0.1,0.1))
   end
end;

 

算額（その816）

算額（その816）

文化三年丙寅正月 藤田貞資門人 石州津和野 二介泉尹
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

直角三角形の田がある。直角を挟む二辺の短い方（鈎），長い方（股）が 6 間，8 間のとき，鈎，股に平行な直線で区切ってできる長方形（長辺（長），短辺（平）で，頂点の一つは斜辺上にある）の田の面積が最大になるとき，長，平はいかほどか。

鈎，股を A, B，鈎，股上にある点を a（平）, b（長） とおくと，長方形の面積 S は S = a*b である。
以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms A, B, a, b, S

eq1 = (A - a)/b - A/B
eq2 = a/(B - b) - A/B
eq3 = a*b - S
res = solve([eq1, eq2, eq3], (a, b, S))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (A*(B - b)/B, b, A*b*(B - b)/B)

S = res[1][3]
S |> println

   A*b*(B - b)/B

問のように，A = 6, B = 8 のとき，b の取り方で，面積 S は b の関数で，S(b) は以下のように変化する。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(S(A => 6, B => 8), xlims=(0, 8), xlabel="b", ylabel="S(b)")

面積が最大になるの上図の曲線の接線の傾きが 0 になるときである。数値的に解くには，S(b) を b で微分し，S'(b) = 0 になるときの b の値を求める。

diff(S, b) |> println

   -A*b/B + A*(B - b)/B

solve(diff(S, b), b)[1] |> println

   B/2

b = B/2，a = A/2 が解である。
すなわち，長方形の短辺が 3 間，長辺が 6 間である。

function draw(more)
    pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B) = (6, 8)
   (a, b) = (3, 4)
   plot([0, B, 0, 0], [0, 0, A, 0])
   rect(0, 0, b, a, :red)
   if more == true
       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)
       plot!(showaxis=false)
       point(0, A, " A（鈎）", :blue, :left, :bottom, delta=delta/2)
       point(B, 0, " B（股）", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a（平）", :red, :left, :bottom, delta=delta/2)
       point(b, 0, " b（長）", :red, :left, :bottom, delta=delta/2)
   end
end;

算額（その815）

算額（その815）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

大円内に 5 個の等円を描く。外円と等円の直径をそれぞれ 16 寸，10 寸としたとき，5 個の円の共通部分（黒積）を求めよ。

求めたい面積 S は
S1 = B を中心とする半径 r1 の円の ∠CBE/360 倍
S2 = 三角形 BCE
S3 = 三角形 OCE
とすると，
S = S1 - S2 + S3

まず，2個の等円の交点座標 C:(x0, y0) を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r0, r1, R, x0, y0, OX1, OY1, OX2, OY2
R = r0 - r1
(OX1, OY1) = (R*cosd(Sym(18)), R*sind(Sym(18)))
(OX2, OY2) = (0, R)
eq1 = (x0 - OX1)^2 + (y0 - OY1)^2 - r1^2
eq2 = (x0 - OX2)^2 + (y0 - OY2)^2 - r1^2
res = solve([eq1, eq2], (x0, y0));

x0 = res[2][1] |> factor
y0 = res[2][2] |> factor;

角度 ∠COE を求める

θ = asind(-x0/r1);

S1, S2, S3 を求める。式は長いが，そのままにしておく。

S1 = π*r1^2*θ/360
S2 = (-x0)*(OY2 - y0)/2
S3 = x0*y0/2
S1 |> println
S2 |> println
S3 |> println

   -r1^2*asin((-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*r1*(-5 + sqrt(5))*sqrt(sqrt(5) + 5)))/2
   -(-5*sqrt(2) + sqrt(10))*(r0 - r1 + (sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(4*(-5 + sqrt(5))))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(16*(-5 + sqrt(5))*sqrt(sqrt(5) + 5))
   -(-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))^2/(64*(-5 + sqrt(5))^2*sqrt(sqrt(5) + 5))

面積を求める。

S = 10(S1 - S2 + S3)
S |> println

   -5*r1^2*asin((-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*r1*(-5 + sqrt(5))*sqrt(sqrt(5) + 5))) + 5*(-5*sqrt(2) + sqrt(10))*(r0 - r1 + (sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(4*(-5 + sqrt(5))))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*(-5 + sqrt(5))*sqrt(sqrt(5) + 5)) - 5*(-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))^2/(32*(-5 + sqrt(5))^2*sqrt(sqrt(5) + 5))

r0 = 16//2, r1 = 10//2 を代入して数値を求める。
面積は 13.6341847764931 平方寸（歩）である。

S(r0 => 16//2, r1 => 10//2).evalf() |> println

   13.6341847764931

S(r0 => 1872//2, r1 => 1466//2).evalf() |> println

   915183.000045518

function transform(x, y; deg=0, dx=0, dy=0)
  return [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [x, y]
end;

function rotatefigure(x2, y2)
   for i in 1:5
       plot!(x2, y2, seriestype=:shape, color=:gray70, fillcolor=:gray70, lw=0)
       i == 5 && return
       (x2, y2) = transform(x2, y2, deg = 72)
   end
end;

function fillblack(r1, θ, OX2, OY2)
   x2 = []
   y2 = []
   for i in 270 - θ:0.5:270
       append!(x2, r1*cosd(i) + OX2)
       append!(y2, r1*sind(i) + OY2)
   end
   x2 = vcat(x2, -reverse(x2), [0, x2[1]])
   y2 = vcat(y2,  reverse(y2), [0, y2[1]])
   rotatefigure(x2, y2)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = Vector{Float64}(undef, 6)
   y = Vector{Float64}(undef, 6)
   (r0, r1) = (16, 10) .// 2
   (r0, r1) = (1872, 1466) .// 2
   R = r0 - r1
   for i = 1:6
       θ = 18 + (i - 1)*72
       x[i] = R*cosd(θ)
       y[i] = R*sind(θ)
   end
   (OX1, OY1) = (R*cosd(18), R*sind(18))
   (OX2, OY2) = (0, R)
   x0 = (-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*(-5 + sqrt(5))*sqrt(sqrt(5) + 5))
   y0 = -(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(4*(-5 + sqrt(5)))
   θ = asind(-x0/r1)
   println("θ = $θ")
   S = 10(π*r1^2*θ/360 - (-x0)*(OY2 - y0)/2 + x0*y0/2)
   println("S = $S")
   plot()
   fillblack(r1, θ, OX2, OY2)
   segment(OX2, OY2, x0, y0, :black, lw=1)
   segment(OX2, OY2, 0, OY2 - r1, :black, lw=1.5)
   segment(0, y0, x0, y0, :black, lw=1)
   segment(0, 0, x0, y0, :black, lw=1)
   circle(0, 0, r0)
   circle(0, 0, r0 - r1, :magenta)
   color = [:blue, :green, :gray80, :gray80, :gray80]
   for i = 1:5
       # segment(x[i], y[i], x[i + 1], y[i + 1], :blue)
       circle(x[i], y[i], r1, color[i])
   end
   if more == true
       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(OX1, OY1, " A:r1,(OX1,OY1)", :blue, :left, :bottom, delta=delta/2)
       point(OX2, OY2, " B:r1,(OX2,OY2)", :green, :left, :bottom, delta=delta/2)
       point(x0, y0, "C:(x0,y0) ", :black, :right, :vcenter)
       point(0, OY2 - r1, " D:(0,OY2 - r1) ", :black, :left, delta=-delta/2)
       point(0, y0, " E:(0,y0) ", :black, :left, :bottom, delta=delta/2)
       point(0, 0, " O", :black, :left, :vcenter)
       point(r0 - r1, 0, " r0-r1", :magenta, :left, :bottom, delta=delta/2)
       point(0, r0, " r0", :red, :left, :bottom)
       # plot!(xlims=(-3, 5), ylims=(-3, 3))
   end
end;

   θ = 26.63149441418678
   S = 915183.0000455183

r0，r1 を与えて黒積を返す関数を定義する。

function func(r0, r1)
   p = √5 + 5
   q = √5 - 5
   s = p*(r0 - r1) - √10*sqrt(q*r0*(r0 - 2r1) + r1^2*(√5 + 3))
   t = (√10 - 5√2)*s/(8q√p)
   5(r0*t - r1*t - r1^2*asin(t/r1))
end;
func(8, 5)

   13.634184776493093

この関数を使って，黒積が整数に極めて近い値になるときの r0, r1 の組を探索することもできる。
そのような候補として，r0 = 1872/2, r1 = 1466/2 のとき，黒積が 915183.0000455179 になる。

func(936, 733)

   915183.0000455179