算額（その710）

算額（その710）

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

一辺の長さが 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）

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

外円内に弦を境界として，上部に甲円 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）

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

外円内に水平な 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;

算額（その707）

算額（その707）

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

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

直径が 2.5 寸の外円内に弦と斜線を入れ，区画された領域に等円を 4 個入れる。等円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
弦と外円の交点座標を (x1, y1); y1 = 2r - R, x1 = sqrt(R^2 - y1^2)
斜線と外円の交点座標を (x2, y2); x2 = sqrt(R^2 - y2^2)
等円の半径と中心座標を r, (0, r - R), (0, R - r), (x, y1 + r)
とおき以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     x1::positive, y1::positive,
     x2::positive, y2::positive
y1 = 2r - R
x1 = sqrt(R^2 - y1^2)
x2 = sqrt(R^2 - y2^2)
eq1 = dist(0, y1, x2, y2, 0, R - r) - r^2
eq2 = dist(0, y1, x2, y2, x, y1 + r) - r^2
eq3 = x^2 + (y1 + r)^2 - (R - r)^2;
# res = solve([eq1, eq2, eq3], (r, x, y2))

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, x, y2) = u
   return [
       -r^2 + (2*R - 3*r)^2*(R^2 - y2^2)*(R - 2*r + y2)^2/(R^2 - y2^2 + (R - 2*r + y2)^2)^2 + (2*R - 3*r - (2*R - 3*r)*(R - 2*r + y2)^2/(R^2 - y2^2 + (R - 2*r + y2)^2))^2,  # eq1
       -r^2 + (r - (r*(R - 2*r + y2) + x*sqrt(R^2 - y2^2))*(R - 2*r + y2)/(R^2 - y2^2 + (R - 2*r + y2)^2))^2 + (x - sqrt(R^2 - y2^2)*(r*(R - 2*r + y2) + x*sqrt(R^2 - y2^2))/(R^2 - y2^2 + (R - 2*r + y2)^2))^2,  # eq2
       x^2 + (-R + 3*r)^2 - (R - r)^2,  # eq3
   ]
end;

R = 25//20
iniv = BigFloat[0.5, 0.8, 1.1]
res = nls(H, ini=iniv)

   ([0.4734152343522919, 0.7576940667754868, 1.0587960931595892], true)

外円の直径が 2.5 寸のとき，等円の直径は 2*0.4734152343522919 = 0.9468304687045838 寸である。

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

r = 0.473415;  x = 0.757694;  x1 = 1.21268;  y1 = -0.30317;  x2 = 0.664418;  y2 = 1.0588

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 25//20
   (r, x, y2) = res[1]
   y1 = 2r - R
   x1 = sqrt(R^2 - y1^2)
   x2 = sqrt(R^2 - y2^2)
   @printf("等円の直径 = %g\n", 2r)
   @printf("r = %g;  x = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n", r, x, x1, y1, x2, y2)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, r - R, r)
   circle(x, y1 + r, r)
   circle(-x, y1 + r, r)
   circle(0, R - r, r)
   segment(-x1, y1, x1, y1, :green)
   segment(0, y1, x2, y2, :green)
   segment(0, y1, -x2, y2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x, y1 + r, " (x,y1+r)", :red, :left, :vcenter)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(0, r - R, " r-R", :red, :left, :vcenter)
       point(0, y1, " y1", :green, :left, delta=-delta/2)
       point(x1, y1, "(x1,y1) ", :green, :right, delta=-delta/2)
       point(x2, y2, "(x2y2)", :green, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その706）

算額（その706）

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

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

キーワード：円5個，外円，二等辺三角形

直径 163 寸の外円内に圭（二等辺三角形），大円 2 個，小円 2 個が入っている。小円の直径はいかほどか。
注：後述するが 163 寸ではなく 162 寸である。

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

include("julia-source.txt");

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

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*R/9, 20*R/81, 8*sqrt(5)*R/27)

小円の半径は 外円の半径の 20/81 倍である。

術では「外径を 20 倍し 80 で割る」とあるが，正しくは「81 で割る」であろう。
また，問にある「外径163寸」というのも「外径162寸」であろう。
この問に限らず，この算額にはミスが多い。

以上 2 点を修正すると，「外円の直径が 162 寸のとき，小円の直径は 40 寸」である。

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

r1 = 36;  r2 = 20;  x2 = 53.6656

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 162//2  # 「問」のミスを修正
   (r1, r2, x2) = (4*R/9, 20*R/81, 8*sqrt(5)*R/27)
   @printf("小円の直径 = %g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g\n", r1, r2, x2)
   y1 = 2r1 - R
   x1 = sqrt(R^2 - y1^2)
   plot([x1, 0, -x1, x1], [y1, R, y1, y1], color=:black, lw=0.5)
   circle(0, 0 , R, :blue)
   circle(0, r1 - R, r1)
   circle(0, 3r1 - R, r1)
   circle(x2, 2r1 - R - r2, r2, :blue)
   circle(-x2, 2r1 - R - r2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1 - R, " 大円:r1\n (0,r1-R)", :red, :left, :vcenter)
       point(0, 3r1 - R, " 大円:r1\n (0,3r1-R)", :red, :left, :vcenter)
       point(x2, 2r1 - R - r2, "小円:r2\n(x2,2r1-R-r2)", :blue, :center, delta=-delta/2)
       point(x1, y1, "(x1,y1)  ", :black, :right, :bottom, delta=delta/2)
   end
end;

 

算額（その705）

算額（その705）

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

外円内に等斜（この場合は弦と矢），圭（二等辺三角形），等円 3 個をいれる。等斜が 1 寸のとき，等円の直径はいかほどか。

