算額（その738）

算額（その738）

千葉県君津市鹿野山 神野寺 明治20年(1887)
山口正義：やまぶき4，第58号
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf

正三角形の中に菱形と円が入っている。正三角形の一辺の長さが 15 寸，円の直径が 4 寸 8 分のとき，菱形の一辺の長さはいかほどか。

正三角形の一辺の長さを 2c，円の半径と中心座標を r, (0, 2b + r)
菱形の対角線の長い方と短い方を 2a, 2b
円の半径と中心座標を r, (0, √3c ‐ 2b ‐ r)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r::positive, a::positive, b::positive, c::positive
eq1 = (sqrt(Sym(3))c - 2b - r)/2 - r
eq2 = b/(c - a) - sqrt(Sym(3))
res = solve([eq1, eq2], (a, b))

   Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}} with 2 entries:
     b => sqrt(3)*c/2 - 3*r/2
     a => c/2 + sqrt(3)*r/2

菱形の一辺の長さ len は sqrt(a^2 + b^2) である。

len = sqrt(res[a]^2 + res[b]^2) |> simplify
len |> println

   sqrt(c^2 - sqrt(3)*c*r + 3*r^2)

@syms 正三角形の一辺の長さ, 円の直径
len2 = len(c => 正三角形の一辺の長さ/2, r => 円の直径/2) |> simplify
len2 |> println

   sqrt(3*円の直径^2 - sqrt(3)*円の直径*正三角形の一辺の長さ + 正三角形の一辺の長さ^2)/2

菱形の一辺の長さは「sqrt(3*円の直径^2 - sqrt(3)*円の直径*正三角形の一辺の長さ + 正三角形の一辺の長さ^2)/2」である。

