算額（その985）

算額（その985）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

楕円の中に2本の斜線を引き，区分された領域に天円 2 個，地円 1 個，人円 2 個を入れる。楕円の長径，短径がそれぞれ 69 寸，23 寸のとき，人円の直径はいかほどか。

楕円の長半径，短半径を a, b
天円の半径と中心座標を r1, (r2 + r1, 0)
地円の半径と中心座標を r2, (0, 0)
人円の半径と中心座標を r3, (r2 - r3, 0)
原点を通る，地円と人円の共通接線と楕円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

1.　天円の半径

天円の半径は，「算法助術」の公式84により求めることができる。
楕円に内接する同じ大きさの 2 円の半径 r は，楕円の長半径 a，短半径 b，原点から円の中心までの距離 dである。
d^2 = (a^2 - b^2)*(b^2 - r^2)/b^2

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     x0::positive, y0::positive,
     r1::positive, r2::positive, r3::positive,
     x01::positive, y01::positive, d
r2 = b
d = r2 + r1
eq1 = d^2 - (a^2 - b^2)*(b^2 - r1^2)/b^2
r1 = solve(eq1, r1)[1]
r1 |> println

   b - 2*b^3/a^2

r1(a => 69/2, b => 23/2).evalf() |> println

   8.94444444444444

r1 は a と b の関数で，b - 2*b^3/a^2 と表すことができる。
a = 69/2, b = 23/2 を代入すると 8.94444444444444 と正しい答えになっている。

2.　人円の半径

次いで，r1 が確定したとして，人円 r3 を求める。これは，人円と天円の中心から斜線までの距離と，原点からそれぞれの中心までの距離の関係から導かれる 3元連立方程式を解くことで求めることができる。

eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = r3/(r2 - r3) - y0/sqrt(x0^2 + y0^2)
eq3 = r1/(r2 + r1) - y0/sqrt(x0^2 + y0^2)
res = solve([eq1, eq2, eq3], (x0, y0, r3))[4];  # 4 of 4

println("x0: ", res[1])
println("y0: ", res[2])
println("r3: ", res[3])

   x0: sqrt(a*sqrt(a + b) - (a^2 - 2*b^2)*sqrt(1/(a - b)))*sqrt(a*sqrt(a + b) + (a^2 - 2*b^2)*sqrt(1/(a - b)))/sqrt(a + b)
   y0: b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a*sqrt(a + b))
   r3: b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a^2*sqrt(1/(a - b)) - 2*b^2*sqrt(1/(a - b)) + 2*sqrt((a^4 - 2*a^2*b^2 + b^4)/(a^2 - b^2))*sqrt(a + b))

res[3](a => 69/2, b => 23/2).evalf() |> println

   3.50000000000000

楕円の長径，短径が それぞれ 69 寸，23 寸のとき，人円の直径は 7 寸である。

3.　連立方程式を一挙に解くには

なお，SymPy の能力的には無理であるが，eq1〜eq6 を連立方程式として解けば，解は一挙に求まる。
x01, y01 は天円と楕円の接点座標である。

eq4 = x01^2/a^2 + y01^2/b^2 - 1
eq5 = (x01 - d)^2 + y01^2 - r1^2
eq6 = -b^2*x01/(a^2*y01) +(x01 - d)/y01;

