算額（その29）

算額（その29）

京都府東山区 安井金比羅宮 平成 7 年
http://www.wasan.jp/kyoto/yasuikonpira3.html

(2) 大阪府豊中市服部元町 服部天神社 天保14年(1843)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

図のように，円の中に大円が 5 個，小円が 1 個ある。それぞれの円の径を求めよ。

外円のないものは「算額（その682）」

服部天神社のものは「答」の誤差が大きすぎる。

外円の半径を 1 とし，以下のように記号を定め，方程式を立てる。

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive, k, m;

以下の 4 方程式を立て，解く。

eq1 = x2 - (r1 + r2)cos(PI/10);
eq2 = y2 - (r1 + r2)sin(PI/10);
eq3 = r1 + 2r2 - 1;
eq4 = x2^2 + (r1 + r2 - y2)^2 - 4r2^2;

res = solve([eq1, eq2, eq3, eq4], (r1, r2, x2, y2));

name = ["r1", "r2", "x2", "y2"]
for i in 1:4
   println("$(name[i]) = $(simplify(res[1][i])) = $(res[1][i].evalf())")
end

   r1 = -4*sqrt(5) - 3*sqrt(10 - 2*sqrt(5)) + sqrt(50 - 10*sqrt(5)) + 11 = 0.259616183682500
   r2 = -5 - sqrt(50 - 10*sqrt(5))/2 + 3*sqrt(10 - 2*sqrt(5))/2 + 2*sqrt(5) = 0.370191908158750
   x2 = -sqrt(10*sqrt(5) + 50)/2 - 3*sqrt(5)/2 + 5/2 + 3*sqrt(2*sqrt(5) + 10)/2 = 0.598983089761037
   y2 = -4 - sqrt(5)*sqrt(10 - 2*sqrt(5))/2 + sqrt(10 - 2*sqrt(5)) + 2*sqrt(5) = 0.194621403573804

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = -4*sqrt(5) - 3*sqrt(10 - 2*sqrt(5)) + sqrt(50 - 10*sqrt(5)) + 11
   r2 = -5 - sqrt(50 - 10*sqrt(5))/2 + 3*sqrt(10 - 2*sqrt(5))/2 + 2*sqrt(5)
   x2 = -sqrt(10*sqrt(5) + 50)/2 - 3*sqrt(5)/2 + 5/2 + 3*sqrt(2*sqrt(5) + 10)/2
   y2 = -4 - sqrt(5)*sqrt(10 - 2*sqrt(5))/2 + sqrt(10 - 2*sqrt(5)) + 2*sqrt(5)
   println("r1 = $r1;  r2 = $r2")
   println("x2 = $x2;  y2 = $y2")
   plot()
   r0 = r1 + r2
   circle(0, 0, r1 + 2r2, :black)
   circle(0, 0, r1, :green)
   circle(0, r0, r2)
   circle(r0*cos(pi/10), r0*sin(pi/10), r2)
   circle(-r0*cos(pi/10), r0*sin(pi/10), r2)
   circle(r0*cos(17pi/10), r0*sin(17pi/10), r2)
   circle(-r0*cos(17pi/10), r0*sin(17pi/10), r2)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r1 + r2, "r1+r2 ", :red, :right)
       point(0, r1, "r1 ", :green, :right, :bottom)
       point(x2, y2, "(x2,y2)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

算額（その28）

算額（その28）

(5) 京都府宮津市天橋立文殊 知恩寺文殊堂 文政元年(1818)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

京都府宮津市天橋立 知恩寺文殊堂 文政元年
http://www.wasan.jp/kyoto/tionji1.html

キーワード：円7個

図のように，4 個の大円と，3 個の小円がある。それぞれの径を求めよ。

見ただけでわかる。大円の径は小円の径の 2 倍である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = 2, 1
   println("r1 = $r1;  r2 = $r2")
   plot()
   circle(0, 0, r1)
   circle(0, r1, r1)
   circle( r1*cos(pi/6), -r1*sin(pi/6), r1)
   circle(-r1*cos(pi/6), -r1*sin(pi/6), r1)
   circle(0, r1 + r2, r2, :blue)
   circle((r1 + r2)*cos(pi/6), -(r1 + r2)*sin(pi/6), r2, :blue)
   circle(-(r1 + r2)*cos(pi/6), -(r1 + r2)*sin(pi/6), r2, :blue)
end;

算額（その27）

算額（その27）

一〇八 北埼玉郡騎西町騎西 玉敷神社 大正4年(1915)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

東京都中央区 福徳神社 令和 2 年
http://www.wasan.jp/tokyo/hukutoku.html

埼玉県加須市騎西 玉敷神社 大正4年(1915)
http://www.wasan.jp/saitama/tamasiki.html

埼玉県北本市本宿 天神社 明治24年
山口正義：やまぶき2，第41号
https://yamabukiwasan.sakura.ne.jp/ymbk41.pdf

キーワード：円7個，外円

図のように，円の中に交わる甲円が 3 個，乙円が 3 個ある。甲円の直径が 5 のとき乙円の直径はいくつか。

以下のように記号を定め，方程式を立てる。

甲円の半径を r1，乙円の半径と中心座標を r2, (x2, y2) とする。

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive;

x2, y2 は以下のようになる。

x2 = (2r1 - r2)*sqrt(Sym(3))/2
y2 = (2r1 - r2)/2

また，甲円と乙円が接することから

eq1 = x2^2 + (r1 - y2)^2 - (r1 + r2)^2
eq1 |> simplify |> println

   r1*(2*r1 - 5*r2)

eq1 を r2 について解く。

res = solve(eq1, r2)

   r2 = 2r1/5

甲円の直径が 5 ならば乙円の直径は 2 である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = 5, 2
   println("r1 = $r1;  r2 = $r2")
   x2, y2 = (sqrt(3)*(10 - r2)/2, 5 - r2/2)
   plot()
   circle(0, 0, 2r1)
   circle(0, r1, r1, :green)
   circle(r1*cos(pi/6), -r1*sin(pi/6), r1, :green)
   circle(-r1*cos(pi/6), -r1*sin(pi/6), r1, :green)
   circle(x2, y2, r2, :magenta)
   circle(-x2, y2, r2, :magenta)
   circle(0, r2 - 2r1, r2, :magenta)
   if more
       point(0, 0, " O", :red)
       point(0, r1, " r1　甲", :green)
       point(x2, y2, "(x2, y2)　乙", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

算額（その26 ）

算額（その26）

福島県田村郡三春町馬場 三春大神宮 明治 9 年
http://www.wasan.jp/fukusima/miharudaijingu.html

図のように，外円の中に 5 種類の半径を持つ円（木，火，土，金，水）が配置されている。「土」の径を 1 寸とすると，「金」の径は如何ほどか。

以下のように記号を定め，方程式を立てる。

外円の半径を 1 とする。「水」の中心の x 座標を x4，「金」の中心座標を (x5, r5) と置く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, r4::positive, r5::positive, x4::positive, x5::positive;

eq1 = x5^2 + (1 - r1 - r5)^2 - (r1 + r5)^2; # 木-金
eq2 = x5^2 + (r5 -(1 - 2r1 - r2))^2 - (r2 + r5)^2; # 火-金
eq3 = (x5 - x4)^2 + r5^2 - (r4 + r5)^2; # 水-金
eq4 = x4^2 + (r2 + 2r3)^2 - (r2 + r4)^2; # 火-水
eq5 = x4^2 + r3^2 - (r3 + r4)^2; # 土-水 
eq6 = 2(r1 + r2 + r3) - 1; # 半径の合計に関する制約
eq7 = x5^2 + r5^2 - (1 - r5)^2; # 金が外円に内接すること

しかし，なぜか solve() では解けない。

そこで，これも nlsolve() で解く。

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=1e-14)
     v = r.zero[1]
  else
     r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
     v = r.zero
  end
  return v, r.f_converged
end;

function h(u)
  r1, r2, r3, r4, r5, x4, x5 = u
  return [
       x5^2 + (1 - r1 - r5)^2 - (r1 + r5)^2, # 木-金
       x5^2 + (r5 -(1 - 2r1 - r2))^2 - (r2 + r5)^2, # 火-金
       (x5 - x4)^2 + r5^2 - (r4 + r5)^2, # 水-金
       x4^2 + (r2 + 2r3)^2 - (r2 + r4)^2, # 火-水
       x4^2 + r3^2 - (r3 + r4)^2, # 土-水 
       2(r1 + r2 + r3) - 1, # 半径の合計に関する制約
       x5^2 + r5^2 - (1 - r5)^2 # 金が外円に内接すること
  ];
end;

iniv = [0.25, 0.2, 0.1, 0.2, 0.4, 0.2, 0.5]
res = nls(h, ini=iniv)[1];
r1, r2, r3, r4, r5, x4, x5 = res;

name = ["r1", "r2", "r3", "x3", "r4", "r5", "x4", "x5"]
for i = 1:7
  println("$(name[i]) = $(res[i])")
end

   r1 = 0.2795863507293534
   r2 = 0.1773223386410461
   r3 = 0.04309131062960051
   x3 = 0.14151590982250398
   r4 = 0.3602068246353233
   r5 = 0.1795075619334039
   x4 = 0.5287592559278309

r5/r3

   8.359152213576166

答え：「金」の径は，約 8.36 寸である。

算額には「9 寸」と描かれている。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2, r3, r4, r5, x4, x5 = [0.27958635072935356, 0.177322338641047, 0.04309131062959942, 0.14151590982250395, 0.3602068246353232, 0.179507561933404, 0.528759255927831]   
   println("r1 = $r1;  r2 = $r2;  r3 = $r3")
   println("r4 = $r4;  r5 = $r5")
   println("x4 = $x4;  x5 = $x5")
   plot()
   circle(0, 1-r1, r1, :green)
   circle(0, r1-1, r1, :green)
   circle(0, r2+2r3, r2, :magenta)
   circle(0, 2r1+r2-1, r2, :magenta)
   circle(0,  r3, r3, :red)
   circle(0, -r3, r3, :red)
   circle( x4, 0, r4, :black)
   circle(-x4, 0, r4, :black)
   circle( x5,  r5, r5, :brown)
   circle(-x5,  r5, r5, :brown)
   circle( x5, -r5, r5, :brown)
   circle(-x5, -r5, r5, :brown)
   if more
       point(0.15, 0.9, "木", :green, mark=false)
       point(0.05, 0.4, "火", :magenta, mark=false)
       point(0, -0.01, "土", :red, :center, mark=false)
       point(0.6, 0.6, "金", :brown, mark=false)
       point(0.2, 0.1, "水", :black, mark=false)
       point(0, 1-r1, "1-r1 ", :green, :right)
       point(0, r2+2r3, "r2+2r3 ", :magenta, :right)
       point(0, r3, "r3  ", :red, :right)
       point(x4, 0, "x4", :black)
       point(x5, r5, " (x5,r5)", :brown)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
   circle(0, 0, 1, :black,)
end;

算額（その25）

算額（その25）

福島県田村郡三春町御木沢 厳島神社 明治 18 年
http://www.wasan.jp/fukusima/miharuitukusima2.html

図のように，外円にひし形の相対する 2 つの頂点が内接し，そのひし形の内部に大小の円が互いに外接している。大円の径が 4 寸のとき，小円の径はいかほどか。

以下のように記号を定め，方程式を立てる。

大円の半径を 1 とする。
x0 は，大円と小円の共通する接線（接点をそれぞれ (x1,y1), (x2, y2) とする）が x 軸と交わるときの x 座標である。

using SymPy

@syms x1::positive, y1::positive, x2::positive, y2::positive, r2::positive, x0::positive;

eq1 = (2 - y1) / (-x1) * y1 / (x1 - 1) + 1; # 接線と半径が直交する
eq2 = (x0-2 - r2) / (x0 - 1) - y2/y1; # 相似関係
eq3 = (x0-2 - r2) / (x0 - 1) - r2; # 相似関係
eq4 = (x1 - 1)^2 + y1^2 - 1; # (x1, y1) が大円の円周上にある
eq5 = y1/(x0-x1) - y2/(x0 - x2); # 相似関係
eq6 = y1 / (x0 - x1) - 2/x0; # 相似関係

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r2, x1, y1, x2, y2, x0))

   1-element Vector{NTuple{6, Sym}}:
    (1/4, 8/5, 4/5, 12/5, 1/5, 8/3)