正三角形の一辺の長さが 15 寸, 円の直径が 4.8 寸のとき，菱形の一辺の長さは 6.50792482007592 寸である。

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

   菱形の一辺の長さ = 6.50792;  対角線の長さ = 11.6569, 5.79038
   c = 7.5;  r = 2.4;  a = 5.82846;  b = 2.89519

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = 15//2
   r = 48//20
   a = (c + √3r)/2
   b = (√3c - 3r)/2
   len = sqrt(a^2 + b^2)
   @printf("菱形の一辺の長さ = %g;  対角線の長さ = %g, %g\n", len, 2a, 2b)
   @printf("c = %g;  r = %g;  a = %g;  b = %g\n", c, r, a ,b)
   plot([c, 0, -c, c], [0, √3c, 0, 0], color=:blue, lw=0.5)
   circle(0, 2b + r, r)
   plot!([0, a, 0, -a, 0], [0, b, 2b, b, 0], color=:green, lw=0.5)
   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, b, " b", :green, :left, :vcenter)
       point(0, √3c, " √3c", :blue, :left, :vcenter)
       point(0, 2b + r, " 2b+r", :red, :left, :vcenter)
       point(0, 2b, " 2b", :red, :center, :bottom, delta=delta/2)
       point(a, b, " (a,b)", :green, :left, :vcenter)
       point(c, 0, "c", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その737）

算額（その737）

千葉県君津市鹿野山 神野寺 明治20年(1887)
山口正義：やまぶき4，第58号
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf

今如圖有鉤股積八百四十零歩只言以股弦差除股弦和得商六個四分之一問鉤股弦各幾何
鈎股弦（直角三角形）において，鈎股積が 840 歩，股弦の差を股弦の和で割ると 6+1/4 になる。鈎股弦はいかほどか。

「鈎股積」は直角三角形の面積のこと，また当然であるが「股弦の差」は「弦 - 股」である。
以下の連立方程式を解く。

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive
eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎*股/2 - 840
eq3 = (股 + 弦)/(弦 - 股) - (6 + 1//4)
res = solve([eq1, eq2, eq3], (鈎, 股, 弦))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (40, 42, 58)

鈎 = 40 寸，股 = 42 寸，弦 = 58 寸である。

算額（その736）

算額（その736）

千葉県君津市鹿野山 神野寺 明治20年(1887)
山口正義：やまぶき4，第58号
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf

半円 2 個と楕円が交差してできる領域に大円と小円をいれる。大円の直径が 6 寸のとき，小円の直径はいかほどか。

明示されていないが楕円の長径は半円の直径と等しいとする。
大円の半径を r1 とすると，半円の直径，楕円の長径は r4 である。
小円の半径と中心座標を r2, (0 2r1 - r2) とおいて，以下の方程式を解く。
見かけは派手な図形であるが，方程式 1 本で，答えもシンプル。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive
a = 2r1
b = 4r1
eq1 = (4r1)^2 + (b - 2r1 - r2)^2 - (r2 + 4r1)^2
res = solve(eq1, r2)

   1-element Vector{Sym{PyCall.PyObject}}:
    r1/3

小円の半径は大円の半径の 1/3 である。
大円の直径が 6 寸のとき，小円の直径は 2 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 6/2
   r2 = r1/3
   @printf("小円の直径 = %g\n", 2r2)
   plot()
   circle(4r1, 0, 4r1, :green, beginangle=90, endangle=270)
   segment(4r1, 4r1, 4r1, -4r1, :green)
   circle(-4r1, 0, 4r1, :green, beginangle=-90, endangle=90)
   segment(-4r1, 4r1, -4r1, -4r1, :green)
   circle42(0, 3r1, r1)
   circle4(r1, 0, r1)
   circle22(0, 2r1 - r2, r2, :magenta)
   ellipse(0, 0, 2r1, 4r1, color=: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(4r1, 0, "4r1 ", :green, :right, :bottom, delta=delta/2)
       point(r1, 0, "r1", :green, :center, :bottom, delta=delta/2)
       point(0, 2r1 - r2, " 小円:r2,(0,2r1-r2)", :black, :left, :vcenter)
       point(0, 3r1, "大円:r1,(0,3r1)", :black, :center, :bottom, delta=delta)
   end
end;

算額（その735）

算額（その735）

福島県田村郡三春町 龍穏院 明治26年(1893)
～落書き帳「○△□」～ 64．手作りの封筒
http://streetwasan.web.fc2.com/math15.5.22.html

一辺の長さが 3 寸の正方形の中に斜線 4 本を描き，区画された領域に宇円4個，宙円 1 個を入れる。それぞれの円の直径を求めよ。

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

include("julia-source.txt");

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

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (a/3, a/5)
    (a, a)

2 組の解が求まるが，最初のものが適解である。

正方形の一辺の長さを 2a としたとき，宇円，宙円の半径は a/3, a/5 である。
正方形の一辺の長さが 3 寸のとき，宇円，宙円の直径は 1 寸，6 分である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3//2
   (r1, r2) = (1/2, 3/10)
   @printf("正方形の一辺の長さ = %g;  宇円の直径 = %g;  宙円の直径 = %g\n", 2a, 2r1, 2r2)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
   circle(0, 0, r2, :blue)
   circle2(a - r1, 0, r1)
   circle22(0, a - r1, r1)
   plot!([-a, 0, a], [a, 2r1 - a, a], color=:green, lw=0.5)
   plot!([-a, 0, a], [-a, a - 2r1, -a], color=:green, lw=0.5)
   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, a, "a", :black, :center, :bottom, delta=delta)
       point(0, r2, "r2", :black, :center, :bottom, delta=delta)
       point(0, a - r1, " a-r1", :red, :center, :bottom, delta=delta)
       point(0, 0, "宙円:r2\n(0,0)", :blue, :center, :bottom, delta=delta)
       point(0, a - 2r1, "a-2r1", :red, :center, :bottom, delta=delta)
       point(a - r1, 0, "宇円:r1,(a-2r1,0)", :red, :center, :bottom, delta=delta)
       point(a, 0, " a", :black, :left, :vcenter)
   end
end;

算額（その734）

算額（その734）

福島県安達郡東和町木幡山（現二本松市） 隠津島神社 明治17年(1884)
～落書き帳「○△□」～ 16．また一寸（○△□のおにぎり弁当）
http://streetwasan.web.fc2.com/math15.5.23.html

直角三角形の中に，正方形，円，正三角形が互いに接して入っている。鈎の長さが 1 寸のとき，「子」の長さはいかほどか。

「鈎」，「股」，「弦」をそのまま変数名として使う。
正三角形の一辺の長さを a，正方形の一辺の長さを b，円の半径を r として以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r::positive, a::positive, b::positive,
     xc::positive, x::positive, y::positive,
     x0::positive, y0::positive, d
弦 = sqrt(鈎^2 + 股^2)
eq1 = 股*y - xc*鈎
eq2 = sqrt(Sym(3))a/(股 - (x + a)) - 鈎/股
eq2 = sqrt(Sym(3))a*股 - 鈎*(股 - (x + a))
eq3 = b*弦 - 股*(鈎 - y)
eq4 = dist(x, 0, x + a, sqrt(Sym(3))a, r, r) - r^2
eq4 = numerator(apart(eq4, d))
eq5 = sqrt((鈎 - y)^2 - b^2) + b + sqrt(3a^2 + (股 - x - a)^2) - 弦
eq6 = y + xc - sqrt(y^2 + xc^2) - 2r
eq7 = dist(x, 0, x + a, sqrt(Sym(3))a, 0, y) - b^2
eq7 = numerator(apart(eq7, d))
eq8 = x0^2 + (y0 - y)^2 - b^2
eq9 = (鈎 - y0)/x0 - 鈎/股;

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, a, b, x, y, xc, x0, y0) = u
   return [
       -xc + y*股,  # eq1
       sqrt(3)*a*股 + a + x - 股,  # eq2
       b*sqrt(股^2 + 1) - 股*(1 - y),  # eq3
       -sqrt(3)*r^2/2 - 3*r*x/2 + sqrt(3)*r*x/2 + 3*x^2/4,  # eq4
       b + sqrt(3*a^2 + (-a - x + 股)^2) + sqrt(-b^2 + (1 - y)^2) - sqrt(股^2 + 1),  # eq5
       -2*r + xc + y - sqrt(xc^2 + y^2),  # eq6
       -b^2 + 3*x^2/4 + sqrt(3)*x*y/2 + y^2/4,  # eq7
       -b^2 + x0^2 + (-y + y0)^2,  # eq8
       -1/股 + (1 - y0)/x0,  # eq9
   ]
end;

鈎 = 1
iniv = BigFloat[1.7, 0.17, 0.37, 0.5, 0.27, 0.46, 0.8, 0.23, 0.87]
res = nls(H, ini=iniv)

   ([1.7320508075688772, 0.16987298107780677, 0.36602540378443865, 0.4641016151377546, 0.2679491924311227, 0.4641016151377546, 0.8038475772933681, 0.2320508075688773, 0.8660254037844386], true)

直角三角形は辺の比が 2:1:√3 のよく見かけるものだ。
鈎が 1 寸のとき，「子」も 1 寸である。

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

   子 = 1.0, 弦 = 1.9999999999999998, b = 0.4641016151377546
   股 = 1.73205;  r = 0.169873;  a = 0.366025;  b = 0.464102;  x = 0.267949;  y = 0.464102;  xc = 0.803848;  x0 = 0.232051;  y0 = 0.866025

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   鈎 = 1
   (股, r, a, b, x, y, xc, x0, y0) = res[1]
   弦 = sqrt(鈎^2 + 股^2)
   println("子 = $(x + 2a), 弦 = $弦, b = $b")
   @printf("股 = %g;  r = %g;  a = %g;  b = %g;  x = %g;  y = %g;  xc = %g;  x0 = %g;  y0 = %g\n", 股, r, a, b, x, y, xc, x0, y0)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r, r, r)
   segment(x, 0, x + a, √3a, :blue)
   segment(x + 2a, 0, x + a, √3a, :blue)
   segment(0, y, xc, 0)
   segment(0, y, x0, y0)
   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, r, "r,(r,r)", :red, :center, delta=-2delta)
       point(x0, y0, " (x0,y0)", :black, :left, :bottom)
       point(x + a, √3a, " (x+a,√3a)", :black, :left, :bottom)
       point(0, 鈎, "鈎 ", :black, :right, :vcenter)
       point(0, y, "y ", :black, :right, :vcenter)
       point(股, 0, "股", :black, :center, delta=-2delta)
       point(x, 0, "x", :black, :center, delta=-2delta)
       point(x + 2a, 0, "x+2a", :black, :center, delta=-2delta)
       point(xc, 0, "xc", :black, :center, delta=-2delta)
       plot!(xlims=(-0.15, 1.9), ylims=(-0.15, 1.1))
       arrow(0, 0, x + 2a, 0, :both)
       # plot!([0, x + 2a], [-0.1, -0.1], c=:red, arrow=:both)#legend=false, arrow=(:both, :open))
       dimension_line(0.0, -0.1, x + 2a, -0.1, "子", color=:red, horizontal=:center, vertical=:bottom, delta=delta/2)
   end