等斜の長さを a，圭の底辺の長さを 2b とする。
外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, a - R + r), (x, a - R + r)
として，いかの連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive, r::positive, R::positive
eq1 = x^2 + (a - R + r)^2 - (R - r)^2
eq2 = (a - R)^2 + (a/2)^2 - R^2
eq3 = dist(0, R, b, a - R, x, a - R + r) - r^2
eq4 = dist(0, R, b, a - R, 0, a - R + r) - r^2
solve([eq1, eq2, eq3, eq4], (r, x, b, R))

   1-element Vector{NTuple{4, Sym}}:
    (a*(3 - sqrt(5))/8, a*sqrt(-2 + sqrt(5))/2, a*(3 - sqrt(5))/(8*sqrt(-2 + sqrt(5))), 5*a/8)

等円の半径は a*(3 - sqrt(5))/8 である。

等斜を 1 寸とすれば，等円の半径は 0.09549150281252627 寸，直径は 0.19098300562505255 寸である。

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

r = 0.0954915;  x = 0.242934;  b = 0.196538;  R = 0.625

a = 1
a*(3 - sqrt(5))/8, 2*a*(3 - sqrt(5))/8

   (0.09549150281252627, 0.19098300562505255)

「答曰 等円径三分八厘二毛」とあるが，これは直径を二倍してしまっているようだ。

「術曰」では，「置五個開平方以滅▢個内餘乗等斜四除之得等径合問」とある。
欠損している文字は「三」で，「三から五の平方根を引き等斜を掛け四で割る」と読める。上述の式の通り a*(3 - sqrt(5))/4 ということで，直径を表している。これで1分9厘0毛9糸なので，それをまた2倍して3分8厘2毛としてしまっている。

a*(3 - sqrt(5))/4, 2*a*(3 - sqrt(5))/4

   (0.19098300562505255, 0.3819660112501051)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (r, x, b, R) = (a*(3 - sqrt(5))/8, a*sqrt(-2 + sqrt(5))/2, a*(3 - sqrt(5))/(8*sqrt(-2 + sqrt(5))), 5*a/8)
   println("等円の直径 = $(2r)")
   @printf("r = %g;  x = %g;  b = %g;  R = %g\n", r, x, b,R)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, a - R + r, r)
   circle(x, a - R + r, r)
   circle(-x, a - R + r, r)
   segment(0, R, b, a - R)
   segment(0, R, -b, a - R)
   segment(-a/2, sqrt(R^2 - a^2/4), a/2, sqrt(R^2 - a^2/4))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, a - R, " a-R", :black, :left, delta=-delta/2)
       point(b, a - R, "b", :black, :center, delta=-delta/2)
       point(x, a - R + r, " r,(x,a-R+r)", :black, :left, :vcenter)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その704）

算額（その704）

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

一辺の長さが 3.5 寸の正方形内に 4 つの四分円を描く。黒積（図に色付した部分の面積の 8 倍）を求めよ。

この問題は，算額（その122）と同趣旨のものである。

正方形の一辺の長さを a とする。

図に色付した部分の面積は，正方形の面積の 1/2 から，半径 a の円の面積の 1/12 と底辺の長さ a/2，高さ √3a/2 の三角形の面積を差し引いたものである（ACDO - ABO - BDO）。
黒積はその 8 倍である。

include("julia-source.txt");

a = 3.5
8(a^2/2 - pi*a^2/12 - √3*a^2/8)

   2.1260376029646118

SymPy の integrate() でも求めてみる。

using SymPy
@syms a, x
a = 3.5
8*integrate(a - sqrt(a^2 - x^2), (x, 0, a/2)) |> println

   2.12603760296460

黒積は約 2.126 平方寸である。

術では 2.195 としているが，円周率として 3.13314827844261 という不思議な値を使っていることになる。

@syms 円周率
eq = 8(a^2/2 - 円周率*a^2/12 - √3*a^2/8) - 2.195
solve(eq, 円周率)[1] |> println

   3.13314827844261

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3.5
   θ = 60:0.01:90
   x = vcat(a .* cosd.(θ), [a/2, a/2])
   y = vcat(a .* sind.(θ), [a, √3a/2])
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(0, 0 , a, beginangle=0, endangle=90)
   circle(a, 0 , a, beginangle=90, endangle=180)
   circle(a, a , a, beginangle=180, endangle=270)
   circle(0, a , a, beginangle=270, endangle=360)
   segment(0, a/2, a, a/2)
   segment(a/2, 0, a/2, a)
   segment(0, 0, a/2, √3a/2)
   plot!(x, y, seriestype=:shape, color=:magenta, fillcolor=:magenta, alpha=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, a, " A", :green, :left, :bottom, delta=delta/2)
       point(a/2, √3*a/2, " B", :green, :left, :bottom, delta=1.5delta)
       point(a/2, a, " C", :green, :center, :bottom, delta=delta/2)
       point(a/2, 0, " D", :green, :left, :bottom, delta=delta/2)
       point(0, 0, "  O", :green, :left, :bottom, delta=delta/4)
   end
end;

算額（その703）

算額（その703）

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

大円 2 個と小円 1 個が直線上に載っている。小円に接する斜線を引き，大円 2 個の弦とする。弦の中心と弧の中心を「矢」というが，直径から矢の長さを引いたものを甲矢，乙矢とする。甲矢，乙矢の長さが与えられたとき，大円の直径を求めよ。

大円と小円が互いに接しまた直線に接していることから，
r1^2 + (r1 - r2)^2 = (r1 + r2) より，r2 = r1/4 である。

大円の中心から斜線までの距離は，甲矢，乙矢の長さから半径を引いたものと一致する。

斜線の x 切片，y 切片を (a, 0), (0, b)
大円の半径と中心座標を r1, (r1, r1), (-r1, r1)
小円の半径と中心座標を r2, (0, r2)
甲矢，乙矢を「甲矢」，「乙矢」
として，以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms 甲矢::positive, 乙矢::positive,
     a::positive, b::positive,
     r1::positive, r2::positive