name = ["r2", "x1", "y1", "x2", "y2", "x0"]
for (name0, res0) in zip(name, res[1])
   println("$name0 = $res0 = $(res0.evalf())")
end

   r2 = 1/4 = 0.250000000000000
   x1 = 8/5 = 1.60000000000000
   y1 = 4/5 = 0.800000000000000
   x2 = 12/5 = 2.40000000000000
   y2 = 1/5 = 0.200000000000000
   x0 = 8/3 = 2.66666666666667

答え：小円の半径は大円の半径の 1/4 である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2, x1, y1, x2, y2, x0 = (1/4, 8/5, 4/5, 12/5, 1/5, 8/3)    
   println("r2 = $r2")
   println("x1 = $x1;  y1 = $y1;  x2 = $x2;  y2 = $y2")
   println("x0 = $x0")
   plot()
   circle(1, 0, 1, :green)
   circle(-1, 0, 1, :green)
   circle(2 + r2, 0, r2, :blue)
   circle(-2 - r2, 0, r2, :blue)
   plot!([0, x0, 0, -x0, 0], [2, 0, -2, 0, 2], color=:black, linewidth=0.5)
   if more
       point(0, 2, "2 ", :black, :right, :bottom)
       point(x1, y1, "(x1,y1) ", :green, :right, :top)
       point(x2, y2, "(x2,y2)", :blue, :left, :bottom)
       point(2 + r2, 0, "2+r2", :blue, :center, :top)
       point(x0, 0, "  x0", :magenta, :center)
       point(1, 0, "1   ", :green, :top)
       point(2, 0, "2   ", :green, :top)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
   circle(0, 0, 2, :red)
