算額（その713）

算額（その713）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード：円1個，楕円3個，正三角形，正六角形

正六角形の中に正三角形と楕円 3 個が入っている。正三角形の中には直径 1 寸の円が入っている。楕円の短径が 0.44 寸のとき，楕円の長径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms r::positive, a::positive, b::positive,
     x0::positive, y0::negative
eq1 = x0^2/a^2 + (y0 + r + b)^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 + r + b)) - 1/√Sym(3)
eq3 = (2r + y0)/x0 - 1/√Sym(3);
solve([eq1, eq2, eq3], (a, x0, y0))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(3)*sqrt(r)*sqrt(-2*b + r), sqrt(3)*r*(2*b - r)/(b - r), r^2/(b - r))

楕円の長半径 a は，r と b の関数である。
sqrt(3r*(r - 2b))
r = 1/2 寸, b = 0.44/2 寸 のとき a = 0.3 である（長径は 0.6 寸）。

r = 1/2
b = 0.44/2
sqrt(3r*(r - 2b))

   0.3

接点の座標は (x0, y0) = (0.185577, -0.892857) である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, b) = (1/2, 0.44/2)
   (a, x0, y0) = (sqrt(3)*sqrt(r)*sqrt(-2*b + r), sqrt(3)*r*(2*b - r)/(b - r), r^2/(b - r))
   @printf("長径 = %g;  x0 = %g;  y0 = %g\n", 2a, x0, y0)
   plot([√3, 0, -√3, √3].*r, [-1, 2, -1, -1].*r, color=:green, lw=0.5)
   plot!([√3, √3, 0, -√3, -√3, 0, √3].*r, [-1, 1, 2, 1, -1, -2, -1].*r, color=:blue, lw=0.5)
   circle(0, 0, r)
   l = r + b
   ellipse(0, -l, a, b)
   ellipse(l*cosd(30), l*sind(30), a, b, φ=120)
   ellipse(l*cosd(150), l*sind(150), a, b, φ=240)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, -r, "-r", :red, :center, :bottom, delta=delta/2)
       point(0, -l, "-r-b")
       point(x0, y0, "(x0,y0)")
   end
end;

算額（その712）

算額（その712）

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

等脚台形の中に正方形と 3 個の等円が入っている。上底の長さが 3 寸のとき，下底の長さはいかほどか。

上底，下底の長さをそれぞれ 2b, 2a とする。
正方形の一辺の長さを 2c
等円の半径と中心座標を r, (0, 2c + r), (c + r, r)
台形の高さは 2c + 2r
とおいて，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     c::positive, r::positive
h = 2c + 2r
eq1 = (c + r)^2 + 4c^2 - 4r^2
eq2 = dist(a, 0, b, h, 0, 2c + r) - r^2
eq3 = dist(a, 0, b, h, c + r, r) - r^2
res = solve([eq1, eq2, eq3], (a, c, r)) 

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (69*b/5, 9*b/5, 3*b)

下底の長さは上底の長さの 69/5 倍，
正方形の一辺の長さは上底の長さの 9/5 倍，
等円の直径は上底の長さの 3 倍である。

b = 3/2
2b .*(69/5, 9/5, 3)

   (41.400000000000006, 5.4, 9.0)

