算額（その182）

算額（その182）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（82）
長野県長野市篠ノ井 布制神社 文化3年(1806)

外円の中に長さ12寸の弦を引く。弦の上に甲乙丙円を置く。甲円の径が2寸5分のとき，丙円の径はいかほどか。

外円，甲円，乙円，丙円の半径を R, r3 = 25, r2, r1 とする。乙円，丙円の x 座標を x2, x1 とする。

using SymPy
@syms R::positive, x1::positive, r1::positive, x2::positive, r2::positive;

r3 = 25
chord = 120

eq0 = (R - 2r3)^2 + chord^2 - R^2
eq1 = (x2 - x1)^2 + (r2 - r1)^2 - (r1 + r2)^2
eq2 = x2^2 + (r3 - r2)^2 - (r3 + r2)^2
eq3 = x1^2 + (r3 - r1)^2 - (r3 + r1)^2
eq4 = x2^2 + (R - 2r3 + r2)^2 - (R - r2)^2

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

   2-element Vector{NTuple{5, Sym}}:
    (169, 144/25, 3600/169, 24, 600/13)
    (169, 3600, 3600/169, 600, 600/13)

元の単位では丙円の直径は = 144/25 = 5.76000000000000

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2, x1, x2) = res[1]
   println("丙円の半径は = $(r1) = $(r1.evalf())")
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r3, r3)
   circle(x1, R - 2r3 + r1, r1, :blue)
   circle(x2, R - 2r3 + r2, r2, :green)
   circle(-x2, R - 2r3 + r2, r2, :green)
   segment(-120, R - 2r3, 120, R - 2r3)
   if more == true
       point(0, R, "R ", :green, :right)
       point(0, R - r3, "R-r3 ", :green, :right)
       point(0, R - 2r3, "R-2r3 ", :green, :right)
       point(x1, R - 2r3, "x1") 
       point(x2, R - 2r3, "x2")
       point(7, R - r3, "甲", :red, mark=false)
       point(x1, R - 2r3, "丙 ", :red, :right, :top, mark=false)
       point(x2, R - 2r3 + r2, "乙", :red, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false, xlims=(-30, 125), ylims=(114, 174))
   end
end;

算額（その181）

算額（その181）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（76）
長野県諏訪郡 諏訪社 文化元年(1801)

一〇 武州不動岡村 不動堂 文化5年(1808)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に天円及び天円の半分，地円，人円の他，甲乙丙丁戊己庚辛壬癸の 10 円を入れる。地円の径を知って各円の径を求めよ。

外円，天円，地円，人円の半径 を 1, r1, r2, r3 とする。地円，人円の中心座標を (x2, y2), (x3, y3) とする。

eq1~ eq7 を解いて r1, r2, r3 および (x2, y2), (x3, y3) を得る。

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive, y2::positive,
       r3::positive, x3::positive, y3::positive;

R = 1
eq1 = x3^2 + (2R - r1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + (R - y3)^2 - (R - r3)^2
eq3 = x2^2 + (R - y2)^2 - (R - r2)^2
eq4 = x3^2 + (y3 - 2R + 3r1)^2 - (r1 + r3)^2
eq5 = x2^2 + (y2 - 2R + 3r1)^2 - (r1 + r2)^2
eq6 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2
eq7 = r1^2 + (3r1 - R)^2 - R^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, x2, y2, r3, x3, y3))

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

次に，人円から始めて，甲円...を得るために eq8~eq10 を解き，漸化式を作る。
2組目が適解である。

@syms r4::positive, x4::positive, y4::positive;
eq8 = x4^2 + (2R - r1 - y4)^2 - (r1 + r4)^2
eq9 = x4^2 + (y4 - R)^2 - (R -r4)^2
eq10 = (x4 - x3)^2 + (y4 - y3)^2 - (r3 + r4)^2
solve([eq8, eq9, eq10], (r4, x4, y4))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-2*r1^(5/2)*x3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*sqrt(r1)*x3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - r1^3*r3^3 - r1^3*r3^2*y3 + r1^3*r3*x3^2 + r1^3*r3*y3^2 - 4*r1^3*r3*y3 + 4*r1^3*r3 + r1^3*x3^2*y3 - 4*r1^3*x3^2 + r1^3*y3^3 - 4*r1^3*y3^2 + 4*r1^3*y3 + r1^2*r3^3 - r1^2*r3^2*y3 + 2*r1^2*r3^2 - r1^2*r3*x3^2 - r1^2*r3*y3^2 + 4*r1^2*r3*y3 - 4*r1^2*r3 + r1^2*x3^2*y3 - 2*r1^2*x3^2 + r1^2*y3^3 - 6*r1^2*y3^2 + 12*r1^2*y3 - 8*r1^2 + r1*r3^3 + r1*r3^2*y3 - r1*r3*x3^2 - r1*r3*y3^2 + 4*r1*r3*y3 - 4*r1*r3 - r1*x3^2*y3 + 4*r1*x3^2 - r1*y3^3 + 4*r1*y3^2 - 4*r1*y3 - r3^3 + r3^2*y3 - 2*r3^2 + r3*x3^2 + r3*y3^2 - 4*r3*y3 + 4*r3 - x3^2*y3 + 2*x3^2 - y3^3 + 6*y3^2 - 12*y3 + 8)/(2*(r1^3*r3^2 + 2*r1^3*r3*y3 + r1^3*y3^2 - r1^2*r3^2 + 2*r1^2*r3*y3 - 4*r1^2*r3 + 4*r1^2*x3^2 + 3*r1^2*y3^2 - 4*r1^2*y3 - r1*r3^2 - 2*r1*r3*y3 + 4*r1*x3^2 + 3*r1*y3^2 - 8*r1*y3 + 4*r1 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4)), (r1^(5/2)*r3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + r1^(5/2)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*r1^(3/2)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - 2*r1^(3/2)*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - sqrt(r1)*r3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + sqrt(r1)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - 2*sqrt(r1)*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*r1^3*r3*x3 + 2*r1^3*x3*y3 - 2*r1^2*r3^2*x3 + 2*r1^2*x3^3 + 2*r1^2*x3*y3^2 - 4*r1^2*x3*y3 + 4*r1^2*x3 - 2*r1*r3^2*x3 - 2*r1*r3*x3 + 2*r1*x3^3 + 2*r1*x3*y3^2 - 6*r1*x3*y3 + 4*r1*x3)/(r1^3*r3^2 + 2*r1^3*r3*y3 + r1^3*y3^2 - r1^2*r3^2 + 2*r1^2*r3*y3 - 4*r1^2*r3 + 4*r1^2*x3^2 + 3*r1^2*y3^2 - 4*r1^2*y3 - r1*r3^2 - 2*r1*r3*y3 + 4*r1*x3^2 + 3*r1*y3^2 - 8*r1*y3 + 4*r1 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4), -sqrt(r1)*x3*sqrt(-(-r3^2 - 2*r3 + x3^2 + y3^2 - 2*y3)*(2*r1*r3 + 2*r1*y3 - 4*r1 - r3^2 + x3^2 + y3^2 - 4*y3 + 4))*(r1 + 1)/(r1^2*r3^2 + 2*r1^2*r3*y3 + r1^2*y3^2 - 2*r1*r3^2 - 4*r1*r3 + 4*r1*x3^2 + 2*r1*y3^2 - 4*r1*y3 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4) + (-r1^2*r3^3 - r1^2*r3^2*y3 + 4*r1^2*r3^2 + r1^2*r3*x3^2 + r1^2*r3*y3^2 + 4*r1^2*r3*y3 + 4*r1^2*r3 + r1^2*x3^2*y3 - 4*r1^2*x3^2 + r1^2*y3^3 + 4*r1^2*y3 - 2*r1*r3^2*y3 - 6*r1*r3^2 - 16*r1*r3 + 2*r1*x3^2*y3 + 10*r1*x3^2 + 2*r1*y3^3 - 2*r1*y3^2 - 8*r1 + r3^3 - r3^2*y3 + 6*r3^2 - r3*x3^2 - r3*y3^2 - 4*r3*y3 + 12*r3 + x3^2*y3 - 2*x3^2 + y3^3 - 2*y3^2 - 4*y3 + 8)/(2*(r1^2*r3^2 + 2*r1^2*r3*y3 + r1^2*y3^2 - 2*r1*r3^2 - 4*r1*r3 + 4*r1*x3^2 + 2*r1*y3^2 - 4*r1*y3 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4)))
    ((2*r1^(5/2)*x3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - 2*sqrt(r1)*x3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - r1^3*r3^3 - r1^3*r3^2*y3 + r1^3*r3*x3^2 + r1^3*r3*y3^2 - 4*r1^3*r3*y3 + 4*r1^3*r3 + r1^3*x3^2*y3 - 4*r1^3*x3^2 + r1^3*y3^3 - 4*r1^3*y3^2 + 4*r1^3*y3 + r1^2*r3^3 - r1^2*r3^2*y3 + 2*r1^2*r3^2 - r1^2*r3*x3^2 - r1^2*r3*y3^2 + 4*r1^2*r3*y3 - 4*r1^2*r3 + r1^2*x3^2*y3 - 2*r1^2*x3^2 + r1^2*y3^3 - 6*r1^2*y3^2 + 12*r1^2*y3 - 8*r1^2 + r1*r3^3 + r1*r3^2*y3 - r1*r3*x3^2 - r1*r3*y3^2 + 4*r1*r3*y3 - 4*r1*r3 - r1*x3^2*y3 + 4*r1*x3^2 - r1*y3^3 + 4*r1*y3^2 - 4*r1*y3 - r3^3 + r3^2*y3 - 2*r3^2 + r3*x3^2 + r3*y3^2 - 4*r3*y3 + 4*r3 - x3^2*y3 + 2*x3^2 - y3^3 + 6*y3^2 - 12*y3 + 8)/(2*(r1^3*r3^2 + 2*r1^3*r3*y3 + r1^3*y3^2 - r1^2*r3^2 + 2*r1^2*r3*y3 - 4*r1^2*r3 + 4*r1^2*x3^2 + 3*r1^2*y3^2 - 4*r1^2*y3 - r1*r3^2 - 2*r1*r3*y3 + 4*r1*x3^2 + 3*r1*y3^2 - 8*r1*y3 + 4*r1 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4)), (-r1^(5/2)*r3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - r1^(5/2)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - 2*r1^(3/2)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*r1^(3/2)*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + sqrt(r1)*r3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - sqrt(r1)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*sqrt(r1)*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*r1^3*r3*x3 + 2*r1^3*x3*y3 - 2*r1^2*r3^2*x3 + 2*r1^2*x3^3 + 2*r1^2*x3*y3^2 - 4*r1^2*x3*y3 + 4*r1^2*x3 - 2*r1*r3^2*x3 - 2*r1*r3*x3 + 2*r1*x3^3 + 2*r1*x3*y3^2 - 6*r1*x3*y3 + 4*r1*x3)/(r1^3*r3^2 + 2*r1^3*r3*y3 + r1^3*y3^2 - r1^2*r3^2 + 2*r1^2*r3*y3 - 4*r1^2*r3 + 4*r1^2*x3^2 + 3*r1^2*y3^2 - 4*r1^2*y3 - r1*r3^2 - 2*r1*r3*y3 + 4*r1*x3^2 + 3*r1*y3^2 - 8*r1*y3 + 4*r1 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4), sqrt(r1)*x3*sqrt(-(-r3^2 - 2*r3 + x3^2 + y3^2 - 2*y3)*(2*r1*r3 + 2*r1*y3 - 4*r1 - r3^2 + x3^2 + y3^2 - 4*y3 + 4))*(r1 + 1)/(r1^2*r3^2 + 2*r1^2*r3*y3 + r1^2*y3^2 - 2*r1*r3^2 - 4*r1*r3 + 4*r1*x3^2 + 2*r1*y3^2 - 4*r1*y3 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4) + (-r1^2*r3^3 - r1^2*r3^2*y3 + 4*r1^2*r3^2 + r1^2*r3*x3^2 + r1^2*r3*y3^2 + 4*r1^2*r3*y3 + 4*r1^2*r3 + r1^2*x3^2*y3 - 4*r1^2*x3^2 + r1^2*y3^3 + 4*r1^2*y3 - 2*r1*r3^2*y3 - 6*r1*r3^2 - 16*r1*r3 + 2*r1*x3^2*y3 + 10*r1*x3^2 + 2*r1*y3^3 - 2*r1*y3^2 - 8*r1 + r3^3 - r3^2*y3 + 6*r3^2 - r3*x3^2 - r3*y3^2 - 4*r3*y3 + 12*r3 + x3^2*y3 - 2*x3^2 + y3^3 - 2*y3^2 - 4*y3 + 8)/(2*(r1^2*r3^2 + 2*r1^2*r3*y3 + r1^2*y3^2 - 2*r1*r3^2 - 4*r1*r3 + 4*r1*x3^2 + 2*r1*y3^2 - 4*r1*y3 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4)))