r2 = r1/4
eq1 = dist(a, 0, 0,  b, r1, r1) - (甲矢 - r1)^2
eq2 = dist(a, 0, 0,  b, -r1, r1) - (乙矢 - r1)^2
eq3 = dist(a, 0, 0,  b, 0, r2) - r2^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (r1, a, b) = u
   return [
       -(-r1 + 甲矢)^2 + (-b*(-a*(-a + r1) + b*r1)/(a^2 + b^2) + r1)^2 + (-a + a*(-a*(-a + r1) + b*r1)/(a^2 + b^2) + r1)^2,  # eq1
       -(-r1 + 乙矢)^2 + (-b*(-a*(-a - r1) + b*r1)/(a^2 + b^2) + r1)^2 + (-a + a*(-a*(-a - r1) + b*r1)/(a^2 + b^2) - r1)^2,  # eq2
       -r1^2/16 + (-a + a*(a^2 + b*r1/4)/(a^2 + b^2))^2 + (-b*(a^2 + b*r1/4)/(a^2 + b^2) + r1/4)^2,  # eq3
   ]
end;

(甲矢, 乙矢) = (15, 10)
iniv = BigFloat[9, 12, 4.7]
res = nls(H, ini=iniv)

   ([8.520833333333334, 14.20138888888889, 4.35848252344416], true)

甲矢 15 寸，乙矢 10 寸のとき，大円の直径は 17.0417 寸である。

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

r1 = 8.52083;  r2 = 2.13021;  a = 14.2014;  b = 4.35848

function 矢(x0, y0, r, a, b)
   x1 = (a*y0 + x0 - a*b)/(a^2 + 1)
   y1 = a*x1 + b
   (x2, y2) = (a*y0 - a*(r/sqrt(a^2 + 1) + y0) + x0, r/sqrt(a^2 + 1) + y0)
   return (x1, y1, x2, y2)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, a, b) = res[1]
   甲矢 = r1 + sqrt(dist(a, 0, 0,  b, r1, r1))
   乙矢 = r1 + sqrt(dist(a, 0, 0,  b, -r1, r1))
   r2 = r1/4
   @printf("甲矢 = %g;  乙矢 = %g;  大円の直径 = %g\n", 甲矢, 乙矢, 2r1)
   @printf("r1 = %g;  r2 = %g;  a = %g;  b = %g\n", r1, r2, a, b)
   plot()
   circle(r1, r1, r1)
   circle(-r1, r1, r1)
   circle(0, r2, r2, :blue)
   abline(0, b, -b/a, -2r1, 2r1)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :top, delta=-2delta)
       point(-r1, r1, " (-r1,r1)", :red, :left, :vcenter)
       point(r1, r1, " (r1,r1)", :red, :left, :vcenter)
       point(0, r2, " r2", :blue, :left, :vcenter)
       (x1, y1, x2, y2) = 矢(r1, r1, r1, -b/a, b)
       segment(x1, y1, x2, y2, :green, lw=1)
       point(10, 12, "甲矢",  mark=false)
       (x3, y3, x4, y4) = 矢(-r1, r1, r1, -b/a, b)
       segment(x3, y3, x4, y4, :green, lw=1)
       point(-10.5, 12, "乙矢",  mark=false)
   end
end;

 

算額（その702）

算額（その702）

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

直径が 25 寸の外円の中に，菱形 1 個，甲円 2 個，乙円 8 個が入っている。左右にある 2 個の乙円は甲円に内接している。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (R/2, 0)
第一象限にある 2 個の乙円の半径と中心座標を r2, (x2, r2), (r2, y2)
菱形の短い方の対角線（菱平）の長さを b とする。菱形の長い方の対角線は外円の周上にある。
第一象限の上方にある中心座標が (r2, y2) の乙円は菱形の斜辺に接している。

この問題は，条件が足りないように思える。
そもそも，答えでは「菱平は四寸四分六厘七毛有奇」とあるが，外円の半径が12.5寸なのに，菱平が 4.467 寸ほどしかないような図は存在しないない。
算額でよくあるように，答えがきれいな数になるような条件を探索すると，乙円の直径が2寸のときに菱平が 22.008990377897288 になることがわかる。

このときの図は以下のようになる。

よって，この問題は条件，「乙円の直径が2寸のとき」を追加して解く。

なお，術にも数文字の欠損が見られ，問・答・術を通じて真実はわからない。

include("julia-source.txt");

using SymPy
@syms R::positive, b::positive, r1::positive,
     r2::positive, x2::positive, y2::positive
R = 25//2
r1 = R/2
r2 = 2//2
eq1 = (x2 - r1)^2 + r2^2 - (r1 - r2)^2
eq2 = dist(0, b, R, 0, r2, y2) - r2^2
eq3 = r2^2 + y2^2 - (R - r2)^2
eq4 = r2/(R - x2) - b/R
res = solve([eq1, eq2, eq3], (x2, y2, b))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (25/4 - 5*sqrt(17)/4, 5*sqrt(42)*sqrt(537 - 92*sqrt(2))/521 + 115*sqrt(21)*sqrt(537 - 92*sqrt(2))/1042, 5*sqrt(179/28 - 23*sqrt(2)/21))
    (25/4 - 5*sqrt(17)/4, -5*sqrt(42)*sqrt(92*sqrt(2) + 537)/521 + 115*sqrt(21)*sqrt(92*sqrt(2) + 537)/1042, 5*sqrt(23*sqrt(2)/21 + 179/28))
    (5*sqrt(17)/4 + 25/4, 5*sqrt(42)*sqrt(537 - 92*sqrt(2))/521 + 115*sqrt(21)*sqrt(537 - 92*sqrt(2))/1042, 5*sqrt(179/28 - 23*sqrt(2)/21))
    (5*sqrt(17)/4 + 25/4, -5*sqrt(42)*sqrt(92*sqrt(2) + 537)/521 + 115*sqrt(21)*sqrt(92*sqrt(2) + 537)/1042, 5*sqrt(23*sqrt(2)/21 + 179/28))