上底の長さが 3 寸のとき，下底の長さは 41.4 寸である。
「答」では 27.6 寸となっている。
「術」では，「上底の 13.8 倍」と書いている。それに従うと 3✕13.8 = 41.4 寸になる。答えは 2✕13.8 = 27.6 寸 と書いてしまったのであろう。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 3//2
   (a, c, r) = (69*b/5, 9*b/5, 3*b)
   @printf("上底の長さ = %g;  下底の長さ = %g;  正方形の一辺の長さ = %g;  等円の直径 = %g\n", 2b, 2a, 2c, 2r)
   h = 2c + 2r
   #@printf("r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g\n", r1, x1, r2, x2, r3)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:green, lw=0.5)
   circle.([1, -1].*(c + r), r, r)
   circle(0, 2c + r, r)
   rect(-c, 0, c, 2c, :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(c, 0, " c", :blue, :left, :bottom, delta=delta/2)
       point(0, 2c, " 2c", :blue, :left, delta=-delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(b, h, " (b,h)", :green, :center, :bottom, delta=delta)
       point(c + r, r, "(c+r,r)", :red, :center, delta=-delta/2)
       point(0, 2c + r, " 2c+r", :red, :left, :vcenter)
   end
end;

算額（その711）

算額（その711）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード：円5個，正方形，直線上
#Julia, #SymPy, #算額, #和算

直線の上に甲円，乙円が 2 個ずつ，丙円が 1 個載っている。丙円の上に正方形が 1 個 載っている。丙円の直径が 1 寸のとき，甲円の直径はいかほどか。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (0, r3)
とおき，以下の連立方程式を解く。
甲円と乙円の中心を結ぶ直線は x 軸と角度 22.5° で交わる。
tan(22.5°) =  tan(PI/8) は，思いの外きれいな数 √2 - 1 である。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive,
     r2::positive, x2::positive,
     r3::positive
eq1 = r2/(x2 + 2r3) - tan(PI/8)
eq2 = r1/(x1 + 2r3) - tan(PI/8)
eq3 = x1 - x2 - 2sqrt(r1*r2)
eq4 = x2 - 2sqrt(r2*r3)
res = solve([eq1, eq2, eq3, eq4], (r1, x1, r2, x2)) 

   1-element Vector{NTuple{4, Sym}}:
    (2*r3*(-4*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 4*sqrt(2 - sqrt(2)) + 7), 2*r3*(-2 + 2*sqrt(4 - 2*sqrt(2)) + 3*sqrt(2)), 2*r3, 2*sqrt(2)*r3)

甲円の半径は 乙円の半径の 2*(-4*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 4*sqrt(2 - sqrt(2)) + 7) = 14 - 8√2 - 4sqrt(4 - 2√2) + 8sqrt(2 - √2) 倍である。

乙円の直径が 1 寸のとき，甲円の直径は 4.479657617687099，4寸4分7厘9毛有奇である。

14 - 8√2 - 4sqrt(4 - 2√2) + 8sqrt(2 - √2)

   4.479657617687099

その他のパラメータは以下のとおりである。正方形の大きさは任意。

r1 = 2.23983;  x1 = 4.40743;  r2 = 1;  x2 = 1.41421;  x3 = 0.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1//2
   a = r3*12 / √2
   (r1, x1, r2, x2) = (2*r3*(-4*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 4*sqrt(2 - sqrt(2)) + 7), 2*r3*(-2 + 2*sqrt(4 - 2*sqrt(2)) + 3*sqrt(2)), 2*r3, 2*sqrt(2)*r3)
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g\n", r1, x1, r2, x2, r3)
   plot([0, 1, 0, -1, 0] .* a, [0, 1, 2, 1, 0] .* a .+ 2r3)
   circle(x1, r1, r1)
   circle(-x1, r1, r1)
   circle(x2, r2, r2, :orange)
   circle(-x2, r2, r2, :orange)
   circle(0, r3, r3, :blue)
   abline(-2r3, 0, tand(45/2), -2.5r3, x1+r1) 
   abline(-2r3, 0, 1, -2.5r3, 1.2x1) 
   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(-2r3, 0, "-2r3", :black, :left, delta=-delta/2)
       point(0, r3, " r3", :blue, :left, delta=-delta/2)
       point(x2, r2, " 乙円:r2\n (x2,r2)", :black, :left, delta=-delta/2)
       point(x1, r1, " 甲円:r1\n (x1,r1)", :red, :left, delta=-delta/2)
   end
end;

算額（その710）

算額（その710）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード：円6個，楕円2個，正方形
#Julia, #SymPy, #算額, #和算

一辺の長さが 2.5 寸の正方形の中に楕円 2 個，楕円の中に大円 2 個，小円 1 個ずつが入っている。小円の直径が 0.4 寸のとき，大円の直径はいかほどか。

正方形の一辺の長さを 2a
大円の半径と中心座標を r1, (x1, a/2)
小円の半径と中心座標を r2, (0, a - r2)
楕円と大円の接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, x1::positive,
     r2::positive,
     x0::positive, y0::positive
b = a/2
eq1 = x1^2 + (a/2 - r2)^2 - (r1 + r2)^2
eq2 = -b^2*x0/(a^2*(y0 - a/2)) + (x0 - x1)/(y0 - a/2)
eq3 = (x0 - x1)^2 + (y0 - a/2)^2 - r1^2
eq4 = x0^2/a^2 + (y0 -a/2)^2/b^2 - 1
res = solve([eq1, eq2, eq3, eq4], (r1, x1, x0, y0)) 

   2-element Vector{NTuple{4, Sym}}:
    (a/2 - r2/2, sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/2, 2*sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/3, a/2 - sqrt(9*a^2 - 24*a*r2 + 12*r2^2)/6)
    (a/2 - r2/2, sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/2, 2*sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/3, a/2 + sqrt(9*a^2 - 24*a*r2 + 12*r2^2)/6)

2 組の解が得られるが，y0 の値が異なる（楕円の長軸に対して対称）だけで，両方とも適解である。

大円の半径は (a - r2)/2なので，直径は a - r2 である。
正方形の一辺の長さが 2.5 寸，小円の直径が 0.4 寸のとき，大円の直径は (2a - 2r2)/2 = (2.5 - 0.4)/2 = 1.05 である。

術は，「正方形の一辺の長さから小円の直径を引き，半分にする」 (2.5 - 0.4)/2 = 1.05

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r2) = (25, 4) .// 20
   b = a/2
   (r1, x1, x0, y0) = (a/2 - r2/2, sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/2, 2*sqrt(3)*sqrt(r2)*sqrt(2*a - r2)/3, a/2 + sqrt(9*a^2 - 24*a*r2 + 12*r2^2)/6)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle4(x1, a/2, r1)
   circle(0, a - r2, r2, :orange)
   circle(0, r2 - a, r2, :orange)
   ellipse(0, a/2, a, b)
   ellipse(0, -a/2, a, b)
   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, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(x1, a/2, " 大円:r1,(x1,a/2)", :red, :center, delta=-delta/2)
       point(0, a - r2, " 小円:r2,(0,a-r2)", :black, :center, :bottom, delta=delta/2)
   end
end;

算額（その709）

算額（その709）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード：円8個，菱形，弦
#Julia, #SymPy, #算額, #和算

外円内に弦を境界として，上部に甲円 4 個，乙円 2 個をいれる。下部には菱形と，菱形に内接する楕円を容れる。
乙円の径を▢寸，楕円の短径を一寸五分としたとき，楕円の長径はいかほどか。

注：欠損した一文字は「答」から推測すると「一」であろう。なお，「答」にも「長径▢寸七分五厘有奇」と欠損文字があるが，こちらは「二」であろう。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (x1, R - 2r1); x1 = √3r1
乙円の半径と中心座標を r2, (x2, ya + r2); ya = R - 4r1 = 2r1 - R
楕円の長半径と短半径を a, b
楕円と菱形の接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::positive,
     a::positive, b::positive,
     x0::positive, y0::positive

x1 = sqrt(Sym(3))r1
yb = r1 - R
xb = sqrt(R^2 - yb^2)
eq1 = x1^2 + (R - 2r1)^2 - (R - r1)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = (x1 - x2)^2 + (2r1 - r2)^2 - (r1 + r2)^2
eq4 = x0^2/a^2 + (y0 - r1 + R)^2/b^2 - 1
eq5 = -b^2*x0/(a^2*(y0 - r1 + R)) + r1/xb
eq6 = (y0 - yb)/(xb - x0) - r1/xb;

図の上半分は R, r1, x2 がわかれば描ける。

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

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-2*sqrt(2)*r2 + 3*r2, r2*(3 - 2*sqrt(2))/3, -2*sqrt(6)*r2/3 + 2*sqrt(3)*r2/3)
    (2*sqrt(2)*r2 + 3*r2, r2*(2*sqrt(2) + 3)/3, 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)

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