end;

算額（その733）

算額（その733）

福島県安達郡東和町木幡山（現二本松市） 隠津島神社 明治17年(1884)
～落書き帳「○△□」～ 15．一寸
http://streetwasan.web.fc2.com/math15.5.22.html

半径が 1 寸の扇の中に，甲半円と乙円が入っている。乙円の直径が最大になるとき，甲円の直径はいかほどか。

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

include("julia-source.txt");

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

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

2 組の解が求まるが，基本的にはどちらも適解である。
どちらの解も，R - 8r2 の平方根を求める操作が必要なので，乙円の半径 r2 は扇の半径 R の 1/8 よりも大きくはなれない。すなわち r2 の最大値は R/8 であり，このときの r1 が求める答えである。

扇の半径が 1 寸ならば，r2 の最大値は 1/8 = 0.125 であり，このとき r1 は 0.5 寸（直径は 1 寸）である。

res[1][1](R => 1, r2 => 0.125) |> println

   0.500000000000000

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   r2 = R/8
   (r1, x2) = ((-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))
   @printf("扇の半径 = %g;  乙円の直径 = %g;  甲円の直径 = %g\n", R, 2r2, 2r1)
   r3 = R - r1
   y = R - r1
   x = sqrt(R^2 - y^2)
   θ = atand(y/ x)
   plot()
   circle(0, 0, R, beginangle=θ, endangle=180-θ, n=500)
   circle(0, 0, r3, :magenta, beginangle=θ, endangle=180-θ, n=500)
   circle(0, y, r1, :gray80)
   circle(0, y, r1, :blue, beginangle=0, endangle=180)
   circle2(x2, y + r2, r2, :green)
   segment(-x, y, x, y)
   segment(0, 0, x, y)
   segment(0, 0, -x, y)
   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)
       point(0, R - r1, "甲半円:r1,(0,R-r1)", :blue, :center, :bottom, delta=delta)
       point(x2, R - r1 + r2, "乙円:r2,(x2,R-r1+r2) ", :green, :right, :vcenter)
   end
end;

 

算額（その732）

算額（その732）

津島神社 明治16年(1883)
～落書き帳「○△□」～ 190．団扇の中のピタゴラス三角形
http://streetwasan.web.fc2.com/math16.10.11.html

団扇の中に弦を引き，大円，小円 2 個ずつをいれる。
大円と小円の直径の和が 4 寸，弦の長さが 4 寸 8 分のとき，団扇の直径はいかほどか。

水平な弦と y 軸の交点座標を (0, y)
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, y - r1)
小円の半径と中心座標を r2, (r2, y + r2)
円弧（扇の骨が見えるところ）の半径と中心座標を r3, (0, -R)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, y::positive, x::positive, R::positive
eq1 = r1^2 + (y - r1)^2 - (R - r1)^2
eq2 = r2^2 + (y + r2)^2 - (R - r2)^2
eq3 = r1^2 + (y - r1 + R)^2 - (r1 + r3)^2
eq4 = 2r1 + 2r2 - 4
eq5 = 2sqrt(R^2 - y^2) - 48//10

solve([eq1, eq2, eq3, eq4, eq5], (R, r1, r2, r3, y))

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

外円の半径は 5/2 = 2 寸 5 分，直径は 5 寸である。

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

R = 2.5;  r1 = 1.2;  r2 = 0.8;  r3 = 1.13238;  y = 0.7

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2, r3, y) = (5/2, 6/5, 4/5, -6/5 + 2*sqrt(34)/5, 7/10)
   @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  y = %g\n", R, r1, r2, r3, y)
   plot()
   circle(0, 0, R)
   circle2(r1, y - r1, r1, :blue)
   circle2(r2, y + r2, r2, :green)
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y)
   (x2, y2) = (r3*sqrt(4R^2 - r3^2)/(2R), (2R^2 - r3^2)/(2R))
   θ = atand(R - y2, x2)
   circle(0, -R, r3, :black, beginangle=θ, endangle=180-θ)
   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(r2, y + r2, "小円:r2,(r1,y+r2)", :green, :center, delta=-delta)
       point(r1, y - r1, "大円:r1,(r1,y-r1)", :blue, :center, delta=-delta)
       point(0, -R, "弧:r3,(0,-R)", :black, :center, delta=10delta)
       point(0, -R, "-R", :black, :center, :bottom, delta=delta)
       point(0, y, "y", :black, :center, :top, delta=-delta/2)
   end
end;

算額（その731）

算額（その731）

～落書き帳「○△□」～ 330．小平潟天満宮（その３） 明治17年
http://streetwasan.web.fc2.com/math17.11.1.html

等脚台形の中に大円が内接しており，大円に 3 個の楕円が内接している。楕円は隣の楕円と大円に一点で接する。台形の頂点近辺にある 6 個の小円の直径が 1 寸のとき，真ん中にある内円の直径はいかほどか。

「問」には「楕円は隣の楕円と大円に一点で接する」としか書いていないが，「大円と一点で接する『短径が最も長い楕円（曲率楕円）』」である。短径が短い（すなわち内円の直径が大きい）場合には解が定まらない。
『算法助術の公式86』で，直径 d の円に 1 点で内接する曲率楕円の長径 a と短径 b の関係が述べられている d = a^2/b