以下の関数は，現在の円の半径，中心座標から，次の円の半径，中心座標を求める。

function nextcircle(r3, x3, y3)
   r1 = 3/5
   return ((2*r1^(5/2)*x3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - 2*sqrt(r1)*x3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - r1^3*r3^3 - r1^3*r3^2*y3 + r1^3*r3*x3^2 + r1^3*r3*y3^2 - 4*r1^3*r3*y3 + 4*r1^3*r3 + r1^3*x3^2*y3 - 4*r1^3*x3^2 + r1^3*y3^3 - 4*r1^3*y3^2 + 4*r1^3*y3 + r1^2*r3^3 - r1^2*r3^2*y3 + 2*r1^2*r3^2 - r1^2*r3*x3^2 - r1^2*r3*y3^2 + 4*r1^2*r3*y3 - 4*r1^2*r3 + r1^2*x3^2*y3 - 2*r1^2*x3^2 + r1^2*y3^3 - 6*r1^2*y3^2 + 12*r1^2*y3 - 8*r1^2 + r1*r3^3 + r1*r3^2*y3 - r1*r3*x3^2 - r1*r3*y3^2 + 4*r1*r3*y3 - 4*r1*r3 - r1*x3^2*y3 + 4*r1*x3^2 - r1*y3^3 + 4*r1*y3^2 - 4*r1*y3 - r3^3 + r3^2*y3 - 2*r3^2 + r3*x3^2 + r3*y3^2 - 4*r3*y3 + 4*r3 - x3^2*y3 + 2*x3^2 - y3^3 + 6*y3^2 - 12*y3 + 8)/(2*(r1^3*r3^2 + 2*r1^3*r3*y3 + r1^3*y3^2 - r1^2*r3^2 + 2*r1^2*r3*y3 - 4*r1^2*r3 + 4*r1^2*x3^2 + 3*r1^2*y3^2 - 4*r1^2*y3 - r1*r3^2 - 2*r1*r3*y3 + 4*r1*x3^2 + 3*r1*y3^2 - 8*r1*y3 + 4*r1 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4)), (-r1^(5/2)*r3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - r1^(5/2)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - 2*r1^(3/2)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*r1^(3/2)*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + sqrt(r1)*r3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) - sqrt(r1)*y3*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*sqrt(r1)*sqrt(2*r1*r3^3 + 2*r1*r3^2*y3 - 2*r1*r3*x3^2 - 2*r1*r3*y3^2 + 8*r1*r3*y3 - 8*r1*r3 - 2*r1*x3^2*y3 + 4*r1*x3^2 - 2*r1*y3^3 + 8*r1*y3^2 - 8*r1*y3 - r3^4 - 2*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 6*r3^2*y3 + 4*r3^2 + 2*r3*x3^2 + 2*r3*y3^2 - 8*r3*y3 + 8*r3 - x3^4 - 2*x3^2*y3^2 + 6*x3^2*y3 - 4*x3^2 - y3^4 + 6*y3^3 - 12*y3^2 + 8*y3) + 2*r1^3*r3*x3 + 2*r1^3*x3*y3 - 2*r1^2*r3^2*x3 + 2*r1^2*x3^3 + 2*r1^2*x3*y3^2 - 4*r1^2*x3*y3 + 4*r1^2*x3 - 2*r1*r3^2*x3 - 2*r1*r3*x3 + 2*r1*x3^3 + 2*r1*x3*y3^2 - 6*r1*x3*y3 + 4*r1*x3)/(r1^3*r3^2 + 2*r1^3*r3*y3 + r1^3*y3^2 - r1^2*r3^2 + 2*r1^2*r3*y3 - 4*r1^2*r3 + 4*r1^2*x3^2 + 3*r1^2*y3^2 - 4*r1^2*y3 - r1*r3^2 - 2*r1*r3*y3 + 4*r1*x3^2 + 3*r1*y3^2 - 8*r1*y3 + 4*r1 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4), sqrt(r1)*x3*sqrt(-(-r3^2 - 2*r3 + x3^2 + y3^2 - 2*y3)*(2*r1*r3 + 2*r1*y3 - 4*r1 - r3^2 + x3^2 + y3^2 - 4*y3 + 4))*(r1 + 1)/(r1^2*r3^2 + 2*r1^2*r3*y3 + r1^2*y3^2 - 2*r1*r3^2 - 4*r1*r3 + 4*r1*x3^2 + 2*r1*y3^2 - 4*r1*y3 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4) + (-r1^2*r3^3 - r1^2*r3^2*y3 + 4*r1^2*r3^2 + r1^2*r3*x3^2 + r1^2*r3*y3^2 + 4*r1^2*r3*y3 + 4*r1^2*r3 + r1^2*x3^2*y3 - 4*r1^2*x3^2 + r1^2*y3^3 + 4*r1^2*y3 - 2*r1*r3^2*y3 - 6*r1*r3^2 - 16*r1*r3 + 2*r1*x3^2*y3 + 10*r1*x3^2 + 2*r1*y3^3 - 2*r1*y3^2 - 8*r1 + r3^3 - r3^2*y3 + 6*r3^2 - r3*x3^2 - r3*y3^2 - 4*r3*y3 + 12*r3 + x3^2*y3 - 2*x3^2 + y3^3 - 2*y3^2 - 4*y3 + 8)/(2*(r1^2*r3^2 + 2*r1^2*r3*y3 + r1^2*y3^2 - 2*r1*r3^2 - 4*r1*r3 + 4*r1*x3^2 + 2*r1*y3^2 - 4*r1*y3 + r3^2 - 2*r3*y3 + 4*r3 + y3^2 - 4*y3 + 4)))
end;

(r1, r2, x2, y2, r3, x3, y3) = (3/5, 1/10, 3*sqrt(5)/10, 2/5, 3/10, 3*sqrt(5)/10, 4/5)
println("地円の半径 = $r2")
(r, x, y) = (r3, x3, y3)
for i in 1:11
   println("$(Char["人甲乙丙丁戊己庚辛壬癸"...][i])円の半径 = $r")
   (r, x, y) = nextcircle(r, x, y)
end

   地円の半径 = 0.1
   人円の半径 = 0.3
   甲円の半径 = 0.1821257430242036
   乙円の半径 = 0.11134091913935494
   丙円の半径 = 0.07243506027318493
   丁円の半径 = 0.050093053005542496
   戊円の半径 = 0.036425148604841105
   己円の半径 = 0.02756521532478412
   庚円の半径 = 0.02153548717194146
   辛円の半径 = 0.01726349566989466
   壬円の半径 = 0.014134327321380202
   癸円の半径 = 0.011777575328980357