--- 結果の検証のための描画プログラム

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r2) = (69, 23) ./ 2
   b = r2
   r1 = b - 2*b^3/a^2
   x0 = sqrt(a*sqrt(a + b) - (a^2 - 2*b^2)*sqrt(1/(a - b)))*sqrt(a*sqrt(a + b) + (a^2 - 2*b^2)*sqrt(1/(a - b)))/sqrt(a + b)
   y0 = b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a*sqrt(a + b))
   r3 = b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a^2*sqrt(1/(a - b)) - 2*b^2*sqrt(1/(a - b)) + 2*sqrt((a^4 - 2*a^2*b^2 + b^4)/(a^2 - b^2))*sqrt(a + b))
   @printf("長径が %g 寸，短径が %g 寸のとき，人円の直径は %g 寸である。\n", 2a, 2b, 2r3)
   plot()
   ellipse(0, 0, a, b, color=:green)
   circle2(r2 + r1, 0, r1)
   circle2(r2 - r3, 0, r3, :magenta)
   circle(0, 0, r2, :blue)
   segment(-x0, -y0, x0, y0, :orange)
   segment(-x0, y0, x0, -y0, :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(r2 + r1, 0, "天円:r1\n(r2+r1,r1)", :red, :center, delta=-2delta)
       point(0, 3b/4, "地円:r2\n(0,0)", :blue, :center, delta=-2delta, mark=false)
       point(r2 - r3, 0, "人円:r3\n(r2-r3,r3)", :black, :center, delta=-2delta, deltax=-4delta)
       point(a, 0, " a", :green, :left, :vcenter)
       point(0, b, "b = r2", :blue, :center, :bottom, delta=2delta)
   end
end;

算額（その984）

算額（その984）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に，甲円 1 個，乙円 2 個，丙円 2 個が入っている。外円，甲円，乙円の直径がそれぞれ 156 寸，52 寸，26 寸のとき，丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, )
甲円の半径と中心座標を r1, (0, y1)
乙円の半径と中心座標を r2, (0, y2)
丙円の半径と中心座標を r3, (x3, y3)
とおいて，以下の連立方程式を解く。
注：丙円の中心は x 軸上にあるわけではない。y3 ≠ 0

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, y1::negative,
     r2::positive, y2::positive,
     r3::positive, x3::positive, y3::negative
eq1 = x3^2 + (y1 - y3)^2 - (r3 - r1)^2
eq2 = x3^2 + (y2 - y3)^2 - (r3 - r2)^2
eq3 = x3^2 + (y3 - r2 + R)^2 - (r2 + r3)^2
eq4 = x3^2 + y3^2 - (R - r3)^2
eq5 = (y2 - y1) - (r1 + r2);
res = solve([eq1, eq2, eq3, eq4, eq5], (r3, y1, y2, x3, y3))[1];

println("r3: ", res[1])
println("y1: ", res[2])
println("y2: ", res[3])
println("x3: ", res[4])
println("y3: ", res[5])

   r3: (-R + r2)*(r1 + r2)*(-2*R*r1 - R*r2 + r1*r2)/((R + r1)*(R*r1 + r2^2))
   y1: -(-R + r1)*(-R + r1 + 2*r2)/(R + r1)
   y2: -(R^2 - 3*R*r1 - 3*R*r2 + r1*r2)/(R + r1)
   x3: 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-(-R + r2)*(-R + r1 + 2*r2)*(-2*R*r1 - R*r2 + r1*r2))*sqrt(r1 + r2)/((R + r1)*(R*r1 + r2^2))
   y3: (-R^3*r1 + R^2*r1^2 + 3*R^2*r1*r2 + R*r1^2*r2 + R*r1*r2^2 + R*r2^3 - r1^2*r2^2 - r1*r2^3)/(R^2*r1 + R*r1^2 + R*r2^2 + r1*r2^2)

外円，甲円，乙円の直径がそれぞれ 156 寸, 52 寸, 26 寸のとき，丙円の直径は 105 寸である。

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