大円の中心を原点に置く。
大円の半径と中心座標を R, (0, 0)
台形の右上と右下の頂点の座標を (xb, R), (xa, -R)
上方の楕円の長径と短径と中心座標を a, b, (0, r1 + b)
内円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (x1, R - r2), (x2, y2), (x3, r2)
とおき以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R, a, b, r1, x0, y0, xa::positive, xb::positive, r2::positive, x1::positive, x2::positive, y2::negative, x3::positive
@syms d
eq1 = 2b + r1 - R
eq2 = R - a^2/b # b = (R - r1)/2
eq3 = x0^2/a^2 + (y0 - r1 - b)^2/b^2 - 1
eq4 = -b^2*x0/(a^2*(y0 - r1 - b)) - 1/sqrt(Sym(3))
eq5 = x0/y0 - sqrt(Sym(3));
eq6 = dist(xa, -R, xb, R, x1, R - r2) - r2^2
eq7 = dist(xa, -R, xb, R, x2, y2) - r2^2
eq8 = dist(0, 0, xa, - R, x2, y2) - r2^2
eq9 = x3^2 + (r2 - R)^2 - (R + r2)^2;
eq10 = x2^2 + y2^2 - (R + r2)^2
eq11 = xa*xb - R^2;
eq12 = x1^2 + (R - r2)^2 - (R + 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)
   (R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3) = u
   return [
       -R + 2*b + r1,  # eq1
       R - a^2/b,  # eq2
       -1 + (-b - r1 + y0)^2/b^2 + x0^2/a^2,  # eq3
       -sqrt(3)/3 - b^2*x0/(a^2*(-b - r1 + y0)),  # eq4
       x0/y0 - sqrt(3),  # eq5
       -r2^2 + (2*R - 2*R*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) - r2)^2 + (x1 - xa - (-xa + xb)*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2,  # eq6
       -r2^2 + (R - 2*R*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) + y2)^2 + (x2 - xa - (-xa + xb)*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2,  # eq7
       -r2^2 + (x2 - xa*(-R*y2 + x2*xa)/(R^2 + xa^2))^2 + (R*(-R*y2 + x2*xa)/(R^2 + xa^2) + y2)^2,  # eq8
       x3^2 + (-R + r2)^2 - (R + r2)^2,  # eq9
       x2^2 + y2^2 - (R + r2)^2,  # eq10
       -R^2 + xa*xb,  # eq11
       x1^2 + (R - r2)^2 - (R + r2)^2,  # eq12
   ]
end;

r2 = 1/2
iniv = BigFloat[7.1, 1.1, 4.6, 3.1, 3.2, 1.8, 8, 6.1, 5.3, 6.6, -4.5, 3]
res = nls(H, ini=iniv)

   ([3.5, 0.5, 2.29128784747792, 1.5, 1.5155444566227676, 0.875, 3.968626966596886, 3.086709862908689, 2.6457513110645907, 3.307189138830738, -2.25, 2.6457513110645907], true)

小円の直径が 1 寸のとき，内円の直径も 1 寸である。
なお，大円の直径は小円の直径の 7 倍，楕円の短径は小円の直径の 3 倍である。

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

R = 3.5;  r1 = 0.5;  a = 2.29129;  b = 1.5;  x0 = 1.51554;  y0 = 0.875;  xa = 3.96863;  xb = 3.08671;  x1 = 2.64575;  x2 = 3.30719;  y2 = -2.25;  x3 = 2.64575

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3) = res[1]
   @printf("R = %g;  r1 = %g;  a = %g;  b = %g;  x0 = %g;  y0 = %g;  xa = %g;  xb = %g;  x1 = %g;  x2 = %g;  y2 = %g;  x3 = %g\n", R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3)
   plot()
   circle(0, 0, R)
   circle(0, 0, r1, :green)
   l = r1 + b
   ellipse(0, r1 + b, a, b, color=:blue)
   x = l*cos(pi/6)
   y = -l*sin(pi/6)
   ellipse(x, y, a, b, color=:blue, φ= 240)
   ellipse(-x, y, a, b, color=:blue, φ= 120)
   circle(x3, r2 - R, r2)
   circle(-x3, r2 - R, r2)
   segment(-xa, -R, xa, -R)
   segment(-xb, R, xb, R)
   segment(xa, -R, xb, R)
   segment(-xa, -R, -xb, R)
   circle(x1, R - r2, r2)
   circle(-x1, R - r2, r2)
   circle(x2, y2, r2)
   circle(-x2, y2, r2)
   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(x0, y0, "(x0,y0)", :blue, :left, delta=-delta/2)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
       point(xa, -R, "(xa,-R)", :black, :right, :bottom, delta=delta/2)
       point(xb, R, "(xb,R)", :black, :right, :bottom, delta=delta/2)
       point(0, r1, " r1", :blue, :center, :bottom, delta=delta/2)
       point(0, R - b, "R-b", :blue, :center, :bottom, delta=delta/2)
       point(a, R - b, "(a,R-b) ", :blue, :right, :vcenter)
       point(x1, R - r2, "(x1,R-r2)", :red, :right, :bottom, delta=delta)
       point(x2, y2, "小円:r2,(x2,y2)", :red, :right, :vcenter)
       point(x3, r2 - R, " (x3,r2-R)", :red, :left, :vcenter)
       point(0, 0, "内円:r1", :green, :center, :vcenter, mark=false)
   end
end;

 

算額（その730）

算額（その730）

福島県福島市 『黒岩虚空蔵堂算額』第9問 明治26年
街角の数学 ～落書き帳「○△□」～ 347．楕円に内外接する円
http://streetwasan.web.fc2.com/math17.12.5.html

外円内に互いに接している楕円と円 2 個が入っている。楕円の長径と短径，乙円の直径が与えられたとき，甲円の直径を求めよ。

まさに「算法助術の公式85」が該当する問題であるが，図形を描くためのパラメータをすべて得るためには公式85は力不足である。

基本的には，以下の 6 元連立方程式を解けばよいのだが，SymPy の数式処理能力に限界があるので，逐次的に解いていく。

eq1 = (x1 - x2)^2 + y1^2 - (r1 + r2)^2
eq2 = x1^2 + y1^2 - (a - r1)^2
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - r1^2
eq4 = (x0 - x2)^2 + y0^2 - r2^2
eq5 = x0^2/a^2 + y0^2/b^2 - 1
eq6 = -b^2*x0/(a^2*y0) - (x2 - x0)/y0

原点を中心とする楕円と，(x, 0) を中心とする円の接点を (x0, y0) とする。
x0，y0 は同時に求めるとかえってややこしい式になるので，x0 を求める。y0 は円周上の点であることから決定できる。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive,
     a::positive, b::positive, x0::positive, y0::positive
a, b, r2 が既知
y0 = sqrt(r2^2 - (x0 - x2)^2)
eq4 = (x0 - x2)^2 + y0^2 - r2^2
eq6 = -b^2*x0/(a^2*y0) - (x2 - x0)/y0
res1 = solve(eq6, x0)

   1-element Vector{Sym{PyCall.PyObject}}:
    a^2*x2/(a^2 - b^2)