res[2][1] |> factor |> println
res[2][2] |> factor |> println
res[2][3] |> factor |> println

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

r2 = 1/2
(2*sqrt(2)*r2 + 3*r2, r2*(2*sqrt(2) + 3)/3, 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)

   (2.914213562373095, 0.9714045207910317, 1.3938468501173515)

下半分も描くには SymPy は力不足なので，数値解を求める。

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

r2 = 1//2
b = 15//20
iniv = BigFloat[2.91, 0.97, 1.39, 0.88, -1.36, 1.38]
res = nls(H, ini=iniv)

   ([2.914213562373095, 0.9714045207910317, 1.3938468501173518, 0.8773124760905445, -1.363750587600455, 1.3804469258418701], true)

楕円の長径は 2.7608938516837402 （長半径は 1.3804469258418701）である。
「答」では「長径は▢寸七分5厘有奇」となっている。

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

   R = 2.91421;  r1 = 0.971405;  x2 = 1.39385;  ya = -0.971405;  x0 = 0.877312;  y0 = -1.36375;  a = 1.38045

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, b) = (1/2, 1.5/2)
   (R, r1, x2, x0, y0, a) = res[1]
   x1 = √3r1
   ya = R - 4r1  # = 2r1 - R
   xa = sqrt(R^2 - ya^2)
   yb = r1 - R
   xb = sqrt(R^2 - yb^2)
   println("長径 = $(2a) （長半径:a = $a）")
   @printf("r2 = %g;  b = %g\n", r2, b)
   @printf("R = %g;  r1 = %g;  x2 = %g;  ya = %g;  x0 = %g;  y0 = %g;  a = %g\n",
       R, r1, x2, ya, x0, y0, a)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, R - r1, r1)
   circle(0, R - 3r1, r1)
   circle(x1, R - 2r1, r1)
   circle(-x1, R - 2r1, r1)
   circle(0, r1 - R, r1, :gray90)
   circle(x2, ya + r2, r2, :orange) 
   circle(-x2, ya + r2, r2, :orange) 
   segment(-xa, ya, xa, ya)
   plot!([xb, 0, -xb, 0, xb], [yb, ya, yb, -R, yb], color=:green, lw=0.5)
   ellipse(0, r1 - R, a, b)
   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 - r1, " 甲円:r1\n (0,R-r1)", :red, :left, :vcenter)
       point(0, R - 3r1, "", :red)
       point(x1, R - 2r1, "(x1,R-2r1)", :red, :center, delta=-delta/2)
       point(x2, ya + r2, " 乙円:r2\n (x2,ya+r2)", :black, :left, :vcenter)
       point(0, ya, " ya=R-4r1", :black, :center, :bottom, delta=delta)
       segment(-xb, yb, xb, yb)
       segment(0, ya, 0, -R)
       point(0, r1 - R, "r1-R ", :black, :right, :bottom, delta=delta)
       point(a, r1 - R, "(a,r1-R) ", :red, :right, :bottom, delta=delta)
       point(0, r1 - R + b, " r1-R+b", :red, :left, delta=-delta)
   end