4組の解が得られるが，3 番目のものが適解である。
これに従うと，菱平は22寸有奇ということになる。

res[3][3] |> sympy.sqrtdenest |> println
res[3][3].evalf()*2 |> println

   5*sqrt(21)*(23 - 2*sqrt(2))/42
   22.0089903778973

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

R = 12.5;  r1 = 6.25;  r2 = 1;  x2 = 11.4039; y2 = 11.4564, b = 11.0045

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 25//2
   r1 = R/2
   r2 = 2//2  # 追加条件
   (x2, y2, b) = (5*sqrt(17)/4 + 25/4, 5*sqrt(42)*sqrt(537 - 92*sqrt(2))/521 + 115*sqrt(21)*sqrt(537 - 92*sqrt(2))/1042, 5*sqrt(179/28 - 23*sqrt(2)/21))
   println("菱平（短い方の対角線） = $(2b)")
   @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g; y2 = %g, b = %g\n", R, r1, r2, x2, y2, b)
   plot([R, 0, -R, 0, R], [0, b, 0, -b, 0], color=:black, lw=0.5)
   circle(0, 0, R, :magenta)
   circle4(r2, y2, r2)
   circle4(x2, r2, r2)
   circle(r1, 0, r1, :blue)
   circle(-r1, 0, r1, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r2, x2, " 乙円:r2,(r2,y2)", :black, :left, :vcenter)
       point(x2, r2, " 乙円:r2,(x2,r2)", :black, :left, :vcenter)
       point(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(r1, 0, "甲円:r1\n(r1,0)", :blue, :center, delta=-delta/2)
       point(0, b, " b", :black, :left, :top, delta=-2delta)
   end
end;

 

算額（その701）

算額（その701）

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

長方形の中に甲円 4 個，乙円 6 個が入っている。
乙円の直径が▢寸のとき，菱形の短い方の対角線の長さはいかほどか。
注：乙円の直径の数が欠損しているが，「一」である。

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

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive,
     r2::positive, x2::positive
r2 = 1//2
eq1 = r2/x2 - r1/(a - r1)
eq2 = (r1 - r2)^2 + (b - r1)^2 - (r1 + r2)^2
eq3 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = -2*a*r2*x2 - b*r2^2 + b*x2^2 + 4*r2^2*x2 # 甲円の中心から菱形の辺までの距離（の分子）

res = solve([eq1, eq2, eq3, eq4], (a, b, r1, x2));
res[4]

   (sqrt(6)/3 + 3, 5/2, 7/2 - sqrt(6), 1/2 + sqrt(6)/3)

4 組の解が得られるが，最後のものが適解である。

菱形の短い方の対角線の長さは長方形の短辺の長さに等しく，5 寸である。
術でも，「乙円の径を5倍する」となっており，問で欠けている乙円の直径の数値は 1 であることがわかる。

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

a = 3.8165;  b = 2.5;  r1 = 1.05051;  x2 = 1.3165

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (a, b, r1, x2) = (sqrt(6)/3 + 3, 5/2, 7/2 - sqrt(6), 1/2 + sqrt(6)/3)
   @printf("菱形の短い方の対角線の長さ = %g\n", 2b)
   @printf("a = %g;  b = %g;  r1 = %g;  x2 = %g\n", a, b, r1, x2)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:magenta, lw=0.5)
   plot!([a - 2r2, 0, 2r2 - a, 0, a - 2r2], [0, b, 0, -b, 0], color=:blue, lw=0.5)
   circle4(a - r1, b - r1, r1, :orange)
   circle4(x2, b - r2, r2)
   circle(a - r2, 0, r2)
   circle(r2 - a, 0, r2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(a - 2r2, 0, "a-2r2  ", :blue, :right, :bottom, delta=delta/2)
       point(a - r2, 0, "a-r2", :blue, :center, delta=-delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a - r1, b - r1, "甲円:r1\n(a-r1,b-r1)", :black, :center, delta=-2delta)
       point(x2, b - r2, " 乙円:r2,(x2,b-r2)", :red, :left, :vcenter)
   end
end;

算額（その700）

算額（その700）

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

正方形の中に四分円 2 個，大円，小円それぞれ 2 個ずつ入っている。小円の直径が 1 寸のとき，大円の直径はいかほどか。

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

include("julia-source.txt");

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

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

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

大円の半径は，小円の半径の (6√3 - 9) 倍である。
小円の直径が 1 寸の場合，大円の直径は 6√3 - 9 = 1.392304845413264 寸である。

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

a = 1.5;  r1 = 0.696152;  y2 = 2.44949

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (a, r1, y2) = (3*r2, 3*r2*(-3 + 2*sqrt(3)), 2*sqrt(6)*r2)
   @printf("大円の直径 = %g;  a = %g;  r1 = %g;  y2 = %g\n", 2r1, a, r1, y2)
   plot([a, a, -a, -a, a], [0, 2a, 2a, 0, 0], color=:magenta, lw=0.5)
   circle(a, 0, 2a, :blue, beginangle=90, endangle=180)
   circle(-a, 0, 2a, :blue, beginangle=0, endangle=90)
   circle(r1, r1, r1)
   circle(-r1, r1, r1)
   circle(a - r2, y2, r2, :green)
   circle(r2 - a, y2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
       point(0, 2a, " 2a", :magenta, :left, :bottom, delta=delta/2)
       point(r1, r1, "大円:r1,(r1, r1)", :red, :center, delta=-delta/2)
       point(a - r2, y2, "小円:r2,(a-r2,y2)", :green, :center, delta=-delta/2)
   end
end;

 

算額（その699）

算額（その699）

山口正義: やまぶき4，第58号，2018/12/06.
千葉県君津市鹿野山 鹿野山神野寺 万延二年(1861)
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf
山口正義：やまぶき，和算と歴史随想 からリンク
https://yamabukiwasan.sakura.ne.jp/page3.html

高さ 64 寸の正三角形内に，短径 6 寸の甲楕円，長径 48 寸の乙楕円が入っている。
甲楕円の長径，乙楕円の短径を求めよ。

正三角形の一辺の長さを 2a とする。
甲楕円の長半径，短半径，中心座標を a1, b1, (0, 2b2 + b1)
乙楕円の長半径，短半径，中心座標を a2, b2, (0, b2)
正三角形の斜辺と甲楕円，乙楕円との接点座標を (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a1::positive, b1::positive, x1::positive, y1::positive,
     a2::positive, b2::positive, x2::positive, y2::positive
@syms a::positive

s3 = sqrt(Sym(3))
(a1, a2) = (6, 48) .// 2
a = 64/s3
y1 = (a - x1)*s3
y2 = (a - x2)*s3
eq1 = x1^2/a1^2 + (y1 - (2b2 + b1))^2/b1^2 - 1
eq3 = x2^2/a2^2 + (y2 - b2)^2/b2^2 - 1
eq2 = -b1^2*x1/(a1^2*(y1 - (2b2 + b1))) + s3
eq4 = -b2^2*x2/(a2^2*(y2 - b2)) + s3
res = solve([eq1, eq2, eq3, eq4], (b1, x1, b2, x2))

   1-element Vector{NTuple{4, Sym}}:
    (13, 9*sqrt(3)/14, 37/2, 1152*sqrt(3)/91)

甲楕円の長径は 26 寸，乙楕円の短径は 37 寸である。

include("julia-source.txt");
function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   s3 = sqrt(3)
   h = 64
   a = h/s3
   (a1, a2) = (6, 48) .// 2
   (b1, x1, b2, x2) = (13, 9*sqrt(3)/14, 37/2, 1152*sqrt(3)/91)
   @printf("甲楕円の長径 = %g;  乙楕円の短径 = %g\n", 2b1, 2b2)
   y1 = (a - x1)*s3
   y2 = (a - x2)*s3
   plot([a, 0, -a, a], [0, h, 0, 0], color=:magenta, lw=0.5)
   ellipse(0, b2, a2, b2, color=:blue)
   ellipse(0, b1 + 2b2, a1, b1, color=:red)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, b1 + 2b2, " 甲楕円:a1,b1,(0,2b2+b1)", :black, :left, :vcenter)
       point(0, b2, " 乙楕円:a2,b2,(0,b2)", :black, :left, :vcenter)
       point(x1, y1, " (x1,y1)", :green, :left, :vcenter)
       point(x2, y2, " (x2,y2)", :green, :left, :vcenter)
   end
end;

 

算額（その698）

算額（その698）

新潟県長岡市上岩井 根立寺 嘉永2年(1849)
http://www.wasan.jp/niigata/konryuji.html
涌田和芳，外川一仁: 三島根立寺の算額，長岡工業高等専門学校研究紀要，第53巻(2017)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_53/53_17wakuta.pdf

直線上に楕円と天円，地円，人円がある。地円の直径が 1 寸のとき，極大となる天円の直径はいかほどか。

楕円の長半径，短半径，中心座標を a, b, (0, b)  注: この問題にいおいては，a < b である
天円の半径と中心座標を r1, (0, y1)
人円の半径と中心座標を r3, (0, 2r2 + r3)
地円の半径と中心座標を r2, (0, r2)
とおく。
算法助術の公式 85, 86, 86 を使用し，以下の連立方程式を解く。
天円は楕円の曲率円にあたるので，その半径は a^2/b である。
なお，楕円の左右の人円はこの問題を解く上では無関係である。

include("julia-source.txt");

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

   1-element Vector{NTuple{4, Sym}}:
    (3/2, 5/2, 3*sqrt(3)/2, 9/2)

天円の半径は 2/3 となり，地円の 3 倍である。すなわち，天円の直径は 3 寸である。
その他のパラメータは以下のとおりである。

r1 = 1.5;  r3 = 2.5;  a = 2.59808;  b = 4.5
y1 = 7.5;  y3 = 3.5

include("julia-source.txt");
function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (r1, r3, a, b) = (3/2, 5/2, 3*sqrt(3)/2, 9/2)
   y3 = 2r2 + r3
   y1 = 2r2 + 2r3 + r1
   @printf("r1 = %g;  r3 = %g;  a = %g;  b = %g\n", r1, r3, a, b)
   @printf("y1 = %g;  y3 = %g\n", y1, y3)
   plot()
   ellipse(0, b, a, b, color=:magenta)
   circle(0, y1, r1)
   circle(0, r2, r2, :blue)
   circle(0, y3, r3, :green)
   d = 2sqrt(r3^2 - r2^2)
   circle(d, r3, r3, :green)
   circle(-d, r3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, b, " b", :black, :left, :vcenter)
       point(0, y3, " 人円:r3\n (0,y3)", :green, :left, :vcenter)
       point(d, r3, "(d,r3)", :green, :center, delta=-delta/2)
       point(-d, r3, "(-d,r3)", :green, :center, delta=-delta/2)
       point(0, y1, " 天円:r1\n (0,y1)", :black, :left, :vcenter)
       point(0, r2, "   地円:r2,(0,r2)", :black, :left, :vcenter)
   end
end;

算額（その697）

算額（その697）

神壁算法 關流清水與市道香門人 上総國長柄郡上之郷 本石與八利重 寛政七年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

不等辺三角形の辺の長さが長い順に，大斜，中斜，小斜がそれぞれ 345 寸，322 寸，299 寸である。三角形内に斜線を 2 本引き，区画された領域に甲円と乙円をいれる。甲円の直径が 115 寸のとき，乙円の直径はいかほどか。

三角形の辺の長さを「大斜」，「中斜」，「小斜」，中斜と小斜の交点座標を (x0, y0)，斜線と小斜，中斜との交点座標を (xa, ya), (xb, yb) とする。
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 大斜::positive, 中斜::positive, 小斜::positive,
     x0::positive, y0::positive,
     xa::positive, xb::positive,
     r1::positive, x1::positive,
     r2::positive, x2::positive, y2::positive

ya = y0/x0 * xa
yb = y0/(大斜 - x0) * (大斜 - xb)
eq1 = x0^2 + y0^2 - 小斜^2
eq2 = (大斜 - x0)^2 + y0^2 - 中斜^2
eq3 = distance(x0, y0, 大斜, 0, x2, y2) - r2^2
eq4 = distance(0, 0, x0, y0, x2, y2) - r2^2
eq5 = distance(0, 0, xb, yb, x1, r1) - r1^2
eq6 = distance(0, 0, xb, yb, x2, y2) - r2^2
eq7 = distance(xa, ya, 大斜, 0, x1, r1) - r1^2
eq8 = distance(xa, ya, 大斜, 0, x2, y2) - r2^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x0, y0, xa, xb, x1, r2, x2, y2) = u
   return [
       x0^2 + y0^2 - 小斜^2,  # eq1
       y0^2 - 中斜^2 + (-x0 + 大斜)^2,  # eq2
       -r2^2 + (x2 - (x0^2*x2 - 2*x0*x2*大斜 + x0*y0*y2 + x2*大斜^2 + y0^2*大斜 - y0*y2*大斜)/(x0^2 - 2*x0*大斜 + y0^2 + 大斜^2))^2 + (-y0*(x0*x2 - x0*大斜 - x2*大斜 + y0*y2 + 大斜^2)/(x0^2 - 2*x0*大斜 + y0^2 + 大斜^2) + y2)^2,  # eq3
       -r2^2 + (-x0*(x0*x2 + y0*y2)/(x0^2 + y0^2) + x2)^2 + (-y0*(x0*x2 + y0*y2)/(x0^2 + y0^2) + y2)^2,  # eq4
       -r1^2 + (r1 - y0*(r1*xb^2*y0 - 2*r1*xb*y0*大斜 + r1*y0*大斜^2 + x0*x1*xb^2 - x0*x1*xb*大斜 - x1*xb^2*大斜 + x1*xb*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2))^2 + (x1 - xb*(r1*x0*xb*y0 - r1*x0*y0*大斜 - r1*xb*y0*大斜 + r1*y0*大斜^2 + x0^2*x1*xb - 2*x0*x1*xb*大斜 + x1*xb*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2))^2,  # eq5
       -r2^2 + (x2 - xb*(x0^2*x2*xb - 2*x0*x2*xb*大斜 + x0*xb*y0*y2 - x0*y0*y2*大斜 + x2*xb*大斜^2 - xb*y0*y2*大斜 + y0*y2*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2))^2 + (-y0*(x0*x2*xb^2 - x0*x2*xb*大斜 - x2*xb^2*大斜 + x2*xb*大斜^2 + xb^2*y0*y2 - 2*xb*y0*y2*大斜 + y0*y2*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2) + y2)^2,  # eq6
       -r1^2 + (r1 - xa*y0*(r1*xa*y0 + x0*x1*xa - x0*x1*大斜 - x0*xa*大斜 + x0*大斜^2)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2))^2 + (x1 - (r1*x0*xa^2*y0 - r1*x0*xa*y0*大斜 + x0^2*x1*xa^2 - 2*x0^2*x1*xa*大斜 + x0^2*x1*大斜^2 + xa^2*y0^2*大斜)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2))^2,  # eq7
       -r2^2 + (x2 - (x0^2*x2*xa^2 - 2*x0^2*x2*xa*大斜 + x0^2*x2*大斜^2 + x0*xa^2*y0*y2 - x0*xa*y0*y2*大斜 + xa^2*y0^2*大斜)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2))^2 + (-xa*y0*(x0*x2*xa - x0*x2*大斜 - x0*xa*大斜 + x0*大斜^2 + xa*y0*y2)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2) + y2)^2,  # eq8
   ]