x0, y0 がわかったとして(x2 を含む式であるが），これを方程式に組み込み，x2 を求める。

@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive,
     a::positive, b::positive, x0::positive, y0::positive

x0 = a^2*x2/(a^2 - b^2)
y0 = sqrt(r2^2 - (x0 - x2)^2)
eq5 = x0^2/a^2 + y0^2/b^2 - 1
res2 = solve(eq5, x2)

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

更に x2 を方程式に組み込み，r1, x1, y1 を求める。

@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive,
     a::positive, b::positive, x0::positive, y0::positive

x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
x0 = a^2*x2/(a^2 - b^2)
y0 = sqrt(r2^2 - (x0 - x2)^2)
eq1 = (x1 - x2)^2 + y1^2 - (r1 + r2)^2
eq2 = x1^2 + y1^2 - (a - r1)^2
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - r1^2;
res3 = solve([eq1, eq2, eq3], (r1, x1, y1))

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

以上をまとめると，
x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
x0 = a^2*x2/(a^2 - b^2)
y0 = sqrt(r2^2 - (x0 - x2)^2)
r1 = r2*(a*r2 - b^2)/(2b^2)
x1 = (2a^2*b^2 - 2a^2*r2^2 + a*b^2*r2 - a*r2^3 - b^4 + b^2*r2^2)/(2b*sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2))
y1 = sqrt((a^6*r2^2 - a^4*b^4 - 2a^4*b^2*r2^2 + 2a^2*b^6 + a^2*b^4*r2^2 - b^8)/(a^6 - 3a^4*b^2 + 3a^2*b^4 - b^6))*(a*r2 + b^2)/(2b^2)
となる。

たとえば，楕円の長径と短径が 100, 50，乙円の直径が 40 のとき，甲円の直径は 12 である。

(a, b, r2) = (50, 25, 20)
r1 = r2*(a*r2 - b^2)/(2*b^2)
println("2r1 = $(2r1)")

   2r1 = 12.0

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

   a = 50;  b = 25;  r2 = 20
   r1 = 6;  x1 = 37.2391;  y1 = 23.4361;  x2 = 25.9808;  x0 = 34.641;  y0 = 18.0278

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r2) = (100, 50, 40) ./ 2
   x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
   x0 = a^2*x2/(a^2 - b^2)
   y0 = sqrt(r2^2 - (x0 - x2)^2)
   r1 = r2*(a*r2 - b^2)/(2b^2)
   x1 = (2a^2*b^2 - 2a^2*r2^2 + a*b^2*r2 - a*r2^3 - b^4 + b^2*r2^2)/(2b*sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2))
   y1 = sqrt((a^6*r2^2 - a^4*b^4 - 2a^4*b^2*r2^2 + 2a^2*b^6 + a^2*b^4*r2^2 - b^8)/(a^6 - 3a^4*b^2 + 3a^2*b^4 - b^6))*(a*r2 + b^2)/(2b^2)
   @printf("楕円の長径 = %g;  楕円の短径 = %g;  乙円の直径 = %g;  甲円の直径 = %g\n", 2a, 2b, 2r2, 2r1)
   @printf("a = %g;  b = %g;  r2 = %g\n", a, b, r2)
   @printf("r1 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", r1, x1, y1, x2, x0, y0)
   plot()
   circle(0, 0, a)
   ellipse(0, 0, a, b, color=:blue)
   circle4(x1, y1, r1, :orange)
   circle(x2, 0, r2, :green)
   circle(-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(a, 0, " a", :blue, :left, :vcenter)
       point(0, b, "b", :blue, :center, :bottom, delta=delta/2)
       point(x2, 0, "乙円:r2,(x2,0)", :green, :center, delta=-delta/2)
       point(x0, y0, "(x0,y0) ", :black, :right, delta=-delta/2)
       point(x1, y1, "甲円:r1,(x1,y1) ", :black, :right, :bottom, delta=delta/2)
   end
end;

 

算額（その729）

算額（その729）

富岡村天王社 明治26年
街角の数学 ～落書き帳「○△□」～ 539. 関根熊吉の算額
http://streetwasan.web.fc2.com/math19.1.22.html

直角三角形内に，内接円と直角の頂点から斜辺（弦）への垂線（中鈎）を入れる。
直角三角形の 3 辺の長さの積が 60 寸（立方寸であるが細かいことは言わない），内接円の直径が 2 寸のとき，短弦（弦と中鈎の交点で区切られる弦の短い方）の長さはいかほどか。

直角を挟む短い方の辺と長い方の辺を，「鈎」，「股」，斜辺を「弦」，中鈎，短弦を「中鉤」，「短弦」という変数にする。
内接円の半径と中心座標を r, (股-r, r)
として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, 中鈎::positive, 短弦::positive

r = 2//2
弦 = sqrt(鈎^2 + 股^2)
eq1 = 鈎*股*弦 - 60
eq2 = 鈎 + 股 - 弦 - 2r  # 有名な公式
eq3 = 弦*中鈎 - 股*鈎
eq4 = 短弦/中鈎 - 鈎/股
res = solve([eq1, eq2, eq3, eq4], (鈎, 股, 中鈎, 短弦))

   2-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (3, 4, 12/5, 9/5)
    (4, 3, 12/5, 16/5)

一般に，鈎 < 股 なので，最初のものが適解である。すなわち，短弦は 9/5 寸，すなわち 1 寸 8 分である。

図を描く必要がなければこれで十分である。

図を描くためには中鉤と弦の交点座標があればよい。

交点の座標を (x, y) とする（y 座標は x*鈎/股 であるが，連立させて求める）。

2 条件を追加して連立方程式を解く。

@syms x::positive, y::positive
eq5 = (鈎 - y)/(股 - x) - 鈎/股
eq6 = sqrt(x^2 + y^2) + 短弦 - 弦
res = solve([eq1, eq2, eq3, eq4, eq6, eq5], (鈎, 股, 中鈎, 短弦, x, y))

   2-element Vector{NTuple{6, Sym{PyCall.PyObject}}}:
    (4, 3, 12/5, 16/5, 27/25, 36/25)
    (3, 4, 12/5, 9/5, 64/25, 48/25)

2 番目のものが適解である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 2//2
   (鈎, 股, 中鈎, 短弦, x, y) = (3, 4, 12/5, 9/5, 64/25, 48/25)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   circle(股 - r, r, r, :blue)
   segment(股, 0, x, y, :red)
   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, " 股", :green, :left, :bottom, delta=delta/2)
       point(股, 鈎, "(股,鈎) ", :green, :right, :vcenter)
       point(x, y, "(x,y) ", :green, :right, :vcenter)
       point(股 - r, r, "r,(股-r,r) ", :blue, :center, delta=-delta/2)
   end