地円の半径は外円の半径の 1/10 である。地円の半径（直径）がわかれば，それぞれの円の半径（直径）がわかる。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x2, y2, r3, x3, y3) = (3/5, 1/10, 3*sqrt(5)/10, 2/5, 3/10, 3*sqrt(5)/10, 4/5)
   println("外円の半径 = 1")
   println("天円の半径 = $r1")
   println("地円の半径 = $r2, 中心座標 = ($x2, $y2)")
   name = "人甲乙丙丁戊己庚辛壬癸"
   R = 1
   plot()
   circle(0, R, R, :black)
   circle(0, 2R - r1, r1)
   circle(0, 2R - 3r1, r1, beginangle=0, endangle=180)
   segment(r1, 2R - 3r1, -r1, 2R - 3r1, :red)
   circle(x2, y2, r2)
   circle(-x2, y2, r2)
   (r, x, y) = (r3, x3, y3)
   for i = 1:11
       println("$(Char[name...][i])円の半径 = $r, 中心座標 = ($x, $y)")
       circle(x, y, r, i)
       circle(-x, y, r, i)
       (r, x, y) = nextcircle(r, x, y)
   end
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

外円の半径 = 1
天円の半径 = 0.6
地円の半径 = 0.1, 中心座標 = (0.6708203932499369, 0.4)
人円の半径 = 0.3, 中心座標 = (0.6708203932499369, 0.8)
甲円の半径 = 0.1821257430242036, 中心座標 = (0.771497027903185, 1.2714970279031843)
乙円の半径 = 0.11134091913935494, 中心座標 = (0.6943295404303229, 1.5546363234425813)
丙円の半径 = 0.07243506027318493, 中心座標 = (0.5965800803642471, 1.71025975890726)
丁円の半径 = 0.050093053005542496, 中心座標 = (0.5127558957652436, 1.799627787977826)
戊円の半径 = 0.036425148604841105, 中心座標 = (0.44570059441936133, 1.8542994055806359)
己円の半径 = 0.02756521532478412, 中心座標 = (0.39242027917805117, 1.8897391387008609)
庚円の半径 = 0.02153548717194146, 中心座標 = (0.34965163365220264, 1.9138580513122398)
辛円の半径 = 0.01726349566989466, 中心座標 = (0.31481828056547234, 1.930946017320418)
壬円の半径 = 0.014134327321380202, 中心座標 = (0.28602320849166196, 1.94346269071449)
癸円の半径 = 0.011777575328980357, 中心座標 = (0.2618869656253128, 1.9528896986840758)

算額（その180）

算額（その180）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（56）
長野県長野市元善町 善光寺 寛政8年(1796)

第6問: 長方形内に 8 個の円を入れる。甲円の径が 17寸7分4厘 のとき，丙円の径を求めよ。

甲円の径を 100 倍して整数で取り扱う。長方形の長辺を 2x，甲円の半径は 887，乙円の半径を r1, 丙円の半径を r2 とする。以下のように連立方程式を立て，解く。

using SymPy
@syms x::positive, r2::positive, r1::positive

eq1 = (x - 887)^2 + (887 - r1)^2 - (r1 + 887)^2
eq2 = (x - 887 - r2)^2 + 887^2 - (887 + r2)^2
eq3 = r2^2 + (1774 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (x, r1, r2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-887*sqrt(5 + 4*sqrt(2)) + 887*sqrt(2) + 887*sqrt(10 + 8*sqrt(2)), -887*sqrt(10 + 8*sqrt(2)) + 887/2 + 2661*sqrt(5 + 4*sqrt(2))/2, -887*sqrt(5 + 4*sqrt(2)) + 887 + 1774*sqrt(2))
    (-887*sqrt(10 + 8*sqrt(2)) + 887*sqrt(2) + 887*sqrt(5 + 4*sqrt(2)), -2661*sqrt(5 + 4*sqrt(2))/2 + 887/2 + 887*sqrt(10 + 8*sqrt(2)), 887 + 1774*sqrt(2) + 887*sqrt(5 + 4*sqrt(2)))

最初の組の解が適解である。
丙円の半径 r2 は -887*sqrt(5 + 4*sqrt(2)) + 887 + 1774*sqrt(2)

直径は

2(-887*sqrt(5 + 4*sqrt(2)) + 887 + 1774*sqrt(2))

   1000.4355208571196

元の単位でいえば「10寸0分0厘あまりあり」ということである。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, r1, r2) = res[1]
   println("丙円の径 = $(2r2) = $(2r2.evalf())")
   println("乙円の径 = $(2r1) = $(2r1.evalf())")
   println("長方形の長辺 = $(2x) = $(2x.evalf())")
   plot([x, x, -x, -x, x], [-1774, 1774, 1774, -1774, -1774], color=:black, lw=0.5)
   circle(x - 887, 887, 887)
   circle(x - 887, -887, 887)
   circle(-x + 887, 887, 887)
   circle(-x + 887, -887, 887)
   circle(0, 1774 - r1, r1, :blue)
   circle(0, -1774 + r1, r1, :blue)
   circle(r2, 0, r2, :orange)
   circle(-r2, 0, r2, :orange)
   if more == true
       point(x - 887, 887, "甲円 (x-887,887)", :green, :top)
       point(0, 1774 - r1, " 1774-r1")
       point(200, 1300, "乙円", mark=false)
       point(r2, 0, " r2")
       point(r2, 300, "丙円", mark=false)
       point(x, 0, " x")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

draw(false)

   丙円の径 = -1774*sqrt(5 + 4*sqrt(2)) + 1774 + 3548*sqrt(2) = 1000.43552085712
   乙円の径 = -1774*sqrt(10 + 8*sqrt(2)) + 887 + 2661*sqrt(5 + 4*sqrt(2)) = 1383.80591988999
   長方形の長辺 = -1774*sqrt(5 + 4*sqrt(2)) + 1774*sqrt(2) + 1774*sqrt(10 + 8*sqrt(2)) = 4907.60603898119

算額（その179）

算額（その179）

長野県上田市別所温泉 北向観音堂 安永7年(1778)
中村信弥「改訂増補 長野県の算額」(29)
http://www.wasan.jp/zoho/zoho.html

大円の中に天円 2 個，地円 2 個，人円 4個を入れる。大円の径は 216 寸である。天, 地, 人, 甲, 乙, 丙, 丁, 戊, 己, 庚, 辛, 壬, 癸, 乾, 兌, 離, 震, 巽, 坎, 艮, 坤の各円の径を求めよ。

大，天，地，人円の半径と中心座標を (0, 0, R), (r1, 0, r1), (0, R-r2, r2), (x3, y3, r3) とする。以下のように連立方程式を立て，解く。

using SymPy

@syms R::positive, r1::positive, r2::positive,
   r3::positive, x3::positive, y3::positive,
   r4::positive, x4::positive, y4::positive,
   r5::positive, x5::positive, y5::positive,
   r6::positive, x6::positive, y6::positive

R = 216 // 2
r1 = R // 2
eq1 = r1^2 + (R - r2)^2 - (r1 + r2)^2  # 天円と地円が外接する
eq2 = (r1 - x3)^2 + y3^2 - (r1 + r3)^2  # 天円と人円が外接する
eq3 = x3^2 + (R - r2 - y3)^2 - (r2 + r3)^2  # 地円と人円が外接する
eq4 = x3^2 + y3^2 - (R - r3)^2;  # 人円が外円に内接する

eq5 = (x4 - x3)^2 + (y3 - y4)^2 - (r3 + r4)^2  # 円と次の円が外接する
eq6 = x4^2 + y4^2 - (R - r4)^2                 # 次の円が外円に内接する
eq7 = (x4 - r1)^2 + y4^2 - (r1 + r4)^2         # 天円と次の円が外接する

eq8 = (x5 - x4)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = x5^2 + y5^2 - (R - r5)^2
eq10 = (x5 - r1)^2 + y5^2 - (r1 + r5)^2

eq11 = (x6 - x5)^2 + (y5 - y6)^2 - (r5 + r6)^2
eq12 = x6^2 + y6^2 - (R - r6)^2
eq13 = (x6 - r1)^2 + y6^2 - (r1 + r6)^2;

eq1 ~ eq4 を解いて r2, r3, x3, y3 を得る。

res = solve([eq1, eq2, eq3, eq4], (r2, r3, x3, y3))

   1-element Vector{NTuple{4, Sym}}:
    (36, 18, 54, 72)

eq1 ~ eq13 を解くと r2, r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6 が得られる（1組目が適解）。

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10,
      eq11, eq12, eq13],
  (r2, r3, x3, y3, r4, x4, y4, r5, x5, y5, r6, x6, y6))

   3-element Vector{NTuple{13, Sym}}:
    (36, 18, 54, 72, 108/11, 864/11, 648/11, 6, 90, 48, 4, 96, 40)
    (36, 18, 54, 72, 108/11, 864/11, 648/11, 6, 90, 48, 108/11, 864/11, 648/11)
    (36, 18, 54, 72, 108/11, 864/11, 648/11, 18, 54, 72, 108/11, 864/11, 648/11)

eq5, eq6, eq7 と eq8, eq9, eq10 を見れば，3つの連立方程式は漸化式になっている。そこで，以下の 3 連立方程式を ri, xi, yi について，現在の円から次の円の半径，と中心座標を求める関数を定義する。