end;

(大斜, 中斜, 小斜, r1) = (345, 322, 299, 115//2)
iniv = BigFloat[152, 258, 105, 219, 165, 42, 162, 182]
res = nls(H, ini=iniv)

   (BigFloat[151.8000000000000000000000000000000000000000000000000000000000000000000000000004, 257.6000000000000000000000000000000000000000000000000000000000000000000000000009, 100.0, 213.75, 160.9999999999999999999999999999999999999999999999999999999999999999999999999956, 42.00000000000000000000000000000000000000000000000000000000000000000000000000111, 156.0, 182.0], true)

乙円の直径は 84 寸である。

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

x0 = 151.8;  y0 = 257.6;  xa = 100;  xb = 213.75;  x1 = 161;  r2 = 42;  x2 = 156;  y2 = 182

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (大斜, 中斜, 小斜, r1) = (345, 322, 299, 115//2)
   (x0, y0, xa, xb, x1, r2, x2, y2) = res[1]
   println("乙円の直径 = $(Float64(2r2))")
   @printf("x0 = %g;  y0 = %g;  xa = %g;  xb = %g;  x1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n",
           x0, y0, xa, xb, x1, r2, x2, y2)
   ya = y0/x0 * xa
   yb = y0/(大斜 - x0) * (大斜 - xb)
   plot([0, 大斜, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, y2, r2, :blue)
   segment(0, 0, xb, yb, :green)
   segment(xa, ya, 大斜, 0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(大斜, 0, "大斜", :black, :left, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :black, :left, :bottom, delta=delta/2)
       point(xa, ya, "(xa,ya)", :black, :right, :bottom, delta=delta/2)
       point(xb, yb, "(xb,yb)", :black, :left, :bottom, delta=delta/2)
       point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, :top, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, :top, delta=-delta/2)
   end
end;

算額（その696）

算額（その696）

神壁算法 關流神谷幸吉定令門人 本田帶刀家士 佐久間常右衛門久豊
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

半円の中に大円および甲円から始まる累円が左右に位置している。大円の直径が与えられ，末円の径がわかったときに累円の個数はいくつあるか。

まず，累円の開始点として，甲円の半径と中心座標を求める。
大円の半径と中心座標を r0, (0, r0)
甲円の半径と中心座標を r1, (x1, y1); y1 = r1

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, x1::positive

eq1 = x1^2 + r1^2 - (2r0 - r1)^2
eq2 = x1^2 + (r0 - r1)^2 - (r0 + r1)^2
res2 = solve([eq1, eq2], (r1, x1))

   1-element Vector{Tuple{Sym, Sym}}:
    (r0/2, sqrt(2)*r0)

甲円の半径は大円の半径を r0 とすると r0/2，中心座標は （√2r0, r0/2) である。
甲円に外接する乙円のパラメータ（半径，中心座標）は，以下のようになる。

@syms y1::positive, r2::positive, x2::positive, y2::positive

eq3 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq4 = x2^2 + (r0 - y2)^2 - (r0 + r2)^2
eq5 = x2^2 + y2^2 - (2r0 - r2)^2
res3 = solve([eq3, eq4, eq5], (r2, x2, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    ((8*r0^3 + 4*r0^2*r1 - 20*r0^2*y1 - 2*r0*r1^2 - 4*r0*r1*y1 + 10*r0*x1^2 + 14*r0*y1^2 - r1^3 + 3*r1^2*y1 + r1*x1^2 + r1*y1^2 - 3*x1^2*y1 - 2*sqrt(2)*x1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 3*y1^3)/(2*(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2)), (8*r0^2*x1 - 4*r0*r1*x1 - 4*r0*x1*y1 + 2*sqrt(2)*r0*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 4*r1^2*x1 + sqrt(2)*r1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) + 4*x1^3 + 4*x1*y1^2 - 3*sqrt(2)*y1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4))/(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2), 3*sqrt(2)*x1*sqrt(-(-4*r0^2 - 4*r0*r1 - r1^2 + x1^2 + y1^2)*(2*r0*r1 - 2*r0*y1 - r1^2 + x1^2 + y1^2))/(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2) + (-8*r0^3 + 4*r0^2*r1 + 12*r0^2*y1 + 10*r0*r1^2 - 12*r0*r1*y1 + 2*r0*x1^2 - 6*r0*y1^2 + 3*r1^3 - 9*r1^2*y1 - 3*r1*x1^2 - 3*r1*y1^2 + 9*x1^2*y1 + 9*y1^3)/(2*(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2)))
    ((8*r0^3 + 4*r0^2*r1 - 20*r0^2*y1 - 2*r0*r1^2 - 4*r0*r1*y1 + 10*r0*x1^2 + 14*r0*y1^2 - r1^3 + 3*r1^2*y1 + r1*x1^2 + r1*y1^2 - 3*x1^2*y1 + 2*sqrt(2)*x1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 3*y1^3)/(2*(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2)), (8*r0^2*x1 - 4*r0*r1*x1 - 4*r0*x1*y1 - 2*sqrt(2)*r0*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 4*r1^2*x1 - sqrt(2)*r1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) + 4*x1^3 + 4*x1*y1^2 + 3*sqrt(2)*y1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4))/(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2), -3*sqrt(2)*x1*sqrt(-(-4*r0^2 - 4*r0*r1 - r1^2 + x1^2 + y1^2)*(2*r0*r1 - 2*r0*y1 - r1^2 + x1^2 + y1^2))/(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2) + (-8*r0^3 + 4*r0^2*r1 + 12*r0^2*y1 + 10*r0*r1^2 - 12*r0*r1*y1 + 2*r0*x1^2 - 6*r0*y1^2 + 3*r1^3 - 9*r1^2*y1 - 3*r1*x1^2 - 3*r1*y1^2 + 9*x1^2*y1 + 9*y1^3)/(2*(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2)))