end;

算額（その728）

算額（その728）

富岡村天王社 明治26年
街角の数学 ～落書き帳「○△□」～ 544. 関根熊吉の算額（その２）
http://streetwasan.web.fc2.com/math19.1.28.html

大円の中に正三角形，中円，小円が入っている。中円の直径が 5 寸のとき，小円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, y2::positive

#r1 = 5//2
eq1 = r2^2 + y2^2 - (R - r2)^2
eq2 = r2^2 + (y2 - (r1 - R))^2 - (r1 + r2)^2
eq3 = (2R - r1)/2 - r1
res = solve([eq1, eq2, eq3], (r2, y2, R))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (12*r1/25, 9*r1/10, 3*r1/2)

小円の半径は中円の半径の 12/25 倍である。
よって，中円の直径が 5 寸のとき，小円の直径は 2 寸 4 分である。
なお，外円の半径は中円の半径の 3/2 倍である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 5//2
   (r2, y2, R) = r1 .* (12/25, 9/10, 3/2)
   plot([√3R/2, 0, -√3R/2, √3R/2], [-R/2, R, -R/2, -R/2], color=:green, lw=0.5)
   circle(0, 0, R, :blue)
   circle(0, r1 - R, r1, :orange)
   circle(r2, y2, r2)
   circle(-r2, y2, r2)
   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", :blue, :left , :bottom, delta=delta/2)
       point(0, r1 - R, " 中円:r1\n (0,r1-R)", :orange, :left , :vcenter)
       point(r2, y2, " 小円:r2\n (r2,y2)", :red, :left , :vcenter)
   end
end;

算額（その22）Ver.2

算額（その22）Ver.2

算額（その22），算額（その225）を，解析解を求めるように書き換えた。

岩手県遠野市附馬牛町 早池峰神社 弘化3年(1846)
http://www.wasan.jp/iwate/hayatine.html

矢祭町 矢祭神社 明治21年(1888)
街角の数学 ～落書き帳「○△□」～ 499. 矢祭神社算額（その２）
http://streetwasan.web.fc2.com/math18.11.17.html

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（205）
長野県千曲市八幡 八幡八幡神社 文久2年(1862)

一辺の長さが 2 の正三角形に 3 種類の半径を持つ円が図のように接している。それぞれの円の半径を求めよ。

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

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, r2::positive, r3::positive, x2::positive, y3::positive;

eq1 = r3 / (sqrt(Sym(3)) - y3) - sin(PI/6);
eq1 = (y3 + 2r3)^2 - 3 |> expand |> simplify
eq2 = r2　/　(1 - x2) - tan(PI/6);
eq2 = 3r2^2 - ( 1 - x2)^2 |> expand |> simplify
eq3 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2 |> expand |> simplify;
eq4 = x2^2 + (y3 - r2)^2 - (r2 + r3)^2 |> expand |> simplify;
eq5 = (x2 - r1)^2 + (r2 - r1)^2 - (r1 + r2)^2 |> expand |> simplify;

@time res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r3, x2, y3))

   214.393407 seconds (4.97 M allocations: 459.084 MiB, 0.07% gc time, 0.08% compilation time)

   1-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    ((-7268914764*sqrt(-12 + 94*sqrt(3))/13 - 51461127444*sqrt(3)/13 + 101195052708/13 + 5137424024*sqrt(-36 + 282*sqrt(3))/13)/(-4327566283*sqrt(3) - 464793774*sqrt(2)*sqrt(-6 + 47*sqrt(3)) + 661062705*sqrt(6)*sqrt(-6 + 47*sqrt(3)) + 15529289325), (1786462829813*sqrt(-36 + 282*sqrt(3))/3 + 392536725370733/13 + 1329217045078973*sqrt(3)/13 + 15848313004378*sqrt(-12 + 94*sqrt(3)))/(-78432136532557*sqrt(6)*sqrt(-6 + 47*sqrt(3)) + 13232389652808*sqrt(3) + 770162769521882 + 200428587172835*sqrt(2)*sqrt(-6 + 47*sqrt(3))), (-8216629030*sqrt(-36 + 282*sqrt(3))/13 - 131473919301/13 + 105598644225*sqrt(3)/13 + 18423011967*sqrt(-12 + 94*sqrt(3))/13)/(-4327566283*sqrt(3) - 464793774*sqrt(2)*sqrt(-6 + 47*sqrt(3)) + 661062705*sqrt(6)*sqrt(-6 + 47*sqrt(3)) + 15529289325), (-62924*sqrt(3)/13 - 3300*sqrt(-36 + 282*sqrt(3))/13 + 11640*sqrt(-12 + 94*sqrt(3))/13 + 182484/13)/(-3398*sqrt(3) - 89*sqrt(6)*sqrt(-6 + 47*sqrt(3)) + 993*sqrt(2)*sqrt(-6 + 47*sqrt(3)) + 16140), -sqrt(-3/169 + 47*sqrt(3)/338) + 3*sqrt(3)/26 + 27/26)

それぞれを表す数式は大きな整数があったりで長いものである。
たとえば，r1 は以下のようになる。

r1 = res[1][1]
r1 |> println

   (-7268914764*sqrt(-12 + 94*sqrt(3))/13 - 51461127444*sqrt(3)/13 + 101195052708/13 + 5137424024*sqrt(-36 + 282*sqrt(3))/13)/(-4327566283*sqrt(3) - 464793774*sqrt(2)*sqrt(-6 + 47*sqrt(3)) + 661062705*sqrt(6)*sqrt(-6 + 47*sqrt(3)) + 15529289325)

しかし，SymPy である程度簡約化できる。

@syms d
apart(res[1][1], d) |> println
apart(res[1][2], d) |> println
apart(res[1][3], d) |> println
apart(res[1][4], d) |> println
apart(res[1][5], d) |> println

   -sqrt(3)*sqrt(-36 + 282*sqrt(3))/78 - sqrt(-36 + 282*sqrt(3))/52 + 5*sqrt(3)/52 + 45/52
   -sqrt(-36 + 282*sqrt(3))/78 - sqrt(3)*sqrt(-36 + 282*sqrt(3))/156 + 11*sqrt(3)/52 + 21/52
   -27/52 + sqrt(3)*sqrt(-36 + 282*sqrt(3))/156 + 23*sqrt(3)/52
   -21*sqrt(3)/52 + 19/52 + sqrt(-36 + 282*sqrt(3))/52 + sqrt(3)*sqrt(-36 + 282*sqrt(3))/78
   -sqrt(2)*sqrt(-6 + 47*sqrt(3))/26 + 3*sqrt(3)/26 + 27/26