end;

算額（その708）

算額（その708）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード：円13個，外円，弦2本
#Julia, #SymPy, #算額, #和算

外円内に水平な 2 本の弦，甲円 1 個，乙円，丙円，丁円を 2 個ずつ，戊円 4 個，己円 1 個を容れる。
己円の直径が 33 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 6r5 - R)
乙円の半径と中心座標を r2, (x2, 2r5 - R)
丙円の半径と中心座標を r3, (x3, y + r3)
丁円の半径と中心座標を r4, (x4, R - 2r6 - r4)
戊円の半径と中心座標を r5, (0, r5 - R), (0, 3r5 - R), (0, 5r5 - R), (0, 7r5 - R)
弦と y 軸の交点座標を (0, y), (0, R - 2r6)
その他の変数の関係は
R = r6 + 4r5
y = 4r5 - R
r1 = 2r5
y4 = R - 2r6 - r4
とおき，以下の連立方程式を(r6 を変数のまま）解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive,
     r4::positive, x4::positive, y4::positive,
     r5::positive, r6::positive
@syms r1, r2, x2, r3, x3, r4, x4, y4, r5, r6
R = r6 + 4r5
y = 4r5 - R
y4 = R - 2r6 - r4
r1 = 2r5
eq1 = x2^2 + (2r5 - R)^2 - (R - r2)^2
eq2 = x3^2 + (y + r3)^2 - (R - r3)^2
eq3 = x4^2 + y4^2 - (R - r4)^2
eq4 = x2^2 + r5^2 - (r2 + r5)^2
eq5 = x3^2 + (6r5 - R - y - r3)^2 - (r1 + r3)^2
eq6 = x4^2 + (6r5 - R - y - r4)^2 - (r1 + r4)^2
eq7 = (x3 - x4)^2 + (y4 - y - r3)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r2, x2, r3, x3, r4, x4, r5));

10 組の解が得られるが，最初の 2 つが適解である。ただし，x2, x3, x4 は符号が負の解が得られる（y 軸で線対称なので適解である）。

res[1]

   (28*r6/11, -2*sqrt(21)*sqrt(r6*(143*r6 + 73*sqrt(r6^2)))/33, r6*(10757*r6 + 11023*sqrt(r6^2))/(66*(59*r6 + 73*sqrt(r6^2))), -sqrt(2)*sqrt(r6^3*(10757*r6 + 11023*sqrt(r6^2)))/(33*r6), 2*r6*(59*r6 + 73*sqrt(r6^2))/(92*r6 + 73*sqrt(r6^2)), -2*sqrt(66)*sqrt(r6/(92*r6 + 73*sqrt(r6^2)))*(59*r6 + 73*sqrt(r6^2))/33, 59*r6/66 + 73*sqrt(r6^2)/66)