2 組の解が得られるが，最初のものが適解である。
乙円の半径はかなり長い式になり，これ以上簡約化もできないが，この式は，「一つ前の累円のパラメータから，それに外接する累円のパラメータを得る」漸化式である。

これを関数にすると，以下のようになる。
現在の累円の半径，中心座標を与えると，次の累円の半径，中心座標を返す。

nextcircle(r0, r1, x1, y1) = ((8*r0^3 + 4*r0^2*r1 - 20*r0^2*y1 - 2*r0*r1^2 - 4*r0*r1*y1 + 10*r0*x1^2 + 14*r0*y1^2 - r1^3 + 3*r1^2*y1 + r1*x1^2 + r1*y1^2 - 3*x1^2*y1 - 2*sqrt(2)*x1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 3*y1^3)/(2*(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2)), (8*r0^2*x1 - 4*r0*r1*x1 - 4*r0*x1*y1 + 2*sqrt(2)*r0*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 4*r1^2*x1 + sqrt(2)*r1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) + 4*x1^3 + 4*x1*y1^2 - 3*sqrt(2)*y1*sqrt(8*r0^3*r1 - 8*r0^3*y1 + 4*r0^2*r1^2 - 8*r0^2*r1*y1 + 4*r0^2*x1^2 + 4*r0^2*y1^2 - 2*r0*r1^3 - 2*r0*r1^2*y1 + 2*r0*r1*x1^2 + 2*r0*r1*y1^2 + 2*r0*x1^2*y1 + 2*r0*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4))/(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2), 3*sqrt(2)*x1*sqrt(-(-4*r0^2 - 4*r0*r1 - r1^2 + x1^2 + y1^2)*(2*r0*r1 - 2*r0*y1 - r1^2 + x1^2 + y1^2))/(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2) + (-8*r0^3 + 4*r0^2*r1 + 12*r0^2*y1 + 10*r0*r1^2 - 12*r0*r1*y1 + 2*r0*x1^2 - 6*r0*y1^2 + 3*r1^3 - 9*r1^2*y1 - 3*r1*x1^2 - 3*r1*y1^2 + 9*x1^2*y1 + 9*y1^3)/(2*(4*r0^2 + 4*r0*r1 - 12*r0*y1 + r1^2 - 6*r1*y1 + 8*x1^2 + 9*y1^2)));