@syms ri, xi, yi, rip1, xip1, yip1
eq01 = (xip1 - xi)^2 + (yi - yip1)^2 - (ri + rip1)^2
eq02 = xip1^2 + yip1^2 - (R - rip1)^2
eq03 = (xip1 - r1)^2 + yip1^2 - (r1 + rip1)^2
res0 = solve([eq01, eq02, eq03], (rip1, xip1, yip1)) 

以下がその関数である。現在注目している円の半径，中心座標から，次の円の半径，中心座標を求める。

nextcircle(ri, xi, yi) = return ((-ri^3 + 3*ri^2*xi - 108*ri^2 + ri*xi^2 - 216*ri*xi + ri*yi^2 + 11664*ri - 3*xi^3 + 756*xi^2 - 3*xi*yi^2 - 58320*xi + 540*yi^2 - 2*sqrt(2)*yi*sqrt(-ri^4 - 108*ri^3 + 2*ri^2*xi^2 - 108*ri^2*xi + 2*ri^2*yi^2 + 11664*ri^2 + 108*ri*xi^2 - 23328*ri*xi + 108*ri*yi^2 + 1259712*ri - xi^4 + 108*xi^3 - 2*xi^2*yi^2 + 11664*xi^2 + 108*xi*yi^2 - 1259712*xi - yi^4 + 11664*yi^2) + 1259712)/(2*(ri^2 - 6*ri*xi + 216*ri + 9*xi^2 - 648*xi + 8*yi^2 + 11664)), 3*(ri^3 - 3*ri^2*xi + 180*ri^2 - ri*xi^2 - 216*ri*xi - ri*yi^2 + 3888*ri + 3*xi^3 - 108*xi^2 + 3*xi*yi^2 + 11664*xi + 36*yi^2 + 2*sqrt(2)*yi*sqrt(-ri^4 - 108*ri^3 + 2*ri^2*xi^2 - 108*ri^2*xi + 2*ri^2*yi^2 + 11664*ri^2 + 108*ri*xi^2 - 23328*ri*xi + 108*ri*yi^2 + 1259712*ri - xi^4 + 108*xi^3 - 2*xi^2*yi^2 + 11664*xi^2 + 108*xi*yi^2 - 1259712*xi - yi^4 + 11664*yi^2) - 419904)/(2*(ri^2 - 6*ri*xi + 216*ri + 9*xi^2 - 648*xi + 8*yi^2 + 11664)), -4*yi*(ri^2 + 54*ri - xi^2 + 54*xi - yi^2 - 5832)/(ri^2 - 6*ri*xi + 216*ri + 9*xi^2 - 648*xi + 8*yi^2 + 11664) + sqrt(2)*sqrt(-(ri^2 - 108*ri - xi^2 + 108*xi - yi^2)*(ri^2 + 216*ri - xi^2 - yi^2 + 11664))*(ri - 3*xi + 108)/(ri^2 - 6*ri*xi + 216*ri + 9*xi^2 - 648*xi + 8*yi^2 + 11664));

人円を初期値として，甲，乙，丙...と求める。ここで，R を外円の半径 108 として順次求められる円の半径 r で割ったもの R/r が整数で，規則があることがわかる。すなわち，外円の半径と i 番目の円の半径の比は i^2+2i+3 である。

R = 108
(r, x, y) = (18, 54, 72)
println("i = 1;  r = $r;  R/r = $(R/r)")
for i = 2:10
   (r, x, y) = nextcircle(r, x, y)
   println("i = $i;  r = $r;  R/r = $(round(Int, R/r)), $(i^2+2i+3)")
end

   i = 1;  r = 18;  R/r = 6.0
   i = 2;  r = 9.818181818181818;  R/r = 11, 11
   i = 3;  r = 5.999999999999991;  R/r = 18, 18
   i = 4;  r = 4.0;  R/r = 27, 27
   i = 5;  r = 2.8421052631578947;  R/r = 38, 38
   i = 6;  r = 2.117647058823516;  R/r = 51, 51
   i = 7;  r = 1.6363636363636425;  R/r = 66, 66
   i = 8;  r = 1.3012048192770946;  R/r = 83, 83
   i = 9;  r = 1.058823529411749;  R/r = 102, 102
   i = 10;  r = 0.8780487804877857;  R/r = 123, 123

人円から始めて，乾，兌，離，震，巽，坎，艮，坤についても同様に行う。

using SymPy
@syms R, r1, r2, r3, x3, y3, r4, x4, y4

R = 216 // 2
r1 = R // 2
eq1 = r1^2 + (R - r2)^2 - (r1 + r2)^2  # 天円と地円が外接する
eq2 = (r1 - x3)^2 + y3^2 - (r1 + r3)^2  # 天円と人円が外接する
eq3 = x3^2 + (R - r2 - y3)^2 - (r2 + r3)^2  # 地円と人円が外接する
eq4 = x3^2 + y3^2 - (R - r3)^2;  # 人円が外円に内接する

eq5 = (x3 - x4)^2 + (y4 - y3)^2 - (r3 + r4)^2  # 円と次の円が外接する
eq6 = x4^2 + y4^2 - (R - r4)^2                 # 次の円が外円に内接する
eq7 = x4^2 + (y4 - R + r2)^2 - (r2 + r4)^2     # 地円と次の円が外接する

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r2, r3, x3, y3, r4, x4, y4))

   4-element Vector{NTuple{7, Sym}}:
    (36, 18, 54, 72, 54/7, 270/7, 648/7)
    (36, 18, 54, 72, 54, 54, 0)
    (36, 54, -54, 0, 18, -54, 72)
    (36, 54, -54, 0, 54, 54, 0)

@syms r3, x3, y3, r4, x4, x5
r2 = 36
eq5 = (x3 - x4)^2 + (y4 - y3)^2 - (r3 + r4)^2  # 円と次の円が外接する
eq6 = x4^2 + y4^2 - (R - r4)^2                 # 次の円が外円に内接する
eq7 = x4^2 + (y4 - R + r2)^2 - (r2 + r4)^2     # 地円と次の円が外接する
solve([eq5, eq6, eq7], (r4, x4, y4))

nextcircle2(r3, x3, y3) = return ((-r3^3 + 2*r3^2*y3 - 108*r3^2 + r3*x3^2 + r3*y3^2 - 216*r3*y3 + 11664*r3 - 2*x3^2*y3 + 324*x3^2 - sqrt(3)*x3*sqrt(-r3^4 - 144*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 144*r3^2*y3 + 7776*r3^2 + 144*r3*x3^2 + 144*r3*y3^2 - 31104*r3*y3 + 1679616*r3 - x3^4 - 2*x3^2*y3^2 + 144*x3^2*y3 + 7776*x3^2 - y3^4 + 144*y3^3 + 7776*y3^2 - 1679616*y3 + 45349632) - 2*y3^3 + 540*y3^2 - 46656*y3 + 1259712)/(2*(r3^2 - 4*r3*y3 + 216*r3 + 3*x3^2 + 4*y3^2 - 432*y3 + 11664)), (-3*r3^2*x3 - 216*r3*x3 + sqrt(3)*r3*sqrt(-r3^4 - 144*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 144*r3^2*y3 + 7776*r3^2 + 144*r3*x3^2 + 144*r3*y3^2 - 31104*r3*y3 + 1679616*r3 - x3^4 - 2*x3^2*y3^2 + 144*x3^2*y3 + 7776*x3^2 - y3^4 + 144*y3^3 + 7776*y3^2 - 1679616*y3 + 45349632) + 3*x3^3 + 3*x3*y3^2 - 216*x3*y3 + 11664*x3 - 2*sqrt(3)*y3*sqrt(-r3^4 - 144*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 144*r3^2*y3 + 7776*r3^2 + 144*r3*x3^2 + 144*r3*y3^2 - 31104*r3*y3 + 1679616*r3 - x3^4 - 2*x3^2*y3^2 + 144*x3^2*y3 + 7776*x3^2 - y3^4 + 144*y3^3 + 7776*y3^2 - 1679616*y3 + 45349632) + 108*sqrt(3)*sqrt(-r3^4 - 144*r3^3 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - 144*r3^2*y3 + 7776*r3^2 + 144*r3*x3^2 + 144*r3*y3^2 - 31104*r3*y3 + 1679616*r3 - x3^4 - 2*x3^2*y3^2 + 144*x3^2*y3 + 7776*x3^2 - y3^4 + 144*y3^3 + 7776*y3^2 - 1679616*y3 + 45349632))/(2*(r3^2 - 4*r3*y3 + 216*r3 + 3*x3^2 + 4*y3^2 - 432*y3 + 11664)), sqrt(3)*x3*sqrt(-(-r3^2 - 216*r3 + x3^2 + y3^2 - 11664)*(-r3^2 + 72*r3 + x3^2 + y3^2 - 144*y3 + 3888))/(r3^2 - 4*r3*y3 + 216*r3 + 3*x3^2 + 4*y3^2 - 432*y3 + 11664) + (r3^3 - 2*r3^2*y3 + 216*r3^2 - r3*x3^2 - r3*y3^2 - 216*r3*y3 + 11664*r3 + 2*x3^2*y3 + 2*y3^3 - 108*y3^2)/(r3^2 - 4*r3*y3 + 216*r3 + 3*x3^2 + 4*y3^2 - 432*y3 + 11664));

R = 108
r2 = 36
(r, x, y) = (54/7, 270/7, 648/7)
for i = 1:10
   println("$i, r = $r, x = $x, y = $y, $(R/r), $(R/(2i^2 + 6i + 6))")
   (r, x, y) = nextcircle2(r, x, y)