end;

2.675 を四捨五入して，小数点以下2桁で表示せよ

2.675 を小数点以下3桁目を四捨五入せよ。

2.68 となることを予想しているのだが，

★ Python

>>> round(2.675, 2)

2.67

で，decimal パッケージの Decimal を使いましょう...ってさ

★ R

> round(2.675, 2)

[1] 2.67

★ Julia

julia> round(2.675, digits=2)

2.68

★ AWK には round 関数はないのだが int(2.675*100+0.5)/100 というような関数を定義すればよい。

$ awk 'BEGIN{print int(2.675*100+0.5)/100}'

2.68

まあ，完全かどうかは虱潰しにでも調べるしかないか?

算額（その24）

算額（その24）

大きな円の中に，半径の異なる 4 種類の円が，図のように配置されている。各円の半径を求めよ。

岩手県東磐井郡藤沢町 藤勢寺 元治2年
http://www.wasan.jp/iwate/fujise.html

外側の円の半径を 2 とし，内部の円の中心座標と半径を小さい順に (0, y1, r1), (0, y2, r2), (x3, y3, r3), (1, 0, 1) と定める。

using SymPy

@syms r1::positive, r2::positive, r3::positive,
       y1::positive, y2::positive, x3::positive, y3::positive;

未知数は r1, r2, r3, x3, y1, y2, y3 の 7 個で，以下の 7 つの方程式を立てて解く。