nextcircle(0.5, 0.25, 0.7071067811865476, 0.25)

   (0.12773958089728293, 0.6167812573081513, 0.6167812573081511)

nextcircle(0.5, 0.12773958089728296, 0.6167812573081511, 0.6167812573081513)

   (0.07322330470336305, 0.4999999999999999, 0.7803300858899107)

この関数を用いて，甲円を初期値として，乙円，丙円のパラメータを求める。

各行の後半 4 項目，r0/rn, f(n), n については後に説明する。

r0 = 1/2
r, x, y = nextcircle(r0, r0/2, sqrt(2)*r0, r0/2)  # 2
for i = 2:10
   println(i, ",", r, ",", x, ",", y, ", r0/rn = ", r0/r, ", f(n) = ", 0.5*i^2 + 0.414213562374*i + 1.08578643763, ", n = ", round(Int, sqrt(2*r0/r - 2) - sqrt(2) + 1))
   r, x, y = nextcircle(r0, r, x, y)
end

   2,0.12773958089728293,0.6167812573081513,0.6167812573081511, r0/rn = 3.914213562373095, f(n) = 3.914213562378, n = 2
   3,0.07322330470336313,0.5,0.7803300858899105, r0/rn = 6.828427124746189, f(n) = 6.828427124752, n = 3
   4,0.04654349098723124,0.41090581831205225,0.8603695270383063, r0/rn = 10.742640687119284, f(n) = 10.742640687126, n = 4
   5,0.03193489522432071,0.34580468567296235,0.9041953143270376, r0/rn = 15.65685424949239, f(n) = 15.6568542495, n = 5
   6,0.023179195594803567,0.2973526214981749,0.9304624132155892, r0/rn = 21.571067811865422, f(n) = 21.571067811874, n = 6
   7,0.017552924734392517,0.26028226525009357,0.9473412257968229, r0/rn = 28.48528137423842, f(n) = 28.485281374247997, n = 7
   8,0.013736454334619978,0.23116292072255745,0.95879063699614, r0/rn = 36.399494936611866, f(n) = 36.399494936622, n = 8
   9,0.011034188473256403,0.20775641354942292,0.9668974345802309, r0/rn = 45.31370849898491, f(n) = 45.313708498996, n = 9
   10,0.009053391497230267,0.18856790503185952,0.9728398255083086, r0/rn = 55.22792206135862, f(n) = 55.22792206137, n = 10