end

   1, r = 7.714285714285714, x = 38.57142857142857, y = 92.57142857142857, 14.0, 7.714285714285714
   2, r = 4.153846153846184, x = 29.076923076923112, y = 99.69230769230774, 25.999999999999808, 4.153846153846154
   3, r = 2.571428571428608, x = 23.142857142857242, y = 102.85714285714286, 41.9999999999994, 2.5714285714285716
   4, r = 1.741935483871028, x = 19.16129032258083, y = 104.51612903225802, 61.99999999999786, 1.7419354838709677
   5, r = 1.2558139534883404, x = 16.32558139534871, y = 105.48837209302322, 86.00000000000217, 1.255813953488372
   6, r = 0.9473684210525465, x = 14.210526315789048, y = 106.10526315789478, 114.00000000001025, 0.9473684210526315
   7, r = 0.7397260273971918, x = 12.575342465752794, y = 106.52054794520554, 146.0000000000135, 0.7397260273972602
   8, r = 0.5934065934064815, x = 11.274725274724002, y = 106.81318681318686, 182.00000000003433, 0.5934065934065934
   9, r = 0.48648648648640563, x = 10.216216216215297, y = 107.0270270270272, 222.0000000000369, 0.4864864864864865
   10, r = 0.40601503759393454, x = 9.338345864660884, y = 107.18796992481218, 266.000000000033, 0.40601503759398494

R/r の数列は 2i^2 + 6i + 6 である。したがって，ri = R/(2i^2 + 6i + 6) である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 216 / 2
   r1 = R / 2
   (r2, r3, x3, y3) = (36, 18, 54, 72)
   println("地円の半径 = $r2, 人円の半径 = $r3")
   plot()
   circle(0, 0, R, :black)
   circle(r1, 0, r1, :red)
   circle(-r1, 0, r1, :red)
   circle(0, R - r2, r2, :blue)
   circle(0, -R + r2, r2, :blue)
   name1 = ["甲", "乙", "丙", "丁", "戊", "己", "庚", "辛", "壬", "癸"]
   (r, x, y) = (r3, x3, y3)
   for i = 1:10
       println("$(name1[i])円の半径　$r, 中心座標 $x, $y")
       circle(x, y, r, i)
       circle(x, -y, r, i)
       circle(-x, y, r, i)
       circle(-x, -y, r, i)
       (r, x, y) = nextcircle(r, x, y)
   end
   name2 = ["乾", "兌", "離", "震", "巽", "坎", "艮", "坤"]
   (r, x, y) = (54/7, 270/7, 648/7)
   for i = 1:8
       println("$(name2[i])円の半径　$r, 中心座標 $x, $y")
       circle(x, y, r, i)
       circle(x, -y, r, i)
       circle(-x, y, r, i)
       circle(-x, -y, r, i)
       (r, x, y) = nextcircle2(r, x, y)
   end
   if more == true
       point(0, 0, " O")
       point(r1, 0, " r1")
       point(R, 0, " R")
       point(0, R-r2, " R-r2")
       point(x3, y3, "(X3,y3,r3)")
       point(50, 25, "天円", mark=false)
       point(10, 85, "地円", mark=false)
       point(50, 85, "人円", mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

draw(false)

   地円の半径 = 36, 人円の半径 = 18
   甲円の半径　18, 中心座標 54, 72
   乙円の半径　9.818181818181818, 中心座標 78.54545454545455, 58.90909090909091
   丙円の半径　5.999999999999991, 中心座標 90.0, 48.0
   丁円の半径　4.0, 中心座標 96.0, 40.0
   戊円の半径　2.8421052631578947, 中心座標 99.47368421052632, 34.10526315789474
   己円の半径　2.117647058823516, 中心座標 101.64705882352943, 29.64705882352941
   庚円の半径　1.6363636363636425, 中心座標 103.09090909090908, 26.181818181818187
   辛円の半径　1.3012048192770946, 中心座標 104.09638554216873, 23.421686746987955
   壬円の半径　1.058823529411749, 中心座標 104.82352941176475, 21.176470588235293
   癸円の半径　0.8780487804877857, 中心座標 105.36585365853662, 19.317073170731703
   乾円の半径　7.714285714285714, 中心座標 38.57142857142857, 92.57142857142857
   兌円の半径　4.153846153846184, 中心座標 29.076923076923112, 99.69230769230774
   離円の半径　2.571428571428608, 中心座標 23.142857142857242, 102.85714285714286
   震円の半径　1.741935483871028, 中心座標 19.16129032258083, 104.51612903225802
   巽円の半径　1.2558139534883404, 中心座標 16.32558139534871, 105.48837209302322
   坎円の半径　0.9473684210525465, 中心座標 14.210526315789048, 106.10526315789478
   艮円の半径　0.7397260273971918, 中心座標 12.575342465752794, 106.52054794520554
   坤円の半径　0.5934065934064815, 中心座標 11.274725274724002, 106.81318681318686

算額（その178）

算額（その178）

愛媛県四国中央市川之江町 八幡神社
川之江高等学校の和算勉強会が奉納した「算額(和算の術)」平成21年
https://www.fureai-cloud.jp/tie/doc/view/10456/

正方形内の一つの対角線上に中心がある三個の等円を図のように描く。等円はそれぞれ外接し，そのうち二つは正方形に内接している。このとき，正方形の一辺の長さが一尺であるとき，等円の直径はいくらか。

和讃家の山口坎山が明和5年（1768）に『道中日記』で 15 通りの解法を著したということのようであるが，最もありふれた解法であろうけれども以下を。

正方形の一辺の長さを 2，等円の半径を r として，図のように記号を定め方程式を解く。
2次方程式なので中学生でも解けるのかな。

using SymPy
@syms r::positive;

eq = 2(1 - r)^2 - 4r^2
res = solve(eq)[1]
println("$res = $(res.evalf())")

   -1 + sqrt(2) = 0.414213562373095

元の単位においては，等円の直径は 4寸1分4厘2毛1糸あまりあり

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = √2 - 1
   println("r = $r")
   plot([1, 1, -1, -1, 1], [-1, 1, 1, -1, -1], color=:black, lw=0.5)
   circle(0, 0, r, :yellow, fill=true)
   circle(1 - r, 1 - r, r, :red, fill=true)
   circle(r - 1, r - 1, r, :blue, fill=true)
   if more == true
       point(1 - r, 1 - r, " (1-r,1-r)", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その177）

算額（その177）

愛媛県四国中央市川之江町 八幡神社
川之江高等学校の和算勉強会が奉納した「算額(大円内等円発展題)」平成21年
https://www.fureai-cloud.jp/tie/doc/view/10455/

径が 1 寸の外円の中に，等円を 3 個入れる場合，等円の径を求めよ。5 個，8 個入れる場合も求めよ。

外円の半径を 1 とし，その中に半径 r の n 個の等円を入れることを考える(n ≧ 2)。
時計盤でいえば12時の位置にある等円と，それに対して反時計回りで左隣にある等円を考える。
隣り合う円が外接するという方程式を解き等円の半径を求めることができる。

以下は，n を与えて r（数式と数値）を求める関数である。

using SymPy
@syms r::positive, x, y

function getradius(n=3)
   θ = 2PI / n + PI / 2
   x = (1 - r)cos(θ)
   y = (1 - r)sin(θ)
   eq = x^2 + (1 - r - y)^2 - 4r^2  # 隣り合う等円が外接する
   res = solve(eq)[1]
   simplify(res), res.evalf()
end;

実際に求めてみる。
等円を 3 個入れる場合

getradius(3)

   (-3 + 2*sqrt(3), 0.464101615137755)

術曰　3の平方根を2倍し3を引く
答曰　4分6厘4毛1糸あまりあり

等円を 5 個入れる場合

getradius(5)

   (-5 - sqrt(50 - 10*sqrt(5))/2 + 3*sqrt(10 - 2*sqrt(5))/2 + 2*sqrt(5), 0.370191908158750)

等円を 8 個入れる場合

getradius(8)

   (-3 - sqrt(4 - 2*sqrt(2)) + 2*sqrt(2 - sqrt(2)) + 2*sqrt(2), 0.276768653914155)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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(n, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   expression, r = getradius(n)
   println("r = $expression\n  = $r")
   plot()
   circle(0, 0, 1, :black)
   for i = 0:n-1
       θ = 2pi/n*i + pi/2
       (x, y) = (1 - r) .* (cos(θ), sin(θ))
       circle(x, y, r, i, fill=true)
       if more == true && i == 1
           point(0, 1 - r, " 1-r", :black)
           point(x, y, "(x,y)", :black)
       end
   end
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

draw(9, false)

   r = (1 - sin(pi/9))*sin(pi/9)/cos(pi/9)^2
     = 0.254854701717149

算額（その176）

算額（その176）

岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第5問: 囲円を 6 個の等円で囲む。楕円の短径の端は囲円に接する。楕円の長径と短径を知って囲円の径を求めよ。

囲円の半径を r，楕円の長径，短径を a, b と置く。
⊿OAB は正三角形である。

楕円の接線 https://blog.goo.ne.jp/r-de-r/e/b19211e3045786a3a9108d7dc4f41b6d を参照。

図で，r + b = sqrt(3a^2 + b^2) である。つまり，r = sqrt(3a^2 + b^2) - b である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function ellipse(ox, oy, ra, rb; φ=0, beginangle=0, endangle=360,
                    color=:black, lty=:solid, lwd=0.5, fcolor="")
"""
(ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。
fcolor を指定すると塗りつぶし。
"""
   θ = beginangle:0.1:endangle
   if φ == 0
       if fcolor == ""
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   else
       x = ra .* cosd.(θ)
       y = rb .* sind.(θ)
       cosine = cosd(φ)
       sine = sind(φ)
       if fcolor == ""
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   end
end;

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = sqrt(3a^2 + b^2) - b
   plot()
   circle(0, 0, r, :aquamarine3, fill=true)
   for i = 120:60:420
       ellipse((r + b)cosd(i - 90), (r + b)sind(i - 90), a, b, φ = i)
   end
   if more == true
       point(0, r, " r")
       point(0, r + b, " r+b")
       point(0, 0, " O")
       segment(-7, r + b, 7, r + b)
       segment(0, 0, 12cosd(60), 12sind(60))
       segment(0, 0, -12cosd(60), 12sind(60))
       point((r + b)/sqrt(3), (r + b), " A")
       point(-(r + b)/sqrt(3), (r + b), "B ", :green, :right)
       point((r + b)cosd(30), (r + b)sind(30), "(x0,y0)")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

算額（その175）

算額（その175）

岐阜県不破郡垂井町宮代 南宮大社奉納算額 天保13年
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第1問: 楕円内に赤円 2 個と最も多くの青等円を奇数個入れる。赤円 2 個と青円 1 個の直径の輪は短径に等しい。また，楕円の左右の端で青円は楕円に接する。赤円と青円の径を知って，青等円の個数を求めよ。

赤円，青円の半径をそれぞれ r1, r2 とおく。以下の方程式を立てる。

using SymPy

@syms n::integer, r1::positive, r2::positive, x::positive, y::positive;
eq1 = x^2/(n*r2)^2 + y^2/(r2 + 2r1)^2 - 1
eq2 = (x -(n - 1)r2)^2 + y^2 - r2^2;

eq2 - eq1*(r2 + 2r1)^2 が (x - n*r2) の因子を持つので

eq3 = eq2 - eq1*(r2 + 2r1)^2 |> factor
eq3 |> println

   -(-n*r2 + x)*(n^3*r2^3 - 2*n^2*r2^3 - n^2*r2^2*x + 4*n*r1^2*r2 + 4*n*r1*r2^2 + n*r2^3 + 4*r1^2*x + 4*r1*r2*x + r2^2*x)/(n^2*r2^2)

除算して簡約化する。

eq4 = eq3  / (x - n*r2) |> simplify
eq4 |> println

   -n*r2 + 2*r2 + x - 4*r1^2/(n*r2) - 4*r1/n - r2/n - 4*r1^2*x/(n^2*r2^2) - 4*r1*x/(n^2*r2) - x/n^2

x = n*r2 を代入して n について解くと解が得られる。

solve(eq4(x => n*r2), n)[1] |> println

   (2*r1 + r2)^2/r2^2

なお，演算順序を変えると解が求まらないことがある。

また，r1, r2 にでたらめな数を与えると，題意に沿わない図になる場合がある。

r1 = 3, r2 = 9 のとき，青円は偶数個になる。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
end;

function ellipse(ox, oy, ra, rb; φ=0, beginangle=0, endangle=360,
                    color=:black, lty=:solid, lwd=0.5, fcolor="")
"""
(ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。
fcolor を指定すると塗りつぶし。
"""
   θ = beginangle:0.1:endangle
   if φ == 0
       if fcolor == ""
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   else
       x = ra .* cosd.(θ)
       y = rb .* sind.(θ)
       cosine = cosd(φ)
       sine = sind(φ)
       if fcolor == ""
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   end
end;

function draw(R1, R2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (R1, R2) ./ 2
   n = floor(Int, (r2 + 2r1)^2/r2^2)
   plot()
   ellipse(0, 0, n*r2, r2 + 2r1)
   circle(0, r2 + r1, r1, :indianred1,fill=true)
   circle(0, -r2 - r1, r1, :indianred1,fill=true)
   for i = -(n - 1)//2:(n - 1)//2
       circle(2i*r2, 0, r2, :slateblue1, fill=true)
   end
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

楕円の接線

楕円の接線

原点を中心とする長径 a，短径 b の楕円がある。y 軸上に中心を持つ円がこの楕円に外接し，さらに円の中心を通る楕円の接線の傾きが slope であるとき，円の半径はいかほどか。

using SymPy

@syms a, b, x1, y1, r, tanθ

eq1 = x1^2/a^2 + y1^2/b^2 - 1
eq2 = -b^2*x1 / (a^2*y1) + (r + b - y1)/x1  #
eq3 = (r + b - y1)/x1 + tanθ;

res = solve([eq1, eq2, eq3], (x1, y1, r))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-a^2*tanθ/sqrt(a^2*tanθ^2 + b^2), b^2/sqrt(a^2*tanθ^2 + b^2), -b + sqrt(a^2*tanθ^2 + b^2))
    (a^2*tanθ/sqrt(a^2*tanθ^2 + b^2), -b^2/sqrt(a^2*tanθ^2 + b^2), -b - sqrt(a^2*tanθ^2 + b^2))

res[1][3] |> println

   -b + sqrt(a^2*tanθ^2 + b^2)

以上の結果に基づき，以下の関数を定義する。

引数
　楕円の長径 a，短径 b，接線の傾き（符号に注意）
戻り値
　接点の座標 (x, y)，楕円に外接する円の半径 r

func(a, b, tanθ) = (-a^2*tanθ/sqrt(a^2*tanθ^2 + b^2), b^2/sqrt(a^2*tanθ^2 + b^2), -b + sqrt(a^2*tanθ^2 + b^2));

使用例

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function ellipse(ox, oy, ra, rb; φ=0, beginangle=0, endangle=360,
                    color=:black, lty=:solid, lwd=0.5, fcolor="")
"""
(ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。
fcolor を指定すると塗りつぶし。
"""
   θ = beginangle:0.1:endangle
   if φ == 0
       if fcolor == ""
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   else
       x = ra .* cosd.(θ)
       y = rb .* sind.(θ)
       cosine = cosd(φ)
       sine = sind(φ)
       if fcolor == ""
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd)
       else
           plot!(cosine .* x .- sine .* y .+ ox,
                 sine .* x .+ cosine .* y .+ oy,
                 linecolor=color, linestyle=lty, linewidth=lwd,
                 seriestype=:shape, fillcolor=fcolor)
       end
   end
end;

using Plots
using Printf

function draw(a, b, c, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, y1, r1) = func(a, b, c)
   println("x1 = $x1, y1 = $y1, r1 = $r1")
   println("slope = $(-(r1 + b - y1)/x1)")
   plot()
   circle(0, r1+b, r1)
   ellipse(0, 0, a, b)
   segment(0, r1+b, x1, y1)
   if more == true
       point(x1, y1, " (x1,y1)", :green, :left, :bottom)
       point(0, b, " b")
       point(a, 0, " a", :green, :left, :top)
       point(0, r1+b, " r1+b")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

function draw3(a, b, num=6)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, y, r) = func(a, b, -cot(pi/num))
   plot()
   circle(0, 0, r, :aquamarine3, fill=true)
   deg = 360/num
   for i = 1:num
       deg2 = i*deg
       ellipse((r + b)cosd(deg2 - 90), (r + b)sind(deg2 - 90), a, b, φ = deg2, color=:lightgoldenrod1, fcolor=:lightgoldenrod1)
   end
   plot!(showaxis=false)
end;

draw3(3, 2, 11)

算額（その174）

算額（その174）

岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第7問: 5 個の正方形，赤，黃，青，黒，白の一辺の長さの差は等しい。赤と黃の面積の和と青と黒と白の面積の和を知って，白の面積を求めよ。
白の一辺の長さを w，それぞれの一辺の差を d とすると，
赤の正方形の面積 = (w + 4d)^2
黃の正方形の面積 = (w + 3d)^2
青の正方形の面積 = (w + 2d)^2
黒の正方形の面積 = (w + d)^2
白の正方形の面積 = w^2
である。
赤と黃の面積の和を RY
青と黒と白の面積の和を BKW
とし，方程式を立てる。

以下は w = d = 1 の場合

using SymPy

@syms w::positive, d::positive, RY::positive, BKW::positive;

eq1 = (w + 4d)^2 + (w + 3d)^2 - RY |> expand
eq2 = (w +2d)^2 + (w + d)^2 + w^2 - BKW|> expand;

連立方程式を BKW, RY について解く。

res = solve([eq1, eq2], (w, d))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-42*BKW + 18*RY - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))*sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365)/(6*(2*BKW - 3*RY)), sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365))
    ((-42*BKW + 18*RY + 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))*sqrt(16*BKW/365 + 21*RY/365 + 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365)/(6*(2*BKW - 3*RY)), sqrt(16*BKW/365 + 21*RY/365 + 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365))

