算額（その1261）

算額（その1261）

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，直角三角形

直角三角形の中に，等円 3 個を容れる。鈎，股が与えられたとき，等円の直径を求めよ。

直角三角形の直角を挟む2辺の短い方を「鈎」，長い方を「股」
斜辺と接する等円の半径と中心座標を r, (2r, (1 + √3)r)
とおき，以下の方程式を解く（斜辺（弦）と等円の距離が等円の半径に等しいとする）。

include("julia-source.txt");

using SymPy
@syms 鈎, 股, r
eq1 = dist2(股, 0, 0, 鈎, 2r, (1 + √Sym(3))r, r)
res = solve(eq1, r)[1]  # 1 of 2
res |> println

   股*鈎*(股 + sqrt(3)*股 + 2*鈎 - sqrt(股^2 + 鈎^2))/(3*股^2 + 2*sqrt(3)*股^2 + 4*股*鈎 + 4*sqrt(3)*股*鈎 + 3*鈎^2)

等円の半径は，鈎，股の関数である。
r = 股*鈎*(股 + sqrt(3)*股 + 2*鈎 - sqrt(股^2 + 鈎^2))/(3*股^2 + 2*sqrt(3)*股^2 + 4*股*鈎 + 4*sqrt(3)*股*鈎 + 3*鈎^2)

鈎 = 3, 股 = 4 のとき，等円の直径は 1.09448091792874 である。

2res(鈎 =>3, 股 => 4).evalf() |> println

   1.09448091792874