実際に計算するときは，更に手を加えた，以下のような計算式のほうがよいだろう。
t = sqrt(282√3 - 36)  # よく出てくるので事前に計算
r1 = (15√3 + 135 - 3t - 2√3t)/156
r2 = (33√3 + 63 - 2t - √3t)/156
r3 = (69√3 - 81 + √3t)/156
x2 = (57 - 63√3 + 3t + 2√3t)/156
y3 = (3√3 + 27 - sqrt(94√3 - 12))/26

正三角形の一辺の長さが 2 寸のとき，それぞれの値は以下のようになる。
r1 = 0.1505478011587474
r2 = 0.2613764437119185
r3 = 0.4830337288182606
x2 = 0.5472827195892904
y3 = 0.765983349932356

矢祭神社の問題は，2r2 が 3 寸 のとき，2r1 を求めることを要求しているので，以下のようになる。
答は 1.72794226236405 寸である。

(3*res[1][1]/res[1][2]).evalf() |> println

   1.72794226236405

八幡神社の問題は，2r3 が 1 寸のとき，2r2 を求めることを要求しているので，以下のようになる。
答えは 0.541114270325127 寸である。

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

  0.541114270325127

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(282√3 - 36)
   r1 = (15√3 + 135 - 3t - 2√3t)/156
   r2 = (33√3 + 63 - 2t - √3t)/156
   r3 = (69√3 - 81 + √3t)/156
   x2 = (57 - 63√3 + 3t + 2√3t)/156
   y3 = (3√3 + 27 - sqrt(94√3 - 12))/26
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  x2 = $x2;  y3 = $y3")
   plot([-1, 1, 0, -1], [0, 0, sqrt(3), 0], color=:black, linewidth=0.25)
   circle(0, sqrt(3)-2r3, r3)
   circle(r1, r1, r1, :blue)
   circle(-r1, r1, r1, :blue)
   circle(1-sqrt(3)*r2, r2, r2, :green)
   circle(sqrt(3)*r2-1, r2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(r1, r1, "(r1,r1)", :blue, :center)
       point(x2, r2, "(x2,r2)", :green, :center)
       point(0, y3, " y3", :red, :left, :vcenter)
       point(0, y3+r3, " y3+r3", :red, :left, :bottom, delta=delta/2)
       point(0, sqrt(3), " √3", :red, :left, :vcenter)
   end
end;

 

算額（その727）

算額（その727）

安富有恒：和算を楽しむ，和算研究所だより，第12号，Vol. 7, No.2, 2004.
https://i-wasan.sakura.ne.jp/inst/component/tayori_12.pdf

福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html

外円内に楕円が 3 個入っている。中央部と楕円の中に直径の同じ円（等円）が 4 個入っている。

問題を簡単にするために，半径 R，中心座標 (0, 0) の外円に内接する正三角形と，正三角形の辺上に中心を持ち相互に外接する楕円を考える。更に，楕円に外接する円の半径 R2 を求める。楕円の長径と円の直径の比を求めれば，任意の楕円の長径に対する求める円の直径がわかる。

求める外円の直径を R2, (0, 0) とする。
仮の外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, -R/2); r = R/4
楕円の長半径と短半径，中心座標を a, b, (0, -R/2); b = r
楕円同士の接点座標を (x0, y0)，楕円と求める外円の接点座標を (x1, y1)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive, b::positive, r::positive,
x0::positive, y0::negative, x1::positive, y1::negative, R2::positive

R = 4r
b = r
eq1 = -b^2*x0/(a^2*(y0 + R/2)) + 1/√Sym(3)
eq2 = x0^2/a^2 + (y0 + R/2)^2/b^2 - 1
eq3 = y0/x0 + 1/√Sym(3)
eq4 = -b^2*x1/(a^2*(y1 + R/2)) + x1/y1
eq5 = x1^2/a^2 + (y1 + R/2)^2/b^2 - 1
eq6 = x1^2 + y1^2 - R2^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R2, a, x0, y0, x1, y1))

   1-element Vector{NTuple{6, Sym{PyCall.PyObject}}}:
    (3*sqrt(6)*r/2, 3*r, 3*sqrt(3)*r/2, -3*r/2, 3*sqrt(15)*r/4, -9*r/4)

等円の半径が r のとき，外円の半径，楕円の長半径がそれぞれ 3*sqrt(6)*r/2, 3*r なので，外円の直径は楕円の長径の √6/2 倍である。

div(res[1][1], res[1][2]) |> println

   sqrt(6)/2

その他のパラメータは以下のとおりである。なお，楕円の長径は等円の直径の 3 倍である。

R2 = 3.674234614174767
R = 4
r = 1
a = 3
b = 1
x0 = 2.598076211353316
y0 = -1.5
x1 = 2.904737509655563
y1 = -2.25

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1
   R = 4r
   b = r
   (R2, a, x0, y0, x1, y1) = (3√6r/2, 3r, 3√3r/2, -3r/2, 3√15r/4, -9r/4)
   @show(R2, R, r, a, b, x0, y0, x1, y1)
   plot()
   circle(0, 0, R, :gray90)
   plot!([√3R/2, 0, -√3R/2, √3R/2], [-R/2, R, -R/2, -R/2], color=:gray90, lw=0.5)
   circle(0, 0, R2, :orange)
   circle(0, 0, r)
   rotate(0, -2r, r)
   ellipse(0, -R/2, a, r, color=:blue)
   ellipse(R/2*cos(pi/6), R/2*sin(pi/6), a, r, φ=120, color=:blue)
   ellipse(R/2*cos(5pi/6), R/2*sin(5pi/6), a, r, φ=240, color=: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, -2r, "(a,-2r) ", :blue, :right, :vcenter)
       point(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta)
       point(x1, y1, " (x1,y1)", :blue, :left, :vcenter)
       point(0, -R/2, "等円:r,(0,-R/2)", :red, :center, delta=-delta)
       point(r, 0, "r ", :red, :right, :bottom, delta=delta/2)
       point(0, -r, " -r=-R/4", :red, :center, :bottom, delta=delta)
       point(0, -R, " -R", :black, :left, :vcenter)
       point(0, -R2, " -R2", :black, :left, :vcenter)
       # plot!(xlims=(2, 4), ylims=(-2.5, -1))
   end