eq1 = x3^2 + (y1 - y3)^2 - (r1 + r3)^2;
eq2 = x3^2 + (y3 - y2)^2 - (r2 + r3)^2;
eq3 = 1 + y2^2 - (1 + r2)^2 |> expand
eq4 = (1 - x3)^2 + y3^2 - (1 + r3)^2;
eq5 = y1 - y2  - r1 - r2;
eq6 = y1 - 2 + r1;
eq7 = x3^2 + y3^2 - (2 - r3)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x3, y1, y2, y3));

解は以下のごとくである。

name = ["r1", "r2", "r3", "x3", "y1", "y2", "y3"]
for i = 1:7
   println("$(name[i]) = $(res[1][i]) = $(res[1][i].evalf())")
end

   r1 = 8/11 - 2*sqrt(5)/11 = 0.320714913181856
   r2 = 16*sqrt(5)/209 + 46/209 = 0.391277931291850
   r3 = 7/11 - sqrt(5)/11 = 0.433084729318201
   x3 = 1/11 + 3*sqrt(5)/11 = 0.700745812045397
   y1 = 2*sqrt(5)/11 + 14/11 = 1.67928508681814
   y2 = 68/209 + 60*sqrt(5)/209 = 0.967292242344437
   y3 = 2/11 + 6*sqrt(5)/11 = 1.40149162409079

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x3, y1, y2, y3) = (8/11 - 2*sqrt(5)/11, 16*sqrt(5)/209 + 46/209, 7/11 - sqrt(5)/11, 1/11 + 3*sqrt(5)/11, 2*sqrt(5)/11 + 14/11, 68/209 + 60*sqrt(5)/209, 2/11 + 6*sqrt(5)/11)
   println("r1 = $r1,  r2 = $r2,  r3 = $r3")
   println("y1 = $y1,  y2 = $y2,  y3 = $y3")
   plot()
   circle(0, y1, r1, :green)
   circle(0, -y1, r1, :green)
   circle(0, y2, r2, :blue)
   circle(0, -y2, r2, :blue)
   circle( x3,  y3, r3, :magenta)
   circle(-x3,  y3, r3, :magenta)
   circle( x3, -y3, r3, :magenta)
   circle(-x3, -y3, r3, :magenta)
   circle( 1, 0, 1)
   circle(-1, 0, 1)
   if more
       point(0, y1, "(0,y1,r1)", :brown, :center)
       point(0, y2, "(0,y2,r2)", :blue, :center)
       point(x3, y3, "(x3,y3,r3)", :magenta, :center)
       point(1, 0, "(1,0,1)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.5)
   end
   circle(0, 0, 2, :black)
end;

nlsolve() の解も付しておく。

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=1e-14)
      v = r.zero[1]
   else
      r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
      v = r.zero
   end
   return v, r.f_converged
end;

function h(u)
   r1, r2, r3, x3, y1, y2, y3 = u
   return [
       x3^2 + (y1 - y3)^2 - (r1 + r3)^2,
       x3^2 + (y3 - y2)^2 - (r2 + r3)^2,
       1^2 + y2^2 - (1 +r2)^2,
       (1 - x3)^2 + y3^2 - (1 + r3)^2,
       (y1 - y2)  - (r1+ r2),
       y1 - 2 + r1,
       x3^2 + y3^2 - (2 - r3)^2
   ]
end;

iniv = [0.354644124017785, 0.33761490117551013, 0.4639474199979425, 0.6081577400061725, 1.3714562740209328, 0.8883769604436443, 1.4105324143047882]  # 初期値

res = nls(h, ini=iniv);

name = ["r1", "r2", "r3", "x3", "y1", "y2", "y3"]
for i = 1:7
   println("$(name[i]) = $(res[1][i])")
end

   r1 = 0.32071491318185646
   r2 = 0.3912779312918499
   r3 = 0.4330847293182009
   x3 = 0.7007458120453972
   y1 = 1.6792850868181435
   y2 = 0.9672922423444372
   y3 = 1.4014916240907944

 

 

算額（その23）

算額（その23）

直角三角形に接する大円，小円の径がそれぞれ 3，2 である。直角三角形の外形を描け。

山形県天童市山元 若松観音堂
http://www.wasan.jp/yamagata/wakamatu2.html

ほかに必要とする座標は図に示すとおりである。

using SymPy

@syms x::positive, x1::positive, y1::positive, x2::positive, y2::positive, x3::positive;

小円の中心座標を(x2, 2) とする。大円と小円が接していることから，

eq1 = 1 + (x2 - 3)^2 - 25
x2 = solve(eq1)[1]
x2 |> println
x2.evalf() |> println # x2

   3 + 2*sqrt(6)
   7.89897948556636