2 組の解が得られるが，BKW = 20, RY = 50 のときを計算してみると

res[1][1](BKW => 20, RY => 50).evalf(), res[1][2](BKW => 20, RY => 50).evalf()

   (1.43525126802703, 1.01117744530322)

res[2][1](BKW => 20, RY => 50).evalf(), res[2][2](BKW => 20, RY => 50).evalf()

   (-4.07737480549386, 2.54644251637035)

2 番目の解では w が負の値になるので不適切解である。

f(RY, BKW) = ((-42*BKW + 18*RY - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))*sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365)/(6*(2*BKW - 3*RY)), sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365));

(w, d) = f(41, 14)
println("w = $w,  d = $d")
println("RY  = $((w + 4d)^2 + (w + 3d)^2)")
println("BKW = $((w +2d)^2 + (w + d)^2 + w^2)")

   w = 1.0000000000000002,  d = 1.0000000000000002
   RY  = 41.000000000000014
   BKW = 14.000000000000007

白の正方形の面積は res[1][1]^2 で，簡約化して以下を得る。

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

   (213*BKW - 17*RY - 16*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))/365

RY = 41
BKW = 14
g(RY, BKW) = (213*BKW - 17*RY - 16*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))/365  # 白の面積
g(41, 14) |> println
g(50, 20) |> println

   1.0
   2.059946202373195

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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 rect(x1, y1, x2, y2, color=:pink)
   plot!([x1, x2, x2, x1, x1], [y1, y1, y2,  y2, y1], color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   colors = [:indianred1, :yellow, :steelblue1, :gray30, :white]
   plot()
   for i = 5:-1:1
       width = w + (i-1)d
       rect(0, 0, width, width, colors[6-i])
   end
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

算額（その173）

算額（その173）

岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第6問: 外円内にこの円と同じ半径の円弧を上下左右に作り，その隙間に 5 個の円を入れる。外円の径を知って個々の累円の径を求めよ。

外円の径を R とし，累円の径を r1, r2, r3, r4, r5 とおき，連立方程式を解く。

using SymPy

@syms R::positive, r1::positive, r2::positive, r3::positive, r4::positive, r5::positive;

eq1 = R^2 + (R - r1)^2 - (R + r1)^2
eq2 = R^2 + (R - 2r1 - r2)^2 - (R + r2)^2
eq3 = R^2 + (R - 2r1 - 2r2 - r3)^2 - (R + r3)^2
eq4 = R^2 + (R - 2r1 - 2r2 - 2r3 - r4)^2 - (R + r4)^2
eq5 = R^2 + (R - 2r1 - 2r2 - 2r3 - 2r4 - r5)^2 - (R + r5)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r3, r4, r5))

   1-element Vector{NTuple{5, Sym}}:
    (R/4, R/12, R/24, R/40, R/60)