甲円を1番目として，n 番目の累円の半径を rn として，r0/rn を n の多項式で表すことを考える。
多項式回帰分析により，r0/rn = f(n) = 0.5*n^2 + 0.414213562374*n + 1.08578643763 が得られる。
この多項式の逆関数を得れば，r0/rn から n を求めることができる。
逆関数は以下のようになる。

@syms n, r0_by_rn
eq = 0.5*n^2 + 0.414213562374*n + 1.08578643763 - r0_by_rn
solve(eq, n)[2] |> println

   1.41421356237502*sqrt(0.99999999999728*r0_by_rn - 1) - 0.414213562374

1.41421356237502 は √2, 0.99999999999728 は 1, 0.414213562374 は √2 - 1 なので，以下のように書ける。

sqrt(Sym(2))*sqrt(r0_by_rn - 1) + 1 - sqrt(Sym(2))|>  simplify |> println

   sqrt(2*r0_by_rn - 2) - sqrt(2) + 1

何番目の累円かを知るには，「大円の半径を n 番目の累円の半径で割ったものを 2 倍して，2 を引いて，平方根を求め，1 を加えて，2 の平方根を引く」

これを関数にすると以下のようになる。累円の半径を与えるとそれが何番目のものかを返す。

nth(rn) = round(Int, sqrt(2r0/rn - 2) + 1 - sqrt(2));

たとえば，10 番目の累円の半径は上の実行結果から,0.009053391497230267 なので，何番目かは nth(0.009053391497230267) = 10 になる。

nth(0.009053391497230267)

   10

function draw(more=false)
   pyplot(size=(800, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   names = Char["甲乙丙丁戊己庚辛壬癸"...]
   s2 = sqrt(2)
   plot()
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   r0 = 1/2
   circlef(0, 0, 2r0, :cyan, beginangle=0, endangle=180)
   segment(-2r0, 0, 2r0, 0)
   circlef(0, r0, r0, 1)
   (r1, x1, y1) = (r0/2, sqrt(2)*r0, r0/2)
   for i in 2:25
       circlef(x1, y1, r1, i)
       circlef(-x1, y1, r1, i)
       i < 8 && point(x1, y1, names[i - 1], :black, :center, :vcenter, mark=false)
       (r1, x1, y1) = nextcircle(r0, r1, x1, y1)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0, " 大円:r0,(0,r0)", :black, :left, :vcenter)
   end
end;