△ABC, △DEC が相似なので，(x3 - x2) / (x3 - 3) = 2//3

x3 = solve((x3 - x2)/ (x3 - 3) - 2//3, x3)[1]
x3 |>  println # x3

   3 + 6*sqrt(6)

大円と (x1, y1) で接し，(x3, 0) を通る直線について，
円周上に存在するので，

eq2 = (x1 - 3)^2 + (y1 - 3)^2 - 9;
eq3 = (x3 - x1)^2 + y1^2 - (x3 - 3)^2;

x1, y1 = solve([eq2, eq3], (x1, y1))[1] # x1, y1

   (12*sqrt(6)/25 + 3, 144/25)

△FGC，△IOC の相似により，

I = y1*x3/(x3-x1) |> simplify |> together |> println # intercept

   12*(sqrt(6) + 12)/23

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x3 = 3 + 6sqrt(6)
   x2 = 3 + 2sqrt(6)
   (x1, y1) = (12*sqrt(6)/25 + 3, 144/25)
   println("x2 = $x2")
   println("x1 = $x1;  x2 = $x2")
   intercept = 12(sqrt(6) + 12)/23
   println("x3 = $x3;  intercept = $intercept")
   plot([0, x3, 0, 0], [0, 0, intercept, 0], color=:black, linewidth=0.5)
   circle(3, 3, 3, :green)
   circle(x2, 2, 2, :magenta)
   if more
       point(3, 3, "A:(3,3)", :green, :center)
       point(3, 0, "\nB:(3,0)", :green, :bottom)
       point(x3, 0, " C:(x3,0)", :black, :left, :bottom)
       point(x2, 2, "D:(x2,2)", :magenta, :center)
       point(x2, 0, "\nE:(x2,0)", :magenta, :bottom)
       point(x1, y1, " F:(x1,y1)", :green, :left, :bottom)
       point(x1, 0, "\n\nG:(x1,y1)", :green, :left, :bottom)
       point(0, 0, " O", :black, :left, :bottom)
       point(0, intercept, " I(0,intercept)", :black, :left, :bottom)
   end
end;

算額（その22）

算額（その22）

バージョンアップ　https://blog.goo.ne.jp/r-de-r/e/b049f677e607ba6a92b04bb805ae3b0f

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

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

東京都 住吉神社 嘉永8年(1851) 小円 2 個のない簡易版
http://www.wasan.jp/tokyo/sumiyosi.html

ほかに必要とする座標は図に示すとおりである。

using SymPy

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

以下の 5 式を立て，solve() で解こうとするも，一向に計算が終了しない。

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

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

そこで，nlsolve() で数値解を求めることにする。

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=1e-14)
      v = r.zero[1]
   else
      r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
      v = r.zero
   end
   return v, r.f_converged
end;

function h(u)
   r1, r2, r3, x2, y3 = u
   return [
       r3 / (sqrt(3) - y3) - sin(pi/6),
       r2　/　(1 - x2) - tan(pi/6),
       r1^2 + (y3 - r1)^2 - (r1 + r3)^2,
       x2^2 + (y3 - r2)^2 - (r2 + r3)^2,
       (x2 - r1)^2 + (r2 - r1)^2 - (r1 + r2)^2
   ]
end;

iniv = [0.20, 0.3, 0.5, 0.5, 0.8]  # 初期値

nls(h, ini=iniv)

   ([0.15054780115874736, 0.2613764437119185, 0.4830337288182606, 0.5472827195892904, 0.7659833499323558], true)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x2, y3) = [0.1505478011587473, 0.2613764437119185, 0.4830337288182606, 0.5472827195892904, 0.7659833499323557]
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  x2 = $x2;  y3 = $y3")
   plot([-1, 1, 0, -1], [0, 0, sqrt(3), 0], color=:black, linewidth=0.25)
   circle(0, sqrt(3)-2r3, r3)
   circle(r1, r1, r1, :blue)
   circle(-r1, r1, r1, :blue)
   circle(1-sqrt(3)*r2, r2, r2, :green)
   circle(sqrt(3)*r2-1, r2, r2, :green)
   if more
       point(r1, r1, "(r1,r1)", :blue, :center)
       point(x2, r2, "(x2,r2)", :green, :center)
       point(0, y3, " y3", :red, :left)
       point(0, y3+r3, "\ny3+r3", :red, :left)
       point(0, sqrt(3), " √3", :red, :left)
       vline!([0], linewidth=0.5)
   end
end;

算額（その21）

算額（その21）

半径が同じ 2 つの半円と，大円，中円，小円がが図のように配置されている。それぞれの円の半径を求めよ。

福島県田村郡三春町御木沢 厳島神社 明治 18 年
http://www.wasan.jp/fukusima/miharuitukusima3.html

半円の半径を 1 とする。そのほかに必要とする座標は図に示すとおりである。

using SymPy

@syms r1::positive, r2::positive, y1::positive;

以下の 3 式を立て，解く

eq1 =(1 - r1)^2 + y1^2 - (1 + r1)^2;
eq2 = r2^2 + y1^2 - (1 - r2)^2;
eq3 = y1^2 + y1^2 - 1;