res[2]

   (28*r6/11, 2*sqrt(21)*sqrt(r6*(143*r6 + 73*sqrt(r6^2)))/33, r6*(10757*r6 + 11023*sqrt(r6^2))/(66*(59*r6 + 73*sqrt(r6^2))), -sqrt(2)*sqrt(r6^3*(10757*r6 + 11023*sqrt(r6^2)))/(33*r6), 2*r6*(59*r6 + 73*sqrt(r6^2))/(92*r6 + 73*sqrt(r6^2)), -2*sqrt(66)*sqrt(r6/(92*r6 + 73*sqrt(r6^2)))*(59*r6 + 73*sqrt(r6^2))/33, 59*r6/66 + 73*sqrt(r6^2)/66)

乙円の半径(r2)は 28*r6/11 となり，己円の半径(r6)の 28/11 倍である。術でも「己径二十八段以十一個除（己円の直径を28倍して11で割る）」としている。
したがって，己円の直径が 33 寸のとき，乙円の直径は 84 寸である。

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

   R = 148.5;  r1 = 66;  r2 = 42;  x2 = -67.3498;  r3 = 41.25;  x3 = -104.355;  r4 = 26.4;  x4 = -83.4841;  y4 = 89.1;  r5 = 33;  r6 = 16.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #(r2, x2, r3, x3, r4, x4, r5) =
   #    (42, 68, 41, 104, 26, 84, 33)
   r6 = 33//2
   (r2, x2, r3, x3, r4, x4, r5)= (28*r6/11, -2*sqrt(21)*sqrt(r6*(143*r6 + 73*sqrt(r6^2)))/33, r6*(10757*r6 + 11023*sqrt(r6^2))/(66*(59*r6 + 73*sqrt(r6^2))), -sqrt(2)*sqrt(r6^3*(10757*r6 + 11023*sqrt(r6^2)))/(33*r6), 2*r6*(59*r6 + 73*sqrt(r6^2))/(92*r6 + 73*sqrt(r6^2)), -2*sqrt(66)*sqrt(r6/(92*r6 + 73*sqrt(r6^2)))*(59*r6 + 73*sqrt(r6^2))/33, 59*r6/66 + 73*sqrt(r6^2)/66)
   R = r6 + 4r5
   y = 4r5 - R
   r1 = 2r5
   y4 = R - 2r6 - r4
   @printf("乙円の直径 = %g;  己円の直径 = %g\n", 2r2, 2r6)
   @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  r6 = %g\n",
       R, r1, r2, x2, r3, x3, r4, x4, y4, r5, r6)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, 6r5 - R, r1, :brown)
   circle(x2, 2r5 - R, r2, :magenta)
   circle(-x2, 2r5 - R, r2, :magenta)
   circle(x3, y + r3, r3, :green)
   circle(-x3, y + r3, r3, :green)
   circle(x4, y4, r4, :orange)
   circle(-x4, y4, r4, :orange)
   x1 = sqrt(R^2 - y^2)
   segment(-x1, y, x1, y, :green)
   circle(0, 7r5 - R, r5)
   circle(0, 5r5 - R, r5)
   circle(0, 3r5 - R, r5)
   circle(0, r5 - R, r5)
   y20 = R - 2r6
   x20 = sqrt(R^2 - y20^2)
   segment(-x20, y20, x20, y20, :green)
   circle(0, R - r6, r6, :gray)
   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, 6r5 - R, "甲円:r1\n(0,6r5-R)", :brown, :center, delta=-delta)
       point(-x2, 2r5 - R, "乙円:r2\n(x2,2r5-R)", :magenta, :center, delta=-delta/2)
       point(-x3, y + r3, "丙円:r3\n(x3,y+r3)", :green, :center, delta=-delta/2)
       point(-x4, y4, "丁円:r4\n(x4,y4)", :orange, :center, delta=-delta/2)
       point(0, r5 - R, "戊円:r5\n(0,r5-R)", :red, :center, delta=-delta/2)
       point.(0, [3, 5, 7].*r5 .- R, "", :red)
       point(0, R - r6, "己円:r6,(0,R-r6)", :black, :center, :bottom, delta=delta/2)
       point(0, y, " y", :green, :left, delta=-delta/2)
       point(0, R - 2r6, " R-2r6", :green, :left, :bottom, delta=delta/2)
   end
end;