end;

算額（その726）

算額（その726）

群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm

外円内に弦を引き，菱形と天円と地円を入れる。地円，天円の直径がそれぞれ 7 寸，2 寸のとき，外円の直径はいかほどか。

菱形の対角線の長い方と短い方の長さをそれぞれ 2a, 2b とする。
外円の半径と中心座標を R, (0, 0)
天円の半径と中心座標を r1, (x1, y1)
地円の半径と中心座標を r2, (a, R - 2b + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive, b::positive, d,
     r1::positive, x1::positive, y1::positive, r2::positive
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = dist(0, R, a, R - b, a, R - 2b + r2) - r2^2
eq2 = div(numerator(apart(eq2, d)), b)  # 分子だけを評価する
eq3 = 4(R^2 - (R - 2r1)^2) - (a^2 + b^2)
eq4 = a^2 + (R - b)^2 - R^2
eq5 = b/a - x1/y1;
res = solve([eq1, eq2, eq3, eq4, eq5], (R, a, b, x1, y1));

12 組の解が得られるが，すべてのパラメータの符号を考慮すると 7 番目のものが適解である。

res[7]   

   (8*r1^2/(4*r1 - r2), r2*sqrt(4*r1 + r2)*sqrt(1/(4*r1 - r2)), 4*r1 + r2, (r1 + r2/4)*sqrt(16*r1^2 - r2^2)/(4*r1 - r2), r2*(4*r1 + r2)/(4*(4*r1 - r2)))

外円の半径は「天円の半径を二乗したものを，天円半径の4倍から地円の半径を引いたもので割り，8 倍する」ことで得られる。
天円，地円の半径がそれぞれ 1 寸, 3.5 寸 のとき，外円の半径は 16 寸である。
天円，地円の直径がそれぞれ 2 寸, 7 寸 のとき，外円の直径は 32 寸である。

res[7][1](r1 => 2//2, r2 => 7//2) |> println

   16

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

R = 16;  a = 13.5554;  b = 7.5;  x1 = 7.26184;  y1 = 13.125

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (2, 7) .// 2
   t = 4*r1 - r2
   u = 4*r1 + r2
   (R, a, b, x1, y1) = (8r1^2/t, r2*sqrt(u/t), u, u*sqrt(t*u)/4t, r2*u/(4t))
   @printf("天円の直径 = %g，地円の直径 = %g のとき，外円の直径 = %g\n", r1, r2, 2R)
   @printf("R = %g;  a = %g;  b = %g;  x1 = %g;  y1 = %g\n", R, a, b, x1, y1)
   plot([a, 0, -a, 0, a], [R - b, R, R - b, R - 2b, R - b], color=:green, lw=0.25)
   circle(0, 0, R, :orange)
   circle(x1, y1, r1)
   circle(-x1, y1, r1)
   circle(a, R - 2b + r2, r2, :blue, beginangle=90, endangle=270)
   circle(-a, R - 2b + r2, r2, :blue, beginangle=270, endangle=450)
   y0 = R - 2b
   x0 = sqrt(R^2 - y0^2)
   segment(-x0, y0, x0, y0, :magenta)
   segment(a, y0, a, y0+2r2, :blue)
   segment(-a, y0, -a, y0+2r2, :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)
   end
end;

算額（その725）

算額（その725）

群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm

外円内に弦と斜線を 2 本引き，区画された領域に等円を 4 個入れる。
外円の直径が与えられたとき，等円の直径を求める方法を問う。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, R - 3r), (0, r - R), (0, R - r)
斜線と円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     x0::positive, y0::negative, d
eq0 = x0^2 + y0^2 - R^2
eq1 = x^2 + (R - 3r)^2 - (R - r)^2
eq2 = dist(0, R - 2r, x0, -sqrt(R^2 - x0^2), x, R - 3r) - r^2
eq2 = dist(0, R - 2r, x0, y0, x, R - 3r) - r^2
eq2 = numerator(apart(eq2, d)) |> simplify  # 分子 = 0 で十分
eq3 = dist(0, R - 2r, x0, -sqrt(R^2 - x0^2), 0, r - R) - r^2;
eq3 = dist(0, R - 2r, x0, y0, 0, r - R) - r^2;
eq3 = numerator(apart(eq3, d)) |> simplify  # 分子 = 0 で十分
res = solve([eq0, eq1, eq2, eq3], (r, x, x0, y0))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (-sqrt(17)*R/34 + R/2, sqrt(2)*R*sqrt(-53 + 13*sqrt(17))/4 + 5*sqrt(34)*R*sqrt(-53 + 13*sqrt(17))/68, R*sqrt(-424/17 + 104*sqrt(17)/17), R*(-4 + 13*sqrt(17)/17))

等円の直径は外円の直径の (17 - √17)/34 倍である。

「術」では，「五を開平して...」と言っているので，間違いであろう。

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

   R*(17 - sqrt(17))/34

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

   r = 46.5841;  x = 74.5571;  x0 = 65.3787;  y0 = -104.186

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 123
   (r, x, x0, y0) = (R*(17 - sqrt(17))/34, sqrt(2)*R*sqrt(-53 + 13*sqrt(17))/4 + 5*sqrt(34)*R*sqrt(-53 + 13*sqrt(17))/68, R*sqrt(-424/17 + 104*sqrt(17)/17), R*(-4 + 13*sqrt(17)/17))
   @printf("r = %g;  x = %g;  x0 = %g;  y0 = %g\n", r, x, x0, y0)
   @show(R*(17 - sqrt(17))/34)
   plot()
   circle(0, 0, R, :blue)
   circle(0, R - r, r)
   circle(x, R - 3r, r)
   circle(-x, R - 3r, r)
   circle(0, r - R, r)
   y1 = R - 2r
   x1 = sqrt(R^2 - y1^2)
   segment(-x1, y1, x1, y1, :green)
   segment(0, y1, x0, y0, :green)
   segment(0, y1, -x0, y0, :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(x0, y0, "(x0,y0)", :green, :left, delta=-delta/2)
       point(x, R - 3r, "等円:r,(x,R-3r)", :red, :center, delta=-delta)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(0, R - 2r, "R-2r", :red, :center, :bottom, delta=delta/2)
       point(0, r - R, " r-R", :red, :left, :vcenter)
   end
end;