res = solve([eq1, eq2, eq3], (r1, r2, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (1/8, 1/4, sqrt(2)/2)

name = ["r1", "r2", "y1"]
for i in 1:3
   println("$(name[i]) = $(simplify(res[1][i])) = $(res[1][i].evalf())")
end

   r1 = 1/8 = 0.125000000000000
   r2 = 1/4 = 0.250000000000000
   y1 = sqrt(2)/2 = 0.707106781186548

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y1) = (1/8, 1/4, sqrt(2)/2)
   println("r1 = $r1,  r2 = $r2,  y1 = $y1")
   plot()
   circle(0, 0, 1, beginangle=0, endangle=180)
   circle(0, 2y1, 1, beginangle=180, endangle=360)
   circle(0, y1, y1, :blue)
   circle( r2, y1, r2, :green)
   circle(-r2, y1, r2, :green)
   circle(1 - r1, y1, r1, :brown)
   circle(r1 - 1, y1, r1, :brown)
   plot!([-1, 1, 1, -1, -1], [0, 0, 2y1, 2y1, 0], color=:black, linewidth=0.25)
   if more
       point(0, y1, "y1", :blue, :right)
       point(r2, y1, "(r2,y1)", :green, :center)
       point(y1, y1, "(y1,y1)", :red, :center)
       point(1 - r1, y1, "(1-r1,y1)", :brown, :left)
       hline!([y1], color=:black, linewidth=0.25)
       vline!([0], color=:black, linewidth=0.25)
   end
end;

算額（その20）

算額（その20）

4 種類の円と 1 つの半円が図のように配置されている。それぞれの円の半径を求めよ。

福島県田村郡三春町御木沢 厳島神社 明治 18 年
http://www.wasan.jp/fukusima/miharuitukusima1.html

半円の半径を 1 とする。そのほかに必要とする座標は図に示すとおりである。

using SymPy

@syms r1::positive, r2::positive, x1::positive, x2::positive, y2::positive;

以下の 5 式を立て，解く

half = 1//2;
eq1 =x1^2 + r1^2 - (1 - r1)^2;
eq2 = x2^2 + y2^2 - (1 + r2)^2;
eq3 = x2^2 +(1 + r1 - y2)^2 - (r1 + r2)^2;
eq4 = x2^2 + (y2 - half - r1)^2 - (half + r1 - r2)^2;
eq5 = x1^2 + half^2 - (half + 2r1)^2;

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x1, x2, y2))

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

name = ["r1", "r2", "x1", "x2", "y2"]
for i in 1:5
   println("$(name[i]) = $(simplify(res[1][i])) = $(res[1][i].evalf())")
end

   r1 = -1/2 + sqrt(2)/2 = 0.207106781186548
   r2 = 1/6 = 0.166666666666667
   x1 = sqrt(2 - sqrt(2)) = 0.765366864730180
   x2 = sqrt(4 - 2*sqrt(2))/3 = 0.360797400097465
   y2 = 1/6 + 2*sqrt(2)/3 = 1.10947570824873

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x1, x2, y2) = (-1/2 + sqrt(2)/2, 1/6, sqrt(2 - sqrt(2)), sqrt(4/9 - 2*sqrt(2)/9), 1/6 + 2*sqrt(2)/3)
   println("r1 = $r1,  r2 = $r2,  x1 = $x1,  x2 = $x2,  y2 = $y2")
   plot()
   circle(0, 0, 1, beginangle=0, endangle=180)
   circle(0, 1 + r1, r1, :blue)
   circle( x1, r1, r1, :blue)
   circle(-x1, r1, r1, :blue)
   circle(0, 0.5, 0.5, :brown)
   circle(0, r1 + 0.5, r1 + 0.5, :magenta)
   circle( x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(x1, r1, "(x1,r1)", :blue, :center)
       point(0, 0.5, "1/2 ", :brown, :right)
       point(0, 0.5+r1, "r1+1/2 ", :magenta, :right)
       point(0, 1, "1 ", :black, :right)
       point(0, 1+r1, "1+r1 ", :brown, :right)
       point(0, 1+2r1, "1+2r1 ", :black, :center, :bottom)
       point(x2, y2, "(x2,y2)", :green, :center)
       plot!([x2, x2], [y2, y2+r2],color=:green, linewidth=0.25)
       annotate!(x2, y2+r2/2, text("r2 ", 10, :green, :right))
       vline!([0], color=:black, linewidth=0.25)
   end
end;

算額（その19）

算額（その19）

円の中に大小 2 種類の円が入っている。外円の径を 17寸としたとき，小円の半径を求めよ。

福島県楢葉町 波倉稲荷神社
http://www.wasan.jp/fukusima/hakura.html

必要とする座標は図に示すとおりである。
大円の半径: r1, 小円の半径: r2, 大円・小円の中心座標 (x, y1), (x, y2)

using SymPy

@syms r1::positive, r2::positive, x::positive, y1::negative, y2::positive;

以下の 5 式を立て，解く。簡単な式もあるがそのまま書く。

eq1 = x^2 + (1 - r1 - y2)^2 - (r1 + r2)^2;
eq2 = x^2 +(1-3r1)^2 - (1 - r1)^2;
eq3 = x^2 + (1-3r1-(r1-1))^2 - 4r1^2;
eq4 = y1 - 1 + 3r1;
eq5 = y1 + r1 + r2 - y2;

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x, y1, y2))

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