[4, 12, 24, 40, 60, ...] は 4[1, 3, 6, 10, 15, ...] で，これは一般項 4(n*(n+1)/2) = 2n*(n+1) であることはよく知られていると思うが，これを SymPy でやるとどうなるかをいかに示しておく。

一般項が何次式になるかわからない場合には次数を多めに取り，一般項を f として定義する。今回の場合は 3次式と仮定する。ついで，使用した次数分の項を terms ベクトルに用意する。f(n => i) - terms[i] を i = 1,2,3,4 とし，連立方程式を解く。

まずはまだるっこしい方法で，記述が面倒くさい。

@syms a, b, c, d, n
f = a*n^3 + b*n^2 + c*n + d
terms = [4, 12, 24, 40]
coeffs = solve([f(n => 1) - terms[1], f(n => 2) - terms[2], f(n => 3) - terms[3], f(n => 4) - terms[4]], (a, b, c, d)) 

   Dict{Any, Any} with 4 entries:
     c => 2
     d => 0
     a => 0
     b => 2

f(coeffs) |> simplify |> println

   2*n*(n + 1)

内包表記を使えば重複する記述を避けることができる。

x = [f(n => i).evalf() - terms[i] for i in 1:4];
x |> println

   Sym[a + b + c + d - 4, 8.0*a + 4.0*b + 2.0*c + d - 12, 27.0*a + 9.0*b + 3.0*c + d - 24, 64.0*a + 16.0*b + 4.0*c + d - 40]

f(solve(x)) |> simplify |> println

   2.0*n*(n + 1)

外円の直径を R とすると n 番目の累円の直径は 1/(2n*(n+1))R/4 = R / (2n*(n+1))

[1/(2n*(n+1)) for n in 1:5]

   5-element Vector{Float64}:
    0.25
    0.08333333333333333
    0.041666666666666664
    0.025
    0.016666666666666666

一方，隣り合う累円について漸化式を立てれば以下のようになる。

@syms R, ri, rip1
a = sqrt(2R*ri + ri^2) - sqrt(2R*rip1 + rip1^2) - ri - rip1
a |> println

   -ri - rip1 + sqrt(2*R*ri + ri^2) - sqrt(2*R*rip1 + rip1^2)

これを rip1 についてとくと，今の累円の半径 ri の次の累円の半径 rip1 を求める関数を定義できる。

solve(a, rip1)[1] |> println

   ri*(R + ri - sqrt(ri*(2*R + ri)))/(R - ri + sqrt(ri*(2*R + ri)))

@syms R, r
f = r*(R + r - sqrt(r*(2R + r)))/(R - r + sqrt(r*(2R + r)))
f |> println

   r*(R + r - sqrt(r*(2*R + r)))/(R - r + sqrt(r*(2*R + r)))

これを Julia の関数定義で書くと以下のようになる。

g(r, R) = r*(R + r - sqrt(r*(2R + r)))/(R - r + sqrt(r*(2R + r)));

外円の半径を R とすると，青円の半径は 1/4 なので，次の赤円の半径は以下のように 0.08333333333333333 となる。円の中心座標も計算できる。

g(1//4, 1)

   0.08333333333333333

関数 g(r, R) を使って，順次累円を描くことができる。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw2(i, r, x, R)
   circle(R*cosd(90+60i), R*sind(90+60i), R, :black, beginangle=270+60i, endangle=330+60i)
   circle(R*cosd(-30+60i), R*sind(-30+60i), R, :black, beginangle=90+60i, endangle=90+60+60i)
   for j = 1:5
       #println("x = $(x[j]*cosd(60i)), y = $(x[j]*sind(60i)), r = $(r[j])")
       circle(x[j]*cosd(60i), x[j]*sind(60i), r[j], j, fill=true)
   end
end

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = 5
   r = zeros(n)
   x = zeros(n)
   R = 1
   r[1] = R / 4
   x[1] = R - r[1]
   plot()
   circle(0, 0, R, :black)
   for i = 2:5
       r[i] = g(r[i-1], R)
       x[i] = x[i-1] - r[i-1] - r[i]
   end
   for i in 0:5
       draw2(i, r, x, R)
   end
   if more == true
       point(R, 0, " R", :black)
       point(R - r[1], 0, " R-r1", :black, :left, :bottom)
       point(R - 2r[1] - r[2], 0, " R-2r1-r2", :black)
       segment(0, R, R - r[1], 0, :red, linestyle=:dash, linewidth=0.5)
       segment(0, R, R - 2r[1] - r[2], 0, :red, linestyle=:dash, linewidth=0.5)
       point(0, R, " R", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

算額（その172）

算額（その172）

岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第4問: 囲円を 8 個の等円で囲む。等円 8 個は互いに外接し，囲円に接している。囲円の径を知って等円の径を求めよ。

この算額は「算額（その60）」の簡単版である。https://blog.goo.ne.jp/r-de-r/e/06710f85cfba141184ddd8f21529fae0

囲円の径 r1，等円の径 r2，右上の等円の中心座標を (x, r2) とする。図のように記号を定め方程式を解く。方程式は 2 本なので，x, r2 を r1 で表す解を求める。

using SymPy

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

eq1 = x^2 + r2^2 - (r1 + r2)^2  # 等円と囲円が外接する
eq2 = 2(x- r2)^2 - 4r2^2;  # 等円同士が外接する
res = solve([eq1, eq2], (x, r2))  # x, r2 は r1 を含む式で得られる

   3-element Vector{Tuple{Sym, Sym}}:
    (-(-12*r1^3*(-2*sqrt(2) - sqrt(20 - 14*sqrt(2)) + 3)^2 + r1^3*(-2*sqrt(2) - sqrt(20 - 14*sqrt(2)) + 3)^3 - r1^3*(-2*sqrt(2) - sqrt(20 - 14*sqrt(2)) + 3) + 2*r1^3)/(2*r1^2), r1*(-2*sqrt(2) - sqrt(20 - 14*sqrt(2)) + 3))
    (-(-12*r1^3*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3)^2 - r1^3*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3) + r1^3*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3)^3 + 2*r1^3)/(2*r1^2), r1*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3))
    (-(-12*r1^3*(2*sqrt(2) + 3 + sqrt(14*sqrt(2) + 20))^2 - r1^3*(2*sqrt(2) + 3 + sqrt(14*sqrt(2) + 20)) + 2*r1^3 + r1^3*(2*sqrt(2) + 3 + sqrt(14*sqrt(2) + 20))^3)/(2*r1^2), r1*(2*sqrt(2) + 3 + sqrt(14*sqrt(2) + 20)))

r2 は正で，r1 より小さくなければならないので，2組目のものが適解である。

res[1][2](r1 => 1).evalf() |> println

   -0.276768653914155

res[2][2](r1 => 1).evalf() |> println

   0.619914404421775

res[3][2](r1 => 1).evalf() |>  println

   12.1370711845441

3 -2√2 + sqrt(20 - 14√2)

   0.6199144044217748

