算額（その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)