for i = 1:5
   simplify(res[1][i]) |> println
end

   3/2 - sqrt(5)/2
   -2 + sqrt(5)
   sqrt(-22 + 10*sqrt(5))
   -7/2 + 3*sqrt(5)/2
   -4 + 2*sqrt(5)

r1 = (3 - √5)/2
r2 = -2 + √5
x  = √(-22 + 10√5)
y1 = (-7 + 3√5)/2
y2 = -4 + 2√5

小円の径は外円の径を 17 寸とすると，17*(√5 - 2) = 4.013155617496427 寸である。

算額では，「答 4 寸」と書いている。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x, y1, y2) = (3/2 - sqrt(5)/2, -2 + sqrt(5), sqrt(-22 + 10*sqrt(5)), -7/2 + 3*sqrt(5)/2, -4 + 2*sqrt(5))
   println("r1 = $r1,  r2 = $r2,  y1 = $y1,  y2 = $y2")
   plot()
   circle(0, 0, 1)
   circle(0, 1 - r1, r1, :blue)
   circle( x, y1, r1, :brown)
   circle(-x, y1, r1, :brown)
   circle( x, y2, r2, :green)
   circle(-x, y2, r2, :green)
   circle(0, r1 - 1, r1, :brown)
   hline!([1 - 2r1], color=:black, linewidth=0.5)
   vline!([0], color=:black, linewidth=0.5)
   if more
       point(0, 1 - r1, "(0,1-r1)", :blue, :right)
       point(x, y2, "(x,y2)", :green)
       point(x, y1, "(x,y1)", :brown)
       point(0, r1 - 1, "(0,-1+r1)", :brown)
       plot!([x, x], [y2, y2 - r2], color=:green)
       point(x, y2 - r2, "(x,y2-r2)", :green, :center)
       hline!([0, 1 - 2r1], linestyle=:dot, linewidth=0.25)
   end
end;

算額（その18）

算額（その18）

大きな 2 つの円が同じ直線に接し，なおかつ互いに接している。更に 3 種類の半径を持つ円が図のように配置されている。
それぞれの円の半径を求めよ。

福島県伊達郡梁川町 八幡神社 文化14年
http://www.wasan.jp/fukusima/fyahata.html

一番大きい 2 つの円の半径を1 とする。3 種類の円の半径の小さい方から r1, r2, r3 とし，半径 r2 の円の中心座標を (x2, r2) とする。

using SymPy

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

以下の 4 式を立て，解く。

eq1 = x2^2 + (r2 - r1)^2 - (r1 + r2)^2;
eq2 = (1 - x2)^2 + (1 - r2)^2 - (1 + r2)^2;
eq3 = 1^2 + (1 - (r3 + 2r1))^2 - (1 + r3)^2;
eq4 = x2^2 + ((r3 + 2r1) - r2)^2 - (r2 + r3)^2;

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

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

r1 = (3 - √5)/8
r2 = (7 - 3√5)/2
r3 = (5 + √5)/40
x2 = -2 + √5

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x2) = (3/8 - sqrt(5)/8, 7/2 - 3*sqrt(5)/2, sqrt(5)/40 + 1/8, -2 + sqrt(5))
   println("r1 = $r1,  r2 = $r2,  r3 = $r3,  x2 = $x2")
   plot()
   circle(1, 1, 1)#, beginangle=170, endangle=280)
   circle(-1, 1, 1)#, beginangle=260, endangle=370)
   circle(0, r1, r1, :blue)
   circle(x2, r2, r2, :brown)
   circle(-x2, r2, r2, :brown)
   circle(0, r3 + 2r1, r3, :green)
   hline!([0])
   if more
       point(1, 1, "(1,1)")
       point(0, r1, "r1", :blue, :right)
       point(0, r3 + 2r1, "r3+2r1", :green, :right)
       point(x2, r2, "(r2,x2)", :brown, :top)
       vline!([0], xlims=(-0.5, 1.05), ylims=(-0.05, 1.05))
   else
       plot!(xlims=(-2.05, 2.05), ylims=(-0.05, 2.05))
   end
end;

 

算額（その17）

算額（その17）

円の内部に大きい円が 6 個，小さい円が 3 個ある。それぞれの円の半径を求めよ。

千葉県木更津市 高柳不動
http://www.wasan.jp/chiba/takayanagi.html

右上の大きな円の半径と中心座標を r2, (r2, y2) とする。
一番上の小さな円の半径と中心座標を r1, (0, y1) とする。

using SymPy

@syms r1::positive, y1::positive, r2::positive, y2::positive;

r2, y2 は方程式を立てなくてもわかるが一応形式的に。

eq1 = y2 - (1 - r2)*sind(Sym(60));
eq2 = r2 - (1 - r2)*cosd(Sym(60));