r1 = 26;  r2 = 13;  R = 78
r3 = 52.5;  y1 = -13;  y2 = 26;  x3 = 25.0998;  y3 = -4.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2) = (156, 52, 26) ./ 2
   (r3, y1, y2, x3, y3) =  ((-R + r2)*(r1 + r2)*(-2*R*r1 - R*r2 + r1*r2)/((R + r1)*(R*r1 + r2^2)), -(-R + r1)*(-R + r1 + 2*r2)/(R + r1), -(R^2 - 3*R*r1 - 3*R*r2 + r1*r2)/(R + r1), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-(-R + r2)*(-R + r1 + 2*r2)*(-2*R*r1 - R*r2 + r1*r2))*sqrt(r1 + r2)/((R + r1)*(R*r1 + r2^2)), (-R^3*r1 + R^2*r1^2 + 3*R^2*r1*r2 + R*r1^2*r2 + R*r1*r2^2 + R*r2^3 - r1^2*r2^2 - r1*r2^3)/(R^2*r1 + R*r1^2 + R*r2^2 + r1*r2^2))
   @printf("外円，甲円，乙円の直径がそれぞれ %g, %g, %g のとき，丙円の直径は %g である。\n", 2R, 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  R = %g\n", r1, r2, R)
   @printf("r3 = %g;  y1 = %g;  y2 = %g;  x3 = %g;  y3 = %g\n", r3, y1, y2, x3, y3)
   plot()
   circle(0, 0, R, :green)
   circle(0, y1, r1)
   circle(0, y2, r2, :blue)
   circle(0, r2 - R, r2, :blue)
   circle2(x3, y3, r3, :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(0, y1, "甲円:r1\n(0,y1)", :brown, :center, delta=-delta/2)
       point(0, y2, "乙円:r2\n(0,y2)", :blue, :center, delta=-delta/2)
       point(0, r2 - R, "乙円:r2\n(0,r2-R)", :blue, :center, :bottom, delta=delta/2)
       point(x3, y3, " 丙円:r3,(x3,y3)", :blue, :left, :vcenter)
       point(R, 0, " R", :green, :left, :bottom, delta=delta/2)
   end
end;

算額（その983）

算額（その983）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

半円の中に甲円 1 個，乙円 4 個が入っている。甲円の直径が 94.9 寸のとき，乙円の直径はいかほどか。

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

include("julia-source.txt");

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

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

   r1*(1 + sqrt(10))/5
   r1*(2*sqrt(10) + 11)/5

乙円の半径は甲円の半径の (1 + √10)/5 倍である。
甲円の直径が 94.9 寸のとき，乙円の直径は 94.9(1 + √10)/5 = 79.00002998999585 寸である。

なお，半円（外円）の直径は甲円の直径の (2√10 + 11)/5 倍の 328.8200599799917 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 94.9/2
   r2 = r1*(1 + √10)/5
   R = r1*(2√10 + 11)/5
   @printf("甲円の直径が %g のとき，乙円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  R = %g\n", r1, r2, R)
   plot()
   circle(0, 0, R, :green, beginangle=0, endangle=180)
   segment(-R, 0, R, 0, :green)
   circle(0, R - r1, r1)
   circle2(r2, r2, r2, :blue)
   circle2(3r2, 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, R, "半円:R,(0,0)", :green, :center, :bottom, delta=delta)
       point(0, R - r1, "甲円:r1\n(0,R-r1)", :brown, :center, delta=-delta/2)
       point(r2, r2, "小円:r2\n(r2,r2)", :blue, :center, delta=-delta/2)
       point(3r2, r2, "小円:r2\n(3r2,r2)", :blue, :center, delta=-delta/2)
   end
end;

算額（その982）

算額（その982）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

直線上に大円，中円，小円が 2 個ずつ積み上がっている。大円の直径が 58 寸のとき，中円の直径はいかほどか。

大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (r2, y2)
小円の半径と中心座標を r3, (0, r3), (0, 3r3); r3 = r1/2
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive,
     r2::positive, y2::positive,
     r3::positive
r3 = r1/2
eq1 = (x1 - r2)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq2 = x1^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = r2^2 + (y2 - 3r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (x1, r2, y2))[1]  # 1 of 2

   (sqrt(2)*r1, r1*(-6 - 2*sqrt(1 + 5*sqrt(2)) + 4*sqrt(2 + 10*sqrt(2)) + 19*sqrt(2))/49, r1*(sqrt(2)/7 + 2*sqrt(1 + 5*sqrt(2))/7 + 10/7))

中円の半径は大円の半径の (4sqrt(2 + 10√2) + 19√2 - 6 - 2sqrt(1 + 5√2))/49 倍である。
大円の直径が 58 寸のとき，中円の直径は 58*(4sqrt(2 + 10√2) + 19√2 - 6 - 2sqrt(1 + 5√2))/49 = 37.00050035355881 寸である。

その他のパラメータは以下のとおりである。
r1 = 29;  x1 = 41.0122;  r2 = 18.5003;  y2 = 70.8269

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 58/2
   r3 = r1/2
   (x1, r2, y2) = (
       √2r1,
       r1*(4sqrt(2 + 10√2) + 19√2 - 6 - 2sqrt(1 + 5√2))/49,
       r1*(√2/7 + 2sqrt(1 + 5√2)/7 + 10/7))
   @printf("大円の直径が %g のとき，中円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  y2 = %g\n", r1, x1, r2, y2)
   plot()
   circle2(x1, r1, r1)
   circle2(r2, y2, r2, :blue)
   circle(0, r3, r3, :green)
   circle(0, 3r3, 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(x1, r1, "大円:r1,(x1,r1)", :green, :center, delta=-delta/2)
       point(r2, y2, "中円:r2\n(r2,y2)", :brown, :center, delta=-delta/2)
       point(0, r3, "小円:r3\n(0,r3)", :magenta, :center, delta=-delta/2)
       point(0, 3r3, "小円:r3\n(0,3r3)", :magenta, :center, delta=-delta/2)
       segment(-x1 - r1, 0, x1 + r1, 0, :black, lw=1)
   end
end;

算額（その981）

算額（その981）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

正方形内に四分円 3 個，甲円 1 個，乙円 1 個，丙円 2 個を入れる。乙円の直径が 3.1 寸（注）のとき，丙円の直径はいかほどか。

注：後述の通り，2.1 寸は 3.1 寸の誤り。

正方形の一辺の長さを 2a
甲円の半径と中心座標を r1, (a, a)
乙円の半径と中心座標を r2, (x2, 2a - x2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

正方形，甲円，乙円を決定する方程式 eq3, eq4, eq5 は独立なので，まず a, r1, x2 を求める。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive, y3::positive
# eq1 = (a - x3)^2 + (y3 - a)^2 - (r1 + r3)^2 |> expand
# eq2 = (2a - x3)^2 + (2a - y3)^2 - (2a - r3)^2 |> expand
eq3 = 2a^2 - (2a - r1)^2 |> expand
eq4 = (2a - x2)^2 + x2^2 - (2a - r2)^2 |> expand
eq5 = 2(x2 - a)^2 - (r1 + r2)^2 |> expand;
# eq6 = dist2(2a, 0, 0, 2a, x3, y3, r3)
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (a, r1, x2, r3, x3, y3))

println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   -2*a^2 + 4*a*r1 - r1^2,  # eq3
   4*a*r2 - 4*a*x2 - r2^2 + 2*x2^2,  # eq4
   2*a^2 - 4*a*x2 - r1^2 - 2*r1*r2 - r2^2 + 2*x2^2,  # eq5

res = solve([eq3, eq4, eq5], (a, r1, x2))[2]  # 2 of 2

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

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

   4.83804653251745
   2.83406204335659
   7.93804653251745

a, r1, x2 が既知として，eq1, eq2, eq6 から 最終的に r3 を求める（r2 を含む式）。

include("julia-source.txt");
using SymPy
@syms a::positive, r1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive, y3::positive
a = r2*(2 + 3*sqrt(Sym(2)))/2
r1 = r2*(-1 + 2*sqrt(Sym(2)))
x2 = r2*(3*sqrt(Sym(2))/2 + 3)
eq1 = (a - x3)^2 + (y3 - a)^2 - (r1 + r3)^2 |> expand
eq2 = (2a - x3)^2 + (2a - y3)^2 - (2a - r3)^2 |> expand
eq6 = dist2(2a, 0, 0, 2a, x3, y3, r3);

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq6, ",  # eq6")

   2*r2^2 + 10*sqrt(2)*r2^2 - 4*sqrt(2)*r2*r3 + 2*r2*r3 - 3*sqrt(2)*r2*x3 - 2*r2*x3 - 3*sqrt(2)*r2*y3 - 2*r2*y3 - r3^2 + x3^2 + y3^2,  # eq1
   12*sqrt(2)*r2^2 + 22*r2^2 + 4*r2*r3 + 6*sqrt(2)*r2*r3 - 6*sqrt(2)*r2*x3 - 4*r2*x3 - 6*sqrt(2)*r2*y3 - 4*r2*y3 - r3^2 + x3^2 + y3^2,  # eq2
   6*sqrt(2)*r2^2 + 11*r2^2 - 3*sqrt(2)*r2*x3 - 2*r2*x3 - 3*sqrt(2)*r2*y3 - 2*r2*y3 - r3^2 + x3^2/2 + x3*y3 + y3^2/2,  # eq6

eq1, eq2 から x3, y3 を求める。

res = solve([eq1, eq2], (x3, y3))[1]  # 1 of 2

   (-r2 + 2*sqrt(2)*r2 - 2*sqrt(r3)*sqrt(6733*sqrt(2)*r2 + 9608*r2 - 2194*r3 - 1513*sqrt(2)*r3)/(58 + 45*sqrt(2)) - sqrt(2)*r3/2 + 2*r3, -r2 + 2*sqrt(2)*r2 + 2*sqrt(r3)*sqrt(6733*sqrt(2)*r2 + 9608*r2 - 2194*r3 - 1513*sqrt(2)*r3)/(58 + 45*sqrt(2)) - sqrt(2)*r3/2 + 2*r3)

res[1](r2=>3.1/2, r3 => 1.0019922445804257).evalf() |> println
res[2](r2=>3.1/2, r3 => 1.0019922445804257).evalf() |> println

   1.51119910595402
   6.74786293740257

eq6 に x3, y3 を代入する。

eq6 = eq6(x3 => res[1], y3 => res[2]) |> simplify;

ans_r3 = solve(eq6, r3)[2]
ans_r3 |> println

   -r2*sqrt(9 - 4*sqrt(2))/4 - sqrt(2)*r2*sqrt(9 - 4*sqrt(2))/8 + sqrt(2)*r2/8 + 5*r2/4

簡約化する。

ans_r3 = ans_r3 |> sympy.sqrtdenest |> simplify
ans_r3 |> println

   r2*(4 - sqrt(2))/4

正しい答えが得られることを確認する。

ans_r3(r2 => 3.1/2).evalf() |> println

   1.00199224458043

丙円の半径は，乙円の半径の (4 - sqrt(2))/4 倍である。
乙円の直径が 3.1 寸のとき，丙円の直径は 3.1 * (4 - sqrt(2))/4 = 2.0039844891608514 である。

逆にたどって，全てのパラメータの解を求める。

r2 = 3.1/2
r3 = r2*(4 - √2)/4
x3 = -r2 + 2√2*r2 - 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
y3 = -r2 + 2√2*r2 + 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
a = r2*(2 + 3√2)/2
r1 = r2*(2√2 - 1)
x2 = r2*(3√2/2 + 3)

   r2 = 1.55; a = 4.83805;  r1 = 2.83406;  x2 = 7.93805;  r3 = 1.00199;  x3 = 1.5112;  y3 = 6.74786

丙円の半径は 0.6787689398770625（直径は 1.357537879754125）
「答」は「丙円径は2寸00有奇」となっている。
「術」は「(1-sqrt(0.5)/2)*乙円径」と言っているのに。計算すると 1.3575378797541253 になる。
乙円の径が 3.1 寸のとき, 丙円の径が 2.0039844891608514 になる。
結局，「答」，「術」は正しく，「問」の「只云乙円径二寸一分」が「只云乙円径三寸一分」の誤記ということだろう。算額が間違えているのか，出版時の校正漏れなのかはわからない。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 3.1/2
   r3 = r2*(4 - √2)/4
   x3 = -r2 + 2√2*r2 - 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
   y3 = -r2 + 2√2*r2 + 2sqrt(r3)*sqrt(6733√2*r2 + 9608r2 - 2194r3 - 1513√2*r3)/(58 + 45√2) - √2r3/2 + 2r3
   a = r2*(2 + 3√2)/2
   r1 = r2*(2√2 - 1)
   x2 = r2*(3√2/2 + 3)
   @printf("乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r2, 2r3)
   @printf("r2 = %g; a = %g;  r1 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r2, a, r1, x2, r3, x3, y3)
   plot([0, 2a, 2a, 0, 0], [0, 0, 2a, 2a, 0], color=:blue, lw=0.5)
   circle(0, 0, 2a, beginangle=0, endangle=90)
   circle(2a, 2a, 2a, beginangle=180, endangle=270)
   circle(0, 2a, 2a, beginangle=270, endangle=360)
   circle(a, a, r1, :green)
   circle(x2, 2a - x2, r2, :brown)
   circle(x3, y3, r3, :magenta)
   circle(2a - y3, 2a - x3, r3, :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(a, a, "甲円:r1,(a,a)", :green, :center, delta=-delta/2)
       point(x2, 2a - x2, "乙円:r2,(x2,2a-x2)", :brown, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3\n(x3,y3)", :magenta, :center, delta=-delta/2)
       point(2a, 2a, "(2a,2a)", :blue, :right, :bottom, delta=delta/2)
   end
end;