function draw(鈎, 股, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 股*鈎*(股 + sqrt(3)*股 + 2*鈎 - sqrt(股^2 + 鈎^2))/(3*股^2 + 2*sqrt(3)*股^2 + 4*股*鈎 + 4*sqrt(3)*股*鈎 + 3*鈎^2)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(r, r, r)
   circle(3r, r, r)
   circle(2r, (1 + √3)r, 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(2r, (1 + √3)r, "等円:r\n(2r,r+√3r)", :red, :center, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta)
       point(股, 0, " 股", :blue, :left, :bottom, delta=delta)
   end
end;

draw(3, 4, true)

 

算額（その1260）

算額（その1260）

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，外円，円弧

外円の一部を折り返し円弧を作り，残りの部分に大円 1 個，小円 2 個を容れる。外円の直径が 1 寸のとき，小円が最大になるときの大円の直径はいかほどか。

本問は図が異なるが，算額（その2053）の「七十八 群馬県甘楽郡下仁田町上小坂 中之嶽神社 安政3年(1856)」と本質的に同じである。
https://blog.goo.ne.jp/r-de-r/e/b7936acdad0ebe7f9934fa840a64eb71

問題の本質をわかりにくくするための細工であろう。外円を折り返した円弧は左側に合同な外円（と円弧）を描いたものと一致する。
「問」の「外円の直径が 1 寸」ということ，「答」の「大円径 4.8584 寸」というのも全く同じである。

include("julia-source.txt");

function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = R*sqrt(-2 + sqrt(5))
   (r2, x2, y2) = (-r1*(-R + r1)*(R + r1)/(R^2 + r1^2), -R*(-R + r1)*(R + r1)/(R^2 + r1^2), 2*R*r1*sqrt(R - r1)*sqrt(R + r1)/(R^2 + r1^2))
   y = sqrt(R^2 - r1^2)
   θ = atand(y, r1)
   plot()
   circle(r1, 0, 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(r1, 0, "外円:R,(r1,0)", :red, :center, delta=-delta/2)
       point(R, 0, "大円:r1,(R,0)", :blue, :center, :bottom, delta=delta/2)
       point(x2, y2, "小円:r2\n(x2,y2)", :green, :center, delta=-delta/2)
   end
   circle(r1, 0, R, :white, beginangle=180-θ, endangle=180+θ)
   circle(-r1, 0, R, beginangle=-θ, endangle=θ)
   circle(R, 0, r1, :blue)
   circle22(x2, y2, r2, :green)
   segment(0, -y, 0, y, :red)
end;

draw(1/2, true)

算額（その1259）

算額（その1259）

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，正方形，長方形

長方形の中に大円 1 個，小円 2 個，正方形 1 個を容れる。小円の直径が 3.432 寸のとき，大円の直径はいかほどか。

大円の半径と中心座標を r1, (-r1, 0)
小円の半径と中心座標を r2, (x2, r1 - r2)
とおき，以下の連立方程式を解く。
長方形の長辺と短辺は 4r1, 2r1 である。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive;
eq1 = dist2(0, 0, r1, r1, x2, r1 - r2, r2)
eq2 = (x2 + r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, x2))[3]  # 3 of 3

   (r2*(sqrt(2) + 3/2), r2/2)

大円の半径 r1 は，小円の半径 r2 の (√2 + 3/2) 倍である。
小円の直径が 3.432余りのとき，大円の直径は 3.432*(√2 + 3/2) = 10.001580946064461 で，ほぼ 10 寸と言いたいのだろう。

整数解にこだわるなら，小円の直径が 1562 寸のとき，大円の直径が 4552.001584426775 になるというのがある。

3.432*(√2 + 3/2)

   10.001580946064461

1562*(√2 + 3/2)

   4552.001584426775

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

draw(3.432/2, true)

算額（その1258）

算額（その1258）

百二十二 群馬県藤岡市東平井 諏訪神社 明治7年(1874)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，斜線2本

斜線 2 本で挟まれた上円，中円，下円とそれらが交差する領域に 2 個の等円を容れる。上円，中円，下円の直径がそれぞれ 3 寸，4 寸，5.4 寸のとき，等円の直径はいかほどか。

斜線（を延長したとき）の交点座標を (a, 0)
等円の半径と中心座標を r4, (x12, 0), (x23, 0); x12 = 2r3 + 2r2 - 3r4; x23 = 2r3 - r4
上円の半径と中心座標を r1, (x1, 0); x1 = 2r3 + 2r2 + r1 - 4r4
中円の半径と中心座標を r2, (x2, 0); x2 = 2r3 + r2 - 2r4
下円の半径と中心座標を r3, (x3, 0); x3 = r3
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive,
     r3::positive, r4::positive,
     x1::positive, x2::positive, x3::positive;
x1 = 2r3 + 2r2 + r1 - 4r4
x2 = 2r3 + r2 - 2r4
x3 = r3
eq1 = r1/(a - x1) - r2/(a - x2)
eq2 = r1/(a - x1) - r3/(a - x3)
res = solve([eq1, eq2], (r4, a));

等円の半径は (r1*r3 - r2^2)/(r1 - 2*r2 + r3) である。
上円，中円，下円の直径がそれぞれ 3 寸，4 寸，5.4 寸のとき，等円の直径は 0.5 寸である。

res[r4] |> println
res[r4](r1 => 3//2, r2 => 4//2, r3 => 54//20) |> println

   (r1*r3 - r2^2)/(r1 - 2*r2 + r3)
   1/4

斜線の交点座標は (18.9, 0) である。

res[a] |> println
res[a](r1 => 3//2, r2 => 4//2, r3 => 54//20) |> println

   2*r3*(-r2 + r3)/(r1 - 2*r2 + r3)
   189/10

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

   r1 = 1.5;  r2 = 2;  r3 = 2.7; , r4 = 0.25;  a = 18.9
   x1 = 9.9;  x2 = 6.9;  x3 = 2.7;  x12 = 8.65;  x23 = 5.15

function draw(r1, r2, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, a) = (1/4, 189/10)
   x1 = 2r3 + 2r2 + r1 - 4r4
   x2 = 2r3 + r2 - 2r4
   x3 = r3
   x12 = 2r3 + 2r2 - 3r4
   x23 = 2r3 - r4
   θ = asind(r1/(a - x1))
   slope = tand(θ)
   @printf("等円の直径は %g である。\n", 2r4)
   @printf("r1 = %g;  r2 = %g;  r3 = %g; , r4 = %g;  a = %g\n", r1, r2, r3, r4, a)
   @printf("x1 = %g;  x2 = %g;  x3 = %g;  x12 = %g;  x23 = %g\n", x1, x2, x3, x12, x23)
   plot()
   circle(x1, 0, r1)
   circle(x2, 0, r2, :blue)
   circle(x3, 0, r3, :green)
   circle(x12, 0, r4, :orange)
   circle(x23, 0, r4, :orange)
   abline(a, 0, slope, 0, 1.05a)
   abline(a, 0, -slope, 0, 1.05a)
   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(x1, 0, " 上円:r1,(x1,0)", :red, :left, :vcenter)
       point(x2, 0, "中円:r2\n(x2,0)", :blue, :center, delta=-3delta)
       point(x3, 0, "下円:r3\n(x3,0)", :green, :center, delta=-3delta)
       point(a, 0, "a", :black, :center, delta=-3delta)
   end
end;

draw(3/2, 4/2, 5.4/2, true)

 

算額（その1257）

算額（その1257）

百二十二 群馬県藤岡市東平井 諏訪神社 明治7年(1874)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，正方形

正方形の中に，中円，東円，西円，南円，北円の 5 個を容れる。東円，西円，南円，北円の直径が 7 寸，6 寸，3 寸，9 寸 のとき，正方形の一辺はいかほどか。

正方形の一辺の長さを a
中円の半径と中心座標を r1, (x1, y1)
東円の半径と中心座標を r2, (a - r2, a - r2)
西円の半径と中心座標を r3, (r3, r3)
南円の半径と中心座標を r4, (a - r4, r4)
北円の半径と中心座標を r5, (r5, a - r5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, x1::positive, y1::positive,
     r2::positive, r3::positive, r4::positive, r5::positive;
eq1 = (a - r2 - x1)^2 + (a - r2 - y1)^2 - (r1 + r2)^2
eq2 = (x1 - r3)^2 + (y1 - r3)^2 - (r1 + r3)^2
eq3 = (a - r4 - x1)^2 + (y1 - r4)^2 - (r1 + r4)^2
eq4 = (x1 - r5)^2 + (a - r5 - y1)^2 - (r1 + r5)^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)
   (a, r1, x1, y1) = u
   return [
-(r1 + r2)^2 + (a - r2 - x1)^2 + (a - r2 - y1)^2,  # eq1
-(r1 + r3)^2 + (-r3 + x1)^2 + (-r3 + y1)^2,  # eq2
-(r1 + r4)^2 + (-r4 + y1)^2 + (a - r4 - x1)^2,  # eq3
-(r1 + r5)^2 + (-r5 + x1)^2 + (a - r5 - y1)^2,  # eq4
   ]
end;

(r2, r3, r4, r5) = [7, 6, 3, 9] ./ 2
iniv = BigFloat[20, 8, 13, 8]

res = nls(H, ini=iniv)

   ([21.0, 7.625, 12.625, 7.5], true)

東円，西円，南円，北円の直径が 7 寸，6 寸，3 寸，9 寸 のとき，正方形の一辺は 21 寸である。

しかし，残念ながら，この算額の問題は欠陥がある。
正方形の一辺の長さは 21 寸になるが，中円が正方形からはみ出す。

北円の直径を 8.9 寸程度にすると，正方形の一辺の長さは19.136 寸程度になり，「群馬の算額」の 140 ページにあるような図になる。

このような齟齬は，得られた解に基づいて正確な図を描いて検証するということがされていなかったことから生まれるのだろう。また，「群馬の算額」も含めて，第三者的立場でのチェックができていないということも根底にある。

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

  r2 = 3.5;  r3 = 3; , r4 = 1.5;  r5 = 4.45;  a = 19.1362;  r1 = 6.37033;  x1 = 11.6481;  y1 = 6.60738

function draw(r2, r3, r4, r5, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1, x1, y1) = res[1]
   @printf("東円，西円，南円，北円の直径が %g，%g，%g，%g のとき，正方形の一辺は %g である。\n",
       2r2, 2r3, 2r4, 2r5, a)
   @printf("r2 = %g;  r3 = %g; , r4 = %g;  r5 = %g;  a = %g;  r1 = %g;  x1 = %g;  y1 = %g\n",
       r2, r3, r4, r5, a, r1, x1, y1)
   plot([0, a, a, 0,  0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(x1, y1, r1)
   circle(a - r2, a - r2, r2, :blue)
   circle(r3, r3, r3, :magenta)
   circle(a - r4, r4, r4, :orange)
   circle(r5, a - r5, r5, :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(a, a, "(a,a)", :green, :right, :bottom, delta=delta/2)
       point(x1, y1, "中円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(a - r2, a - r2, "東円:r2,(a-r2,a-r2)", :blue, :center, delta=-delta/2)
       point(r3, r3, "西円:r3,(r3,r3)", :magenta, :center, delta=-delta/2)
       point(a - r4, r4, "南円:r4,(a-r4,r4)", :black, :right,:bottom, delta=delta/2, deltax=5delta)
       point(r5, a - r5, "北円:r5,(r5,a-r5)", :brown, :center, delta=-delta/2)
   end
end;

draw(7/2, 6/2, 3/2, 9/2, true)

算額（その1256）

算額（その1256）

百十四 群馬県富岡市神成 宇芸神社 明治3年(1870)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，長方形

長方形の中に大円と小円を 2 個ずつ容れる。大円と小円の直径の和が 31.25 寸のとき，小円の直径はいかほどか。

大円の半径と中心座標を r1, (r1, 0)
小円の半径と中心座標を r2, (0, r1 - r2)
とおき，以下の連立方程式を解く。
長方形の長辺と短辺はそれぞれ，4r1, 2r1 である。

include("julia-source.txt");

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

大円の半径 r1 は，径和の 2/5 倍である。
径和が 31.25 寸のとき，大円の直径は 25 である。

res[1] |> println
res[1](径和 => 31.25) |> println

   2*径和/5
   12.5000000000000

小円の半径 r2 は，径和の 1/10 倍である。
径和が 31.25 寸のとき，小円の直径は 6.25 寸である。

res[2] |> println
res[2](径和 => 31.25) |> println

   径和/10
   3.12500000000000

function draw(径和, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (2*径和/5, 径和/10)
   @printf("大円，小円の直径の和が %g のとき，大円，小円の直径は %g, %g である。\n", 径和, 2r1, 2r2)
   plot([2r1, 2r1, -2r1, -2r1, 2r1], [-r1, r1, r1, -r1, -r1], color=:green, lw=0.5)
   circle2(r1, 0, r1)
   circle22(0, r1 - 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,(r1,0)", :red, :center, delta=-delta)
       point(0, r1 - r2, " 小円:r2,(0,r1-r2)", :blue, :left, :vcenter)
       point(2r1, r1, "(2r1,r1)", :green, :right, :bottom, delta=delta)
   end
end;

draw(31.25, true)

算額（その1255）

算額（その1255）

二十五 群馬県高崎市木部町 木部村鎮守社 文化10年(1813)
百十　群馬県高崎市山名町 八幡宮 慶応3年(1867)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，弦，矢

大小 2 個の円が交差している。共通弦が 8 寸，大円，小円の直径の和が 21.6 寸，大円，小円の矢の和が 3.6 寸のとき，大円の直径はいかほどか。

共通弦を「弦」，直径の和を「径和」，矢の和を「矢和」
大円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (0, r1 + r2 - 矢和)
大円と小円の交点座標を (x, y)
とおき，以下の連立方程式を解く。
ただし，一度に解くと SymPy では簡約化できないほど複雑な式になるので，逐次解いていく。

include("julia-source.txt");

using SymPy
@syms x::positive, y::positive,
     r1::positive, r2::positive,
     弦::positive, 径和::positive, 矢和::positive;
x = 弦/2
eq1 = x^2 + y^2 - r1^2
eq2 = x^2 + (y - (r1 + r2 - 矢和))^2 - r2^2
eq3 = r1 + r2 - 径和/2;

まず，eq1 から y，eq3 から r2 を求め，eq2 に代入し eq12 を作る。

ans_y = solve(eq1, y)[1];
ans_y |> println

   sqrt(4*r1^2 - 弦^2)/2

ans_r2 = solve(eq3, r2)[1]
ans_r2 |> println

   -r1 + 径和/2

eq12 = eq2(y => ans_y, r2 => ans_r2);

eq12 を解いて r1 を求める。

弦 = 8, 径和 = 21.6, 矢和 = 3.6 のとき，r1 = 5.8 である。

ans_r1 = solve(eq12, r1)[2]
ans_r1 |> println
ans_r1(弦 => 8, 径和 => 21.6, 矢和 => 3.6) |> println

   (径和*sqrt(矢和)*(径和 - 矢和) + sqrt(-(径和 - 矢和)*(弦^2 - 径和*矢和 + 矢和^2))*(径和 - 2*矢和))/(4*sqrt(矢和)*(径和 - 矢和))
   5.80000000000000

y, r2 は既知の値を代入すれば，求まる。

ans_y |> println
ans_y(弦 => 8, r1 => ans_r1(弦 => 8, 径和 => 21.6, 矢和 => 3.6).evalf()) |> println

   sqrt(4*r1^2 - 弦^2)/2
   4.20000000000000

ans_r2 |> println
ans_r2(径和 => 21.6, r1 => ans_r1(弦 => 8, 径和 => 21.6, 矢和 => 3.6).evalf()) |> println

   -r1 + 径和/2
   5.00000000000000

function draw(弦, 径和, 矢和, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = (径和*sqrt(矢和)*(径和 - 矢和) + sqrt(-(径和 - 矢和)*(弦^2 - 径和*矢和 + 矢和^2))*(径和 - 2*矢和))/(4*sqrt(矢和)*(径和 - 矢和))
   y = sqrt(4*r1^2 - 弦^2)/2
   r2 = -r1 + 径和/2
   x = 弦/2
   plot()
   circle(0, 0, r1)
   circle(0, r1 + 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(0, 0, "大円:r1,(0,0)", :red, :center, delta=-delta/2)
       point(0, r1 + r2 - 矢和, "小円:r1,(0,r1+r2-矢和)", :blue, :center, delta=-delta/2)
       point(x, y, " (x,y)", :green, :left, :vcenter)
       dimension_line(-x, y, x, y, "弦", delta=-2delta, deltax = -3delta)
       dimension_line(0, r1, 0, r1 - 矢和, "矢和", :green, delta=3delta, deltax = 3delta)
       point(0, r1, " r1", :red, :left, :bottom, delta=delta/2)
       point(0, r1 - 矢和, " r1-矢和", :blue, :left, delta=-delta/2)
   end
end;

draw(8, 21.6, 3.6, true)

 

算額（その1254）

算額（その1254）

百八 群馬県邑楽郡板倉町板倉 雷電神社 慶応3年(1867)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，半円2個，長方形

長方形の中に等円 4 個，半円 2 個を容れる。長方形の短辺が 3 寸，長辺が 6 寸のとき，等円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive;
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x2^2 + r2^2 - (r1 - r2)^2
res = solve([eq1, eq2], (r2, x2))[1]
res |> println

   (r1/6, sqrt(6)*r1/3)

等円の半径 r2 は，半円の半径 r1（長方形の短辺） の 1/6 である。
長方形の長辺が 6 寸のとき，等円の直径は 1 寸である。

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = a
   (r2, x2) = (r1/6, sqrt(6)*r1/3)
   @printf("長方形の長辺，短辺が %g, %g のとき，等円の直径は %g である。\n", b, a, 2r2)
   plot([a, a, -a, -a, a], [0, a, a, 0, 0], color=:blue, lw=0.5)
   circle2(x2, r1 - r2, r2)
   circle2(x2, r2, r2)
   circle(0, 0, r1, :green, beginangle=0, endangle=180)
   circle(0, r1, r1, :green, beginangle=180, endangle=360)
   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(x2, r2, "等円:r2\n(x2,r2)", :red, :center, delta=-delta/2)
       point(r1, r1, "(r1,r1)", :blue, :right, :bottom, delta=delta/2)
   end
end;

draw(3, 6, true)

算額（その1253）

算額（その1253）

七十八 群馬県甘楽郡下仁田町上小坂 中之嶽神社 安政3年(1856)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円8個

外円 2 個が交わり，区画された領域に大円 2 個，小円 4 個を容れる。外円の直径が 1 寸のとき，小円の直径が最大になるときの大円の直径はいかほどか。

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

include("julia-source.txt");

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

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

ans_r2 = res[1]
ans_r2 |> println

   -r1*(-R + r1)*(R + r1)/(R^2 + r1^2)

小円の半径 r2 は，大円の半径 r1 と外円の半径 R の関数である。

たとえば，R = 5/2 のときには，下図のように，r1 が 0.2 〜 0.3 の間で r2 が最大になることがわかる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(ans_r2(R => 1/2), xlims=(0, 0.5), xlabel="r1", ylabel="r2")

r1 がどのような値をとると r2 が最大になるかを知るためには，ans_r2 の導関数をとり，それが 0 になるときの r1 を求めればよい。

ans_r1 = solve(diff(ans_r2, r1), r1)[1]
ans_r1 |> println

   R*sqrt(-2 + sqrt(5))

ans_r2 = ans_r2(r1 => ans_r1) |> simplify |> factor
ans_r2 |> println

   R*sqrt(-2 + sqrt(5))*(-1 + sqrt(5))/2

r1 が R*sqrt(-2 + sqrt(5)) のとき，r2 が最大値 R*sqrt(-2 + sqrt(5))*(-1 + sqrt(5))/2 になる。

外円の直径が 1 寸のとき，大円の直径が 0.485868271756646 寸のときに，小円の直径は最大値 0.300283106000778 寸になる。

2ans_r1(R => 1/2).evalf() |> println
2ans_r2(R => 1/2).evalf() |> println

   0.485868271756646
   0.300283106000778

function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = R*sqrt(-2 + sqrt(5))
   (r2, x2, y2) = (-r1*(-R + r1)*(R + r1)/(R^2 + r1^2), -R*(-R + r1)*(R + r1)/(R^2 + r1^2), 2*R*r1*sqrt(R - r1)*sqrt(R + r1)/(R^2 + r1^2))
   @printf("外円の直径が %g のとき，大円の直径が %g のときに小円の直径は最大値 %g になる。\n", 2R, 2r1, 2r2)
   plot()
   circle2(r1, 0, R)
   circle2(R, 0, r1, :blue)
   circle4(x2, y2, 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, 0, "外円:R,(r1,0)", :red, :center, delta=-delta/2)
       point(R, 0, "大円:r1,(R,0)", :blue, :center, :bottom, delta=delta/2)
       point(x2, y2, "小円:r2\n(x2,y2)", :green, :center, delta=-delta/2)
   end
end;

draw(1/2, true)

算額（その1252）

算額（その1252）

七十四 群馬県甘楽郡妙義町菅原 菅原神社 嘉永4年(1851)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円11個

5 個の大円が交差する隙間に 6 個の小円を容れる。大円の直径が 5 寸のとき，小円の直径はいかほどか。



大円の半径と中心座標を r1, (0, r1 + r2)
小円の半径と中心座標を r2, (0, 0), (x2, y2); y2 = x2*tand(54)
とおき，以下の連立方程式を解く。

一度に解くと r2 を表す式がとてつもなく複雑になるので，まず eq2 から x2 を求め，その解を eq1 に代入して r2 を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive, y2::positive;
y2 = x2*tand(Sym(54))
eq1 = x2^2 + (r1 + r2 - y2)^2 - (r1 - r2)^2
eq2 = x2 - (r1 + r2)*cosd(Sym(18))/2;

ans_x2 = solve(eq2, x2)[1];
eq11 = eq1(x2 => ans_x2);

@syms d
ans_r2 = solve(eq11, r2)[1]
ans_r2 = solve(eq11, r2)[1]/r1 |> x -> apart(x, d)*r1 |> simplify
ans_r2 |> println
(ans_r2.evalf()) |> println
ans_r2(r1 => 5/2).evalf() |> println

   r1*(-2*sqrt(10*sqrt(5) + 50) - 4*sqrt(5) + 11 + 4*sqrt(2*sqrt(5) + 10))
   0.259616183682498*r1
   0.649040459206249

小円の半径 r2 は，大円の半径 r1 の (-2*sqrt(10*sqrt(5) + 50) - 4*sqrt(5) + 11 + 4*sqrt(2*sqrt(5) + 10)) = 0.259616183682498 倍である。
大円の直径が 5 寸のとき，小円の直径は 1.29808091841249 寸である。

ans_x2 = ans_x2(r2 => ans_r2) |> simplify
ans_x2 |> println
ans_x2.evalf() |> println
ans_x2(r1 => 5/2).evalf() |> println

   r1*(-sqrt(10*sqrt(5) + 50) - 3*sqrt(5) + 5 + 3*sqrt(2*sqrt(5) + 10))/2
   0.598983089761036*r1
   1.49745772440259

大円の直径が 5 寸のとき，小円の半径の中心座標は (3.743644311006475 寸, 5.152684346344076 寸) である。

1.49745772440259*5/2, 1.49745772440259*5/2*tand(54)

   (3.743644311006475, 5.152684346344076)

function draw(r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = r1*(11 + 4*sqrt(2√5 + 10) - 2sqrt(10√5 + 50) - 4√5)
   x2 = r1*( 5 + 3*sqrt(2√5 + 10) -  sqrt(10√5 + 50) - 3√5)/2
   y2 = x2*tand(54)
   @printf("大円の直径が %g のとき，小円の直径は %g である。\n", 2r1, 2r2)
   plot()
   rotate(0, r1 + r2, r1, angle=72)
   rotate(x2, y2, r2, :blue, angle=72)
   circle(0, 0, 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(x2, y2, " 小円:r2,(x1,y2)", :blue, :left, :vcenter)
       point(0, r1 + r2, "大円:r1,(0,r1+r2)", :red, :center, :bottom, delta=delta/2)
   end
end;

draw(5/2, true)

算額（その1251）

算額（その1251）

七十四 群馬県甘楽郡妙義町菅原 菅原神社 嘉永4年(1851)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：長方形，正五角形3個

長方形の中に正五角形 3 個を容れる。長方形の短辺（直平）が与えられたとき，長方形の長辺（直長）を得る術を述べよ。

問題としては，正五角形の一辺の長さを与えて，直平，直長を求めさせるほうが素直な気がする。

1.　単純に x-y 座標を順に求める方法

正五角形の一辺の長さを a
点 α の座標を (αx, αy)
とおく。

include("julia-source.txt");

using SymPy
@syms a::positive;
cosd2(θ) = 2cosd(θ/2)^2 - 1
sind2(θ) = 2sind(θ/2)*cosd(θ/2)
(s36, s72, s18) = (Sym(36), Sym(72), Sym(18))
(Ax, Ay) = (0, a*sind(s36))
(Bx, By) = (Ax + a*cosd(s36), 0)
(Cx, Cy) = (Bx + a*cosd(s36), By + a*sind(s36))
(Dx, Dy) = (Cx + a, Cy)
(Ex, Ey) = (Dx +a*cosd2(s72), Dy + a*sind2(s72))
(Fx, Fy) = (Ex + a*cosd(s36), Ey + a*sind(s36))
(Gx, Gy) = (Fx - a*cosd2(s72), Fy + a*sind2(s72))
直長 = Fx |> simplify
直平 = Gy |> simplify |> sympy.sqrtdenest |> factor;

直長 |> println
直長(a => 1).evalf() |> println
直平 |> println
直平(a => 1).evalf() |> println

   a*(3 + 2*sqrt(5))/2
   3.73606797749979
   a*sqrt(5 - sqrt(5))*(sqrt(10) + 3*sqrt(2))/4
   3.07768353717525

@syms d
f = apart(直長/直平, d) |> simplify
f |> println
f.evalf() |> println

   sqrt(5 - sqrt(5))*(5*sqrt(2) + 7*sqrt(10))/40
   1.21392207235472

直長は直平の sqrt(5 - √5)*(5√2 + 7√10)/40 = 1.21392207235472 倍である。

3.07768353717525 * sqrt(5 - √5)*(5√2 + 7√10)/40

   3.7360679774997854

2.　正五角形の2点間の距離の累和を求める方法

正五角形の一辺の長さを a，正五角形の外接円の半径を r とする。
CD = a = 2r*sind(36)
PD = r = a/sind(36)
PL = r*cosd(36)
IL = r*(1 + cosd(36))
長辺は AC + CL + IF = EJ + a/2 + EJ = 2(2a*cosd(36)) + a/2
短辺は BN + EM = 2IL = 2(IP + PL) = 2(r + r*cosd(36)) = 2r(1 + cosd(36)) = 2a/2sind(36)*(1 + cosd(36)) = a/sind(36)*(1 + cosd(36))

@syms a::positive, d;
直長 = 2(2a*cosd(Sym(36))) + a/2 |> simplify 
直平 = a/sind(Sym(36))*(1 + cosd(Sym(36))) |> simplify |> factor;

直長 |> println
直長(a => 1).evalf() |> println
直平 |> println
直平(a => 1).evalf() |> println

   a*(3 + 2*sqrt(5))/2
   3.73606797749979
   sqrt(2)*a*(sqrt(5) + 5)/(2*sqrt(5 - sqrt(5)))
   3.07768353717525

直平の式の見かけは前述のものと異なるが，直長/直平の比は同じである。

@syms d
f = apart(直長/直平, d) |> simplify
f |> println
f.evalf() |> println

   sqrt(5 - sqrt(5))*(5*sqrt(2) + 7*sqrt(10))/40
   1.21392207235472

直長は直平の sqrt(5 - √5)*(5√2 + 7√10)/40 = 1.21392207235472 倍である。

3.07768353717525 * sqrt(5 - √5)*(5√2 + 7√10)/40

   3.7360679774997854

3.　正五角形の一辺の長さを a としたときの図(a = 1)

cosd2(θ) = 2cosd(θ/2)^2 - 1
sind2(θ) = 2sind(θ/2)*cosd(θ/2)
function draw(a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   EJ = 2a*cosd(36)
   r = a/2sind(36)
   (Ax, Ay) = (0, a*sind(36))
   (Bx, By) = (Ax + a*cosd(36), 0)
   (Cx, Cy) = (Bx + a*cosd(36), By + a*sind(36))
   (Dx, Dy) = (Cx + a, Cy)
   (Ex, Ey) = (Dx +a*cosd2(72), Dy + a*sind2(72))
   (Fx, Fy) = (Ex + a*cosd(36), Ey + a*sind(36))
   (Gx, Gy) = (Fx - a*cosd2(72), Fy + a*sind2(72))
   (Hx, Hy) = (Gx - a, Gy)
   (Ix, Iy) = (Fx - EJ, Fy)
   (Jx, Jy) = (Ex - EJ, Ey)
   (Kx, Ky) = (Jx - a, Jy)
   直長 = a*(3 + 2√5)/2
   直平 = √2*a*(√5 + 5)/(2*sqrt(5 - √5))
   plot([Ax, Bx, Cx, Dx, Ex, Fx, Gx, Hx, Ix, Jx, Kx, Ax],
        [Ay, By, Cy, Dy, Ey, Fy, Gy, Hy, Iy, Jy, Ky, Ay],
        color=:blue, lw=0.5)
   plot!([0, 直長, 直長, 0, 0], [0, 0, 直平, 直平, 0], color=:green, lw=0.5)
   segment(Ex, Ey, Ix, Iy, :blue)
   segment(Cx, Cy, Jx, Jy, :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)
       circle(EJ + r*sind(36), Iy - r, r, :red)
       point(Ax, Ay, " A", :green, :left, :vcenter)
       point(Bx, By, " B")
       point(Cx, Cy, " C")
       point(Dx, Dy, " D")
       point(Ex, Ey, " E")
       point(Fx, Fy, " F")
       point(Gx, Gy, " G")
       point(Hx, Hy, " H")
       point(Ix, Iy, "I ", :green, :right, :bottom)
       point(Jx, Jy, " J")
       point(Kx, Ky, " K")
       point(EJ + r*sind(36), Iy - r, " P", :red)
       point(EJ + r*sind(36), Cy, " L", :green)
   else
       plot!(showaxis=false)
   end
end;

draw(1, true)

draw(1, false)

算額（その1250）

算額（その1250）

七十四 群馬県甘楽郡妙義町菅原 菅原神社 嘉永4年(1851)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，正方形3個

2 個の円が交わっており，区画された領域に正方形 3 個を容れる。円の直径が与えられたとき，正方形の一辺の長さを求めよ。

算額（その429）から小円を除いた問題である。小円の有無は円と正方形の大きさには無関係である。

正方形の一辺の長さを 2a
円の半径と中心座標を r1, (x1, 0), (-x1, 0)
として，以下の連立方程式を解く。
コメントアウトしたのは算額（その429）でのもの。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     a::positive;

# eq1 = x1^2 + (a + r2)^2 - (r1 - r2)^2
eq2 = a^2 + (x1 + a)^2 - r1^2
eq3 = (r1 - 2x1 + 2a)^2 + a^2 - r1^2
# res = solve([eq1, eq2, eq3], (r2, x1, a))
res = solve([eq2, eq3], (x1, a))[1]

   ((-919*r1^3*(-4/25 + 6*sqrt(6)/25)^2 - 1700*r1^3*(-4/25 + 6*sqrt(6)/25)^3 + 180*r1^3 + 484*r1^3*(-4/25 + 6*sqrt(6)/25))/(180*r1^2), r1*(-4/25 + 6*sqrt(6)/25))

x1, a は簡約化できる。

res[1] |> simplify |> println

   r1*(2*sqrt(6) + 7)/25

res[2] |> simplify |> println

   2*r1*(-2 + 3*sqrt(6))/25

正方形の一辺の長さ 2a は 円の直径 2r1 の 2(3√6 - 2)/25 倍である。
円の直径が 6.17 のとき，正方形の一辺の長さは 2.64000441111333 である。

6.17*2(3√6 - 2)/25

   2.64000441111333

function draw(r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x1 = r1*(2√6 + 7)/25
   a  = 2r1*(3√6 - 2)/25
   @printf("円の直径が %g のとき，正方形の一辺の長さは %g である。\n", 2r1, 2a)
   @printf("x1 = %.15g;  a = %.15g\n", x1, a)
   plot()
   circle(x1, 0, r1, :red)
   circle(-x1, 0, r1, :red)
   rect(-a, -a, a, a, :blue)
   # # circle(0, a + r2, r2, :green)
   rect(r1 - x1, -a, r1 - x1 + 2a, a, :blue)
   rect(-r1 + x1, -a, -r1 + x1 - 2a, a, :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 - x1, 0, " r1-x1", :red, :left, delta=-delta)
       point(x1 - r1, 0, "x1-r1", :red, :right, delta=-delta)
       point(x1, 0, "x1 ", :red, :right, :bottom, delta=delta/2)
       point(-a, a, " (-a,a)", :blue, :left, delta=-delta/2)
       point(r1 - x1 + 2a, a, "(r1-x1+2a,a) ", :blue, :right, delta=-delta/2)
       # # point(0, a + r2, " a+r2", :green)
   end
end;

draw(6.17/2, true)

 

算額（その1249）

算額（その1249）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個

天円 3 個，地円 2 個，人円 1 個が互いに接し合っている。天円，人円の直径を与えたときに地円の直径を求めるすべを述べよ。

天円の半径と中心座標を r1, (x1, y1)
地円の半径と中心座標を r2, (r2, 0)
人円の半径と中心座標を r3, (0, y3)
とおき，以下の連立方程式を解く。
SymPy では連立方程式を一度に解くことはできないので，順に解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, x1::positive, y1::positive,
     r2::positive, r3::positive, y3::positive
eq1 = x1^2 + (y3 + r3 + r1 - y1)^2 - 4r1^2 |> expand
eq2 = (x1 - r2)^2 + y1^2 - (r1 + r2)^2 |> expand
eq3 = x1^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq4 = r2^2 + y3^2 - (r2 + r3)^2 |> expand;

まず，eq1, eq3 から，x1, x2 を求める。

res = solve([eq1, eq3], (x1, y1))[1]  # 1 of 1

   (2*r1*sqrt(r3)*sqrt(2*r1 + r3)/(r1 + r3), (-r1^2 + 2*r1*r3 + r1*y3 + r3^2 + r3*y3)/(r1 + r3))

eq4 に x1, y1 を代入し y3 を求める。

eq14 = eq4(x1 => res[1], y1 => res[2])
ans_y3 = solve(eq14, y3)[1]  # 1 of 1
ans_y3 |> println

   sqrt(r3)*sqrt(2*r2 + r3)

eq2 に x1, y1, y3 を代入し r2 を求める。

eq12 = eq2(x1 => res[1], y1 => res[2], y3 => ans_y3) |> simplify |> numerator;
ans_r2 = solve(eq12, r2)[2]  # 2 of 2
ans_r2 |> println

   2*r1*(r1^5*r3 + r1^4*r3^(3/2)*sqrt(2*r1 + r3) + r1^4*r3^2 + 4*r1^3*r3^(5/2)*sqrt(2*r1 + r3) + 6*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 2*r1^2*r3^4 + 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 3*r1*r3^5 + r3^(11/2)*sqrt(2*r1 + r3) + r3^6 + sqrt(r1^10*r3^2 + 2*r1^9*r3^(5/2)*sqrt(2*r1 + r3) + 2*r1^9*r3^3 + 2*r1^8*r3^(7/2)*sqrt(2*r1 + r3) - 7*r1^8*r3^4 - 16*r1^7*r3^(9/2)*sqrt(2*r1 + r3) - 24*r1^7*r3^5 - 32*r1^6*r3^(11/2)*sqrt(2*r1 + r3) - 10*r1^6*r3^6 + 12*r1^5*r3^(13/2)*sqrt(2*r1 + r3) + 52*r1^5*r3^7 + 92*r1^4*r3^(15/2)*sqrt(2*r1 + r3) + 102*r1^4*r3^8 + 112*r1^3*r3^(17/2)*sqrt(2*r1 + r3) + 88*r1^3*r3^9 + 64*r1^2*r3^(19/2)*sqrt(2*r1 + r3) + 41*r1^2*r3^10 + 18*r1*r3^(21/2)*sqrt(2*r1 + r3) + 10*r1*r3^11 + 2*r3^(23/2)*sqrt(2*r1 + r3) + r3^12))/(r1^6 + 4*r1^5*sqrt(r3)*sqrt(2*r1 + r3) + 10*r1^5*r3 + 8*r1^4*r3^(3/2)*sqrt(2*r1 + r3) + 19*r1^4*r3^2 + 12*r1^3*r3^3 - 8*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 3*r1^2*r3^4 - 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 2*r1*r3^5 + r3^6)

式が長く，SymPy では簡約化できないが，天円，人円の半径を与えれば，地円の半径が求まる。
天円，人円の直径が 15 寸，7 寸のとき，地円の直径は 5.72928861346286 寸である。

2ans_r2(r1 => 15/2, r3 => 7/2) |> println

   5.72928861346286

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

  r1 = 7.5;  r3 = 3.5;  r2 = 2.86464;  y3 = 5.68353;  x1 = 10.9728;  y1 = 6.45626

function draw(r1, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2*r1*(r1^5*r3 + r1^4*r3^(3/2)*sqrt(2*r1 + r3) + r1^4*r3^2 + 4*r1^3*r3^(5/2)*sqrt(2*r1 + r3) + 6*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 2*r1^2*r3^4 + 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 3*r1*r3^5 + r3^(11/2)*sqrt(2*r1 + r3) + r3^6 + sqrt(r1^10*r3^2 + 2*r1^9*r3^(5/2)*sqrt(2*r1 + r3) + 2*r1^9*r3^3 + 2*r1^8*r3^(7/2)*sqrt(2*r1 + r3) - 7*r1^8*r3^4 - 16*r1^7*r3^(9/2)*sqrt(2*r1 + r3) - 24*r1^7*r3^5 - 32*r1^6*r3^(11/2)*sqrt(2*r1 + r3) - 10*r1^6*r3^6 + 12*r1^5*r3^(13/2)*sqrt(2*r1 + r3) + 52*r1^5*r3^7 + 92*r1^4*r3^(15/2)*sqrt(2*r1 + r3) + 102*r1^4*r3^8 + 112*r1^3*r3^(17/2)*sqrt(2*r1 + r3) + 88*r1^3*r3^9 + 64*r1^2*r3^(19/2)*sqrt(2*r1 + r3) + 41*r1^2*r3^10 + 18*r1*r3^(21/2)*sqrt(2*r1 + r3) + 10*r1*r3^11 + 2*r3^(23/2)*sqrt(2*r1 + r3) + r3^12))/(r1^6 + 4*r1^5*sqrt(r3)*sqrt(2*r1 + r3) + 10*r1^5*r3 + 8*r1^4*r3^(3/2)*sqrt(2*r1 + r3) + 19*r1^4*r3^2 + 12*r1^3*r3^3 - 8*r1^2*r3^(7/2)*sqrt(2*r1 + r3) + 3*r1^2*r3^4 - 4*r1*r3^(9/2)*sqrt(2*r1 + r3) + 2*r1*r3^5 + r3^6)
   y3 = sqrt(r3)*sqrt(2*r2 + r3)
   (x1, y1) = (2*r1*sqrt(r3)*sqrt(2*r1 + r3)/(r1 + r3), (-r1^2 + 2*r1*r3 + r1*y3 + r3^2 + r3*y3)/(r1 + r3))
   @printf("天円，人円の直径が %g，%g のとき，地円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  r3 = %g;  r2 = %g;  y3 = %g;  x1 = %g;  y1 = %g\n", r1, r3, r2, y3, x1, y1)
   plot()
   circle(0, y3 + r3 + r1, r1)
   circle2(x1, y1, r1)
   circle2(r2, 0, r2, :blue)
   circle(0, y3, r3, :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, y3 + r3 + r1, "天円:r1,(0,y3+r3+r1)", :red, :center, delta=-delta/2)
       point(x1, y1, "天円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(r2, 0, "地円:r2\n(r2,0)", :blue, :center, delta=-delta/2)
       point(0, y3, "人円:r3\n(0,y3)", :green, :center, delta=-delta/2)
   end
end;

draw(15/2, 7/2, true)

算額（その1248）

算額（その1248）

四十九 群馬県安中市板鼻 鷹巣神社 文政11年(1828)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，楕円2個

等楕円と等円が 2 個，互いに接し合っている。楕円の長径，短径が与えられたとき，等円の直径を求める術を述べよ。

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

include("julia-source.txt");

using SymPy

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

等円の半径は，以下の式のようになる。計算の順序を考えると見た目ほどは複雑ではないといえるかもしれない。

res[1]

なお，当然とも言えるが，a, b は勝手にどのような値でも取れるわけではない。
また，算額では想定していない，円と楕円の接点が別の場所になる解もある。

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x0, y0) = (sqrt(2)*sqrt((a^2 + 3*b^2 - sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(a^2 - b^2))*(3*a^2 - 3*b^2 + sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(8*a), sqrt(2)*a*sqrt((a^2 + 3*b^2 - sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(a^2 - b^2))/2, b*(3*a^2 - 3*b^2 + sqrt(a^4 - 10*a^2*b^2 + 9*b^4))/(2*(a^2 - b^2)))
   @printf("a = %g;  b = %g;  r = %g;  x0 = %g;  y0 = %g\n", a, b, r, x0, y0)
   plot()
   #ellipse(0, b, a, b, color=:red)
   #ellipse(0, -b, a, b, color=:red)
   circle2(r, 0, r, :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(r, 0, " 等円:r,(r,0)", :blue, :left, :vcenter)
       point(0, b, "   等楕円:a,b,(0,b)", :red, :left, :vcenter)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta/2)
   end
   ellipse(0, b, a, b, color=:red)
   ellipse(0, -b, a, b, color=:red)
end;

draw(1, 0.25, true)

算額（その1247）

算額（その1247）

四十九 群馬県安中市板鼻 鷹巣神社 文政11年(1828)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

キーワード：円5個，長方形

長方形の中に，甲，乙，丙，丁，戊の 5 円を容れる。長方形の長辺が 38 寸のとき， 短辺はいかほどか。

長方形の長編，短辺を a, b
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, b - r2)
丙円の半径と中心座標を r3, (r3, r3)
丁円の半径と中心座標を r4, (x4, b - r4)
戊円の半径と中心座標を r5, (a - r5, b - r5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive,
     r3::positive, r4::positive, x4::positive, r5::positive
@syms a, b, r1, r2, r3, r4, x4, r5
a = Sym(38)
eq1 = (a - r1 - r2)^2 + (b - r2 -r1)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (a - r1 - x4)^2 + (r1 - b + r4)^2 - (r1 + r4)^2
eq4 = (r1 - r5)^2 + (b - r5 - r1)^2 - (r1 + r5)^2
eq5 = (r2 - r3)^2 + (b - r2 - r3)^2 - (r2 + r3)^2
eq6 = (x4 - r2)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq7 = (a - r5 - x4)^2 + (r5 - r4)^2 - (r4 + r5)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (b, r1, r2, r3, r4, x4, r5))

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)
   (b, r1, r2, r3, r4, x4, r5) = u
   return [
       -(r1 + r2)^2 + (a - r1 - r2)^2 + (b - r1 - r2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - r3)^2,  # eq2
       -(r1 + r4)^2 + (a - r1 - x4)^2 + (-b + r1 + r4)^2,  # eq3
       (r1 - r5)^2 - (r1 + r5)^2 + (b - r1 - r5)^2,  # eq4
       (r2 - r3)^2 - (r2 + r3)^2 + (b - r2 - r3)^2,  # eq5
       (-r2 + x4)^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq6
       (-r4 + r5)^2 - (r4 + r5)^2 + (a - r5 - x4)^2,  # eq7
   ]
end;

a = 38
iniv = BigFloat[33, 12, 9, 7.5, 4.5, 23, 5.5]
res = nls(H, ini=iniv)

   ([33.00311463837703, 11.814579503859783, 9.106235741252021, 7.437509234825292, 4.897027068931824, 22.461906135655383, 5.325015237401601], true)

長方形の長辺が 38 寸のとき， 短辺は 33.00311463837703 寸である。

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

  a = 38;  b = 33.0031;  r1 = 11.8146;  r2 = 9.10624;  r3 = 7.43751;  r4 = 4.89703;  x4 = 22.4619;  r5 = 5.32502

function draw(a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, r1, r2, r3, r4, x4, r5) = res[1]
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  x4 = %g;  r5 = %g\n",
       a, b, r1, r2, r3, r4, x4, r5)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, b - r2, r2, :green)
   circle(r3, r3, r3, :blue)
   circle(x4, b - r4, r4, :magenta)
   circle(a - r5, b - r5, r5, :orange)
   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 - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-delta)
       point(r2, b - r2, "乙円:r2,(r2,b-r2)", :green, :center, delta=-delta)
       point(r3, r3, "丙円:r3,(r3,r3)", :blue, :center, delta=-delta)
       point(x4, b - r4, "丁円:r4\n(x4,b-r4)", :magenta, :center, delta=-delta)
       point(a - r5, b - r5, "戊円:r5\n(a-r5,b-r5)", :orange, :center, delta=-delta)
       point(a, b, "(a,b)", :black, :right, :bottom, delta=delta)
   end
end;

draw(38, true)