先の円と接する小さな塩の関係から以下の式を得る。

eq3 = (y2 - y1)^2 + r2^2 - (r1 + r2)^2;

一番上にある小さな円の半径を 1，中心座標を (0, y1) とすると，

eq4 = y1 * sind(Sym(60)) - r1;

方程式を解き，解を求める。

res = solve([eq1, eq2, eq3, eq4], (r1, y1, r2, y2))

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

r1 = 2*sqrt(2) + 3 は外円より大きいので不適切解である。

回答：

大きな円の半径 1/3
小さな円の半径 3 - 2√2

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, y1, r2, y2) = (3 - 2*sqrt(2), -4*sqrt(6)/3 + 2*sqrt(3), 1/3, sqrt(3)/3)
   println("r1 = $r1,  y1 = $y1,  r2 = $r2,  y2 = $y2")
   plot()
   circle(0, 0, 1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(r2, -y2, r2, :blue)
   circle(-r2, -y2, r2, :blue)
   circle(1-r2, 0, r2, :blue)
   circle(r2-1, 0, r2, :blue)
   circle(0, y1, r1, :magenta)
   circle(r1, y1 - 2r1*sind(60), r1, :magenta)
   circle(-r1, y1 - 2r1*sind(60), r1, :magenta)
   if more
       point(r2, y2, " (r2,y2)", :blue)
       plot!([0, r2], [y2, y2])
       point(0, y2, "y2 ", :magenta, :right)
       point(0, y1, "y1 ", :magenta, :right)
       point(r1, 0, " r1", :magenta, :left, :top)
       plot!([0, r1], [y1 - 2r1*sind(60), y1 - 2r1*sind(60)])
       hline!([0])
       vline!([0])
   end
end;

 

算額（その16）

算額（その16）

図のように円が十字に仕切られており，内部に 4 つの円がある。大きい順に甲円，乙円，丙円，丁円とする。乙円，丙円，丁円の径はそれぞれ 4 寸，3 寸，2 寸である。甲円の径はいくつか。

福島県福島市上湯町 薬師堂
http://www.wasan.jp/fukusima/kamiyuyakusi.html

「甲径は 8 寸 4 分である」としているが，もう少し正確には 8.4748 寸（図にすると，外接円との間に隙間がある）。

• 十字の仕切り線の好転を原点とする。
• 甲円の半径を r とすれば，中心座標は (r, -r) である。
• 外円の半径と中心座標を R, (x, y) とする。

以下のように方程式を立て，それを解く。

using SymPy

@syms x::positive, y::negative, r::positive, R::positive;

eq1 = (-4 - x)^2 + (-4 - y)^2 - (R - 4)^2;
eq2 = ( 3 - x)^2 + ( 3 - y)^2 - (R - 3)^2;
eq3 = (-2 - x)^2 + ( 2 - y)^2 - (R - 2)^2;
eq4 = ( r - x)^2 + (-r - y)^2 - (R - r)^2;

res = solve([eq1, eq2, eq3, eq4], (r, x, y, R))

   1-element Vector{NTuple{4, Sym}}:
    (-5*sqrt(97)/24 + 3/8 + sqrt(2)*sqrt(1395*sqrt(97) + 15941)/24, 5*sqrt(97)/24 + 19/8, -59/16 - 5*sqrt(97)/16, 91/16 + 35*sqrt(97)/48)

rem = ["r", "x", "y", "R"]
for i in 1:4
   println("$(rem[i]) = $(res[1][i].evalf())")
   println(res[1][i], "\n")
end

   r = 8.47480963363268
   -5*sqrt(97)/24 + 3/8 + sqrt(2)*sqrt(1395*sqrt(97) + 15941)/24

   x = 4.42684537537419
   5*sqrt(97)/24 + 19/8

   y = -6.76526806306128
   -59/16 - 5*sqrt(97)/16

   R = 12.8689588138097
   91/16 + 35*sqrt(97)/48

r を取り出すのは res[1][1]。sympy.sqrtdenest で二重根号を外してから expand すると簡単な式になる。

res[1][1] |> sympy.sqrtdenest |> expand |> println

   sqrt(97)/6 + 41/6

(41 + √97) / 6

   8.474809633632683

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, y, r, R) = (5*sqrt(97)/24 + 19/8,
                   -59/16 - 5*sqrt(97)/16,
                   -5*sqrt(97)/24 + 3/8 + sqrt(2)*sqrt(1395*sqrt(97) + 15941)/24,
                   91/16 + 35*sqrt(97)/48)
   println("rx = $x,  y = $y,  r = $r,  R = $R")
   plot()
   circle(-4, -4, 4, :blue)
   circle( 3,  3, 3, :blue)
   circle(-2,  2, 2, :blue)
   circle( r, -r, r, :blue)
   circle( x,  y, R, :green)
   if more
       point(-4, -4, " 乙", :blue)
       point( 3,  3, " 丙", :blue)
       point(-2,  2, " 丁", :blue)
       point( r, -r, " 甲(r,-r)", :blue)
       point( x,  y, " 外円(x,y)", :green)
   end
   hline!([0])
   vline!([0])
end;