等円の径は囲円の径の 3 -2√2 + sqrt(20 - 14√2) =　0.619914404421775 倍である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (x, r2) = (-(-12*r1^3*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3)^2 - r1^3*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3) + r1^3*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3)^3 + 2*r1^3)/(2*r1^2), r1*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3))
   @printf("囲円（赤）の径は等円（青）の径の 3 -2√2 + sqrt(20 - 14√2) = %.6f 倍である。", r2/r1)
   plot()
   circle(0, 0, r1, :indianred1, fill=true)
   circle(x, r2, r2, :steelblue1, fill=true)
   circle(-x, r2, r2, :steelblue1, fill=true)
   circle(r2, x, r2, :steelblue1, fill=true)
   circle(-r2, x, r2, :steelblue1, fill=true)
   circle(x, -r2, r2, :steelblue1, fill=true)
   circle(-x, -r2, r2, :steelblue1, fill=true)
   circle(r2, -x, r2, :steelblue1, fill=true)
   circle(-r2, -x, r2, :steelblue1, fill=true)
   if more == true
       point(x, r2, "(x,r2)", :black)
       point(r2, x, "(r2,x)", :black)
       point(r1, 0, " r1", :black, :left, :bottom)
       point(0, 0, " 0", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

   囲円（赤）の径は等円（青）の径の 3 -2√2 + sqrt(20 - 14√2) = 0.619914 倍である。

 

算額（その171）

算額（その171）

岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第3問: 半円内に直角二等辺三角形を作り，その隙間に青，赤の 6 個の円を入れる，半円の直径を知って，赤円の径を求めよ。

半円の径を R とする。後々 R は 8 の倍数にしておくほうが方程式の係数が簡単になことがわかる。

白円，青円，赤円の半径を r0, r1,  r2 とし，赤円の中心座標を (x2, y2) とする。

r0, r1 はすぐに計算できる。r2, x2, y2 が未知数であるが，solve() では以下の 4 個の方程式のどの 3 つをとっても，有限の時間内には解が見つからないようである。
そこでいつものように nlsolve で解を求める。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

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

R = 16
r0 = R/(1 + sqrt(Sym(2)))
r1 = R*(2 - sqrt(Sym(2)))/4
x1 = (R - r1) / sqrt(Sym(2))
r2 = R // 8
eq1 = distance(R, 0, 0, R, x2, y2) - r2^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = (x2 - x1)^2 + (y2 - x1)^2 - (r1 + r2)^2
eq4 = x2^2 + (y2 - r0) - (r0 + r2)^2;

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")

   2*(x2/2 + y2/2 - 8)^2 - 4,
   x2^2 + y2^2 - 196,
   (x2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 + (y2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 - (10 - 4*sqrt(2))^2,
   x2^2 + y2 - (2 + 16/(1 + sqrt(2)))^2 - 16/(1 + sqrt(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=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)
   (r2, x2, y2) = u
   return [
       2*(x2/2 + y2/2 - 8)^2 - 4,
       x2^2 + y2^2 - 196,
       (x2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 + (y2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 - (10 - 4*sqrt(2))^2,
       # x2^2 + y2 - (2 + 16/(1 + sqrt(2)))^2 - 16/(1 + sqrt(2)),
   ]
end;

iniv = [2.0, 6, 12]
res = nls(H, ini=iniv)
println(res)

   ([2.0, 6.352746103452378, 12.475681021293815], true)

経験則で r2 は R の 1/8 となる。
そこで r2 = R // 8 を条件に加えて eq1, eq2 の二元連立方程式を解く。

res = solve([eq2, eq1], (x2, y2))

   2-element Vector{Tuple{Sym, Sym}}:
    (-(sqrt(2) + sqrt(32 - 16*sqrt(2)) + 8)^3/34 - 30*sqrt(32 - 16*sqrt(2))/17 - 30*sqrt(2)/17 - 32/17 + 8*(sqrt(2) + sqrt(32 - 16*sqrt(2)) + 8)^2/17, sqrt(2) + sqrt(32 - 16*sqrt(2)) + 8)
    (-(-sqrt(32 - 16*sqrt(2)) + sqrt(2) + 8)^3/34 - 30*sqrt(2)/17 - 32/17 + 30*sqrt(32 - 16*sqrt(2))/17 + 8*(-sqrt(32 - 16*sqrt(2)) + sqrt(2) + 8)^2/17, -sqrt(32 - 16*sqrt(2)) + sqrt(2) + 8)

x2 が特に複雑な式になり，これが solve() で解が求まらない原因と思われる。

赤円の半径は 半円の径の 1/8 である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   R = 16
   r0 = R / (1 + √2)
   r1 = (R - R/√2)/2
   x1 = (R -r1)/√2
   r2 = R // 8
   (x2, y2) = res[1]
   (x22, y22) = res[2] 
   circle(0, 0, R, :yellow, fill=true, beginangle=0, endangle=180)
   plot!([R, 0, -R, R], [0, R, 0, 0], linecolor=:palegreen2, linewidth=0.5, seriestype=:shape, fillcolor=:palegreen2)
   circle(0, r0, r0, :white, fill=true)
   circle(x2, y2, r2, :indianred1, fill=true)
   circle(x22, y22, r2, :indianred1, fill=true)
   circle(x1, x1, r1, :steelblue1, fill=true)
   circle(-x2, y2, r2, :indianred1, fill=true)
   circle(-x22, y22, r2, :indianred1, fill=true)
   circle(-x1, x1, r1, :steelblue1, fill=true)
   if more == true
       point(x2, y2, "(x2,y2)", :black)
       point(x1, x1, "((R-r1)/√2,(R-r1)/√2)", :black, :center)
       point(0, r0, " r0", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その170）

算額（その170）

岐阜県大垣市西外側町 大垣八幡神社 天保年間
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第29問: 直行するに直線に接し互いに外接する 3 個の等円（赤）を作り，青円と黒円を置く。青円の直径を知って黒円の直径を求めよ。

赤円，青円，黒円の直径をそれぞれ r1, r2, r3 として，更に，青円，黒円の中心の y 座標を y2, y3 として連立方程式を解く。なお，2 個の青円の中心が x 軸になるように座標を定めないと，連立方程式の係数が複雑になり解けないので注意。

using SymPy

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

eq1 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = (2r1 - r3)^2 + (sqrt(Sym(3))r1 - y3)^2 - (r1 + r3)^2
eq4 = (r2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = (2r1 - r2)^2 + (y2 - sqrt(Sym(3))r1)^2 - (r1 + r2)^2;

res = solve([eq1, eq3, eq4, eq5], (r1, r3, y2, y3))

   1-element Vector{NTuple{4, Sym}}:
    (4*sqrt(2)*r2/9 + r2, -2*sqrt(2)*r2 + 11*r2/3, 4*sqrt(6)*r2/9 + 16*sqrt(3)*r2/9, -49*sqrt(3)*r2*(116/441 - 580*sqrt(2)/441)/58)

黒円の半径は -2*sqrt(2)*r2 + 11*r2/3 である。
つまり，「黒円の半径 = (11 - 6√2)/3 × 青円の半径」である。

-2*sqrt(Sym(2))*r2 + 11*r2/3 |> simplify |> println

   r2*(11 - 6*sqrt(2))/3

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1  # 青円の半径を 1 として,左右反転した図を描く
   (r1, r3, y2, y3) = (4*sqrt(2)*r2/9 + r2, -2*sqrt(2)*r2 + 11*r2/3, 4*sqrt(6)*r2/9 + 16*sqrt(3)*r2/9, -49*sqrt(3)*r2*(116/441 - 580*sqrt(2)/441)/58)
   plot()
   circle(r1, 0, r1, :indianred1, fill=true)
   circle(3r1, 0, r1, :indianred1, fill=true)
   circle(2r1, sqrt(3)r1, r1, :indianred1, fill=true)
   circle(r2, y2, r2, :steelblue1, fill=true)
   circle(r3, y3, r3, :snow4, fill=true)
   hline!([0], color=:black, lw=0.5)
   vline!([0], color=:black, lw=0.5)
   if more == true
       point(r2, y2, "(r2,y2)", :black)
       point(r3, y3, "(r3,y3)", :black)
       point(r1, 0, "r1", :black)
       point(2r1, sqrt(3)r1, "(2r1,√3r1)", :black)
   end
end;

算額（その169）

算額（その169）

岐阜県大垣市西外側町 大垣八幡神社 天保年間
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第28問: 互いに中心を通る 2 等円（黃）に 3 個の等円（青）を入れ，更に赤円 2 個，黒円 4 個を入れる。青円径を知って黒円径を求めよ。

算額（その85）を一段階複雑にしたものである。

青円，赤円，黒円の半径を r1, r2, r3，黒円の中心座標を (x3, y3) とする。黃円の半径は 2r1 である。未知数 5 個に対して方程式は 4 本なので，r2, r3, x3, y3 が r1 を含む解が求まる。

using SymPy

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

eq1 = x3^2 + y3^2 - (r1 + r3)^2  # 青円と黒円が外接
eq2 = x3^2 + (r1 + r2 - y3)^2 - (r2 + r3)^2  # 赤円と黒円が外接
eq3 = r1^2 + (r1 + r2)^2 - (2r1 - r2)^2  # 赤円が黃円に内接
eq4 = (r1 + x3)^2 + y3^2 - (2r1 - r3)^2  # 黒円が黃円に内接
res = solve([eq1, eq2, eq3, eq4], (r2, r3, x3, y3))

   1-element Vector{NTuple{4, Sym}}:
    (r1/3, 2*r1/11, 5*r1/11, 12*r1/11)

黒円の半径は青円の半径の 2/11 である。
赤円の半径は青円の半径の 1/3 である。
黒円の中心座標は青円の半径の (5/11, 12/11) である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1  # 青円の半径を 1 として図を描く
   (r2, r3, x3, y3) = (r1/3, 2*r1/11, 5*r1/11, 12*r1/11)
   plot()
   circle(r1, 0, 2r1, :yellow1, fill=true)
   circle(-r1, 0, 2r1, :yellow1, fill=true)
   circle(r1, 0, 2r1, :snow4)
   circle(-r1, 0, 2r1, :snow4)
   circle(0, 0, r1, :steelblue1, fill=true)
   circle(2r1, 0, r1, :steelblue1, fill=true)
   circle(-2r1, 0, r1, :steelblue1, fill=true)
   circle(0, r1 + r2, r2, :indianred1, fill=true)
   circle(0, -r1 - r2, r2, :indianred1, fill=true)
   circle(x3, y3, r3, :snow4, fill=true)
   circle(-x3, y3, r3, :snow4, fill=true)
   circle(x3, -y3, r3, :snow4, fill=true)
   circle(-x3, -y3, r3, :snow4, fill=true)
   if more == true
       point(0, r1, " r1", :black)
       point(0, r1+r2, " r1+r2", :black, :left, :bottom)
       point(x3, y3, "(x3,y3,r3)", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;