算額（その1117）

算額（その1117）

四 岩手県花巻市南笹間 東光寺 慶応2年(1866)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：累円，外円

大円の中に小円と初円を容れる。大円，小円と外接する円を初円と呼ぶ。初円が決まれば次の円（二円），また次の円（三円）と云うように，順次円（累円）が決まる。大円，小円，初円の直径がそれぞれ 168 寸，88 寸，14 寸のとき，次々に円を決めていくとき，円の直径が 3 寸になるのは初円を 1 番目として数えると何番目の円か。

算額（その676）も参照のこと。

小円の中心が x 軸上に来るようにして図を描く。
大円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r1, (R - r1, 0)
初円の半径と中心座標を r2, (x2, y2)
二円の半径と中心座標を r3, (x3, y3)
とおき，まず初円の中心座標を決定する。

include("julia-source.txt")

using SymPy
@syms R, r1, r2, x2, y2
#(R, r1, r2) = (168, 88, 14) .// 2
eq1 = (R - r1 - x2)^2 + y2^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
res2 = solve([eq1, eq2], (x2, y2))[2]

   ((-R^2 + R*r1 + R*r2 + r1*r2)/(-R + r1), 2*sqrt(-R*r1*r2*(-R + r1 + r2))/(-R + r1))

二円以降は漸化式のようにして，前の円の情報を元にして次の円の半径，中心座標を決定することができる。
以下の連立方程式は，前の円のパラメータから次の円のパラメータを求める。

@syms r3, x3, y3
(x2, y2) = res2
eq3 = (R - r1 - x3)^2 + y3^2 - (r1 + r3)^2
eq4 = x3^2 + y3^2 - (R - r3)^2
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
res3 = solve([eq3, eq4, eq5], (r3, x3, y3));

これを関数にすれば，累円のパラメータを順に決めることができる。注意点としては x 座標，y 座標の符号が切り替わるときに，用いる関係式が変わるということである。

nextcircle(r2, x2, y2, R=168/2, r1=88/2) = (R*r1*(R^4*r1^2*r2 - R^4*r2^3 - 2*R^3*r1^3*r2 - 2*R^3*r1^2*r2^2 + 2*R^3*r1*r2^3 + R^2*r1^4*r2 + 2*R^2*r1^3*r2^2 - 2*R^2*r1^2*r2^3 + 2*R*r1^3*r2^3 - r1^4*r2^3 + 2*sqrt(r2^2*(R^6*r1^2 - 2*R^6*r1*r2 + R^6*r2^2 - 4*R^5*r1^3 + 6*R^5*r1^2*r2 - 2*R^5*r1*r2^2 + 6*R^4*r1^4 - 8*R^4*r1^3*r2 + 3*R^4*r1^2*r2^2 - 4*R^3*r1^5 + 8*R^3*r1^4*r2 - 4*R^3*r1^3*r2^2 + R^2*r1^6 - 6*R^2*r1^5*r2 + 3*R^2*r1^4*r2^2 + 2*R*r1^6*r2 - 2*R*r1^5*r2^2 + r1^6*r2^2))*sqrt(R*r1*r2*(R - r1 - r2)))/(R^5*r1^3 - 3*R^5*r1^2*r2 + 3*R^5*r1*r2^2 - R^5*r2^3 - 2*R^4*r1^4 + 7*R^4*r1^3*r2 - 4*R^4*r1^2*r2^2 - R^4*r1*r2^3 + R^3*r1^5 - 7*R^3*r1^4*r2 + 2*R^3*r1^3*r2^2 + 3*R^2*r1^5*r2 - 4*R^2*r1^4*r2^2 + 3*R*r1^5*r2^2 + R*r1^4*r2^3 + r1^5*r2^3), R*(-R^6*r1^3 + 3*R^6*r1^2*r2 - 3*R^6*r1*r2^2 + R^6*r2^3 + 3*R^5*r1^4 - 9*R^5*r1^3*r2 + 7*R^5*r1^2*r2^2 - R^5*r1*r2^3 - 3*R^4*r1^5 + 13*R^4*r1^4*r2 - 8*R^4*r1^3*r2^2 + R^3*r1^6 - 11*R^3*r1^5*r2 + 6*R^3*r1^4*r2^2 + 4*R^2*r1^6*r2 - 5*R^2*r1^5*r2^2 - R^2*r1^4*r2^3 + 3*R*r1^6*r2^2 + R*r1^5*r2^3 + 2*R*r1*sqrt(r2^2*(R^6*r1^2 - 2*R^6*r1*r2 + R^6*r2^2 - 4*R^5*r1^3 + 6*R^5*r1^2*r2 - 2*R^5*r1*r2^2 + 6*R^4*r1^4 - 8*R^4*r1^3*r2 + 3*R^4*r1^2*r2^2 - 4*R^3*r1^5 + 8*R^3*r1^4*r2 - 4*R^3*r1^3*r2^2 + R^2*r1^6 - 6*R^2*r1^5*r2 + 3*R^2*r1^4*r2^2 + 2*R*r1^6*r2 - 2*R*r1^5*r2^2 + r1^6*r2^2))*sqrt(R*r1*r2*(R - r1 - r2)) + 2*r1^2*sqrt(r2^2*(R^6*r1^2 - 2*R^6*r1*r2 + R^6*r2^2 - 4*R^5*r1^3 + 6*R^5*r1^2*r2 - 2*R^5*r1*r2^2 + 6*R^4*r1^4 - 8*R^4*r1^3*r2 + 3*R^4*r1^2*r2^2 - 4*R^3*r1^5 + 8*R^3*r1^4*r2 - 4*R^3*r1^3*r2^2 + R^2*r1^6 - 6*R^2*r1^5*r2 + 3*R^2*r1^4*r2^2 + 2*R*r1^6*r2 - 2*R*r1^5*r2^2 + r1^6*r2^2))*sqrt(R*r1*r2*(R - r1 - r2)))/(-R^6*r1^3 + 3*R^6*r1^2*r2 - 3*R^6*r1*r2^2 + R^6*r2^3 + 3*R^5*r1^4 - 10*R^5*r1^3*r2 + 7*R^5*r1^2*r2^2 - 3*R^4*r1^5 + 14*R^4*r1^4*r2 - 6*R^4*r1^3*r2^2 - R^4*r1^2*r2^3 + R^3*r1^6 - 10*R^3*r1^5*r2 + 6*R^3*r1^4*r2^2 + 3*R^2*r1^6*r2 - 7*R^2*r1^5*r2^2 - R^2*r1^4*r2^3 + 3*R*r1^6*r2^2 + r1^6*r2^3), -2*R*r1*sqrt(r2*(R*r1*(R - r1 - r2)*(-R^2*r1 + R^2*r2 + R*r1^2 - 2*R*r1*r2 + r1^2*r2)^2 + (-R^2*r1 + R^2*r2 + R*r1^2 - R*r1*r2 + r1^2*r2)*(R^4*r1^2 - 2*R^4*r1*r2 + R^4*r2^2 - 2*R^3*r1^3 + 6*R^3*r1^2*r2 + R^2*r1^4 - 6*R^2*r1^3*r2 - 2*R^2*r1^2*r2^2 + 2*R*r1^4*r2 + r1^4*r2^2)))/(R^4*r1^2 - 2*R^4*r1*r2 + R^4*r2^2 - 2*R^3*r1^3 + 6*R^3*r1^2*r2 + R^2*r1^4 - 6*R^2*r1^3*r2 - 2*R^2*r1^2*r2^2 + 2*R*r1^4*r2 + r1^4*r2^2) + 2*R*r1*sqrt(R*r1*r2*(R - r1 - r2))*(R*r1 - R*r2 + r1*r2)/(-R^3*r1^2 + 2*R^3*r1*r2 - R^3*r2^2 + R^2*r1^3 - 4*R^2*r1^2*r2 - R^2*r1*r2^2 + 2*R*r1^3*r2 + R*r1^2*r2^2 + r1^3*r2^2));

nextcircle(14/2, 61.6, -46.2)

   (12.157894736842104, 45.09473684210526, -55.92631578947368)

nextcircle(12.157894736842104, 45.09473684210526, -55.92631578947368)

   (23.099999999999962, 10.079999999999686, -60.05999999999997)

nextcircle2(23.099999999999962, 10.079999999999686, -60.05999999999997)

   (38.50000000000007, -39.200000000000394, -23.100000000000193)

nextcircle2(38.50000000000007, -39.200000000000394, -23.100000000000193)

   (32.99999999999988, -21.59999999999963, 46.20000000000031)

3 番目の符号の切り替わり。

@syms r7, x7, y7
@syms r8, x8, y8
eq6 = (R - r1 - x8)^2 + y8^2 - (r1 + r8)^2
eq7 = x8^2 + y8^2 - (R - r8)^2
eq8 = (x7 - x8)^2 + (y7 - y8)^2 - (r7 + r8)^2
res3 = solve([eq6, eq7, eq8], (r8, x8, y8))[1];

nextcircle3(r7, x7, y7, R=168/2, r1=88/2) = ((R^5 - 2*R^4*r1 + R^4*r7 - 3*R^4*x7 + R^3*r1^2 - 2*R^3*r1*r7 + 4*R^3*r1*x7 - R^3*r7^2 - 2*R^3*r7*x7 + 3*R^3*x7^2 + R^3*y7^2 + R^2*r1^2*r7 - R^2*r1^2*x7 + 2*R^2*r1*r7^2 + 4*R^2*r1*r7*x7 - 2*R^2*r1*x7^2 + 2*R^2*r1*y7^2 - R^2*r7^3 + R^2*r7^2*x7 + R^2*r7*x7^2 + R^2*r7*y7^2 - R^2*x7^3 - R^2*x7*y7^2 - R*r1^2*r7^2 - 2*R*r1^2*r7*x7 - R*r1^2*x7^2 - 3*R*r1^2*y7^2 + 2*R*r1*r7^3 - 2*R*r1*r7*x7^2 - 2*R*r1*r7*y7^2 - 2*R*y7*sqrt(R*r1*(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 - 2*R^2*r1*r7 + 2*R^2*r1*x7 - 4*R^2*r7*x7 + 2*R*r1*r7^2 + 4*R*r1*r7*x7 + 2*R*r1*x7^2 + 2*R*r1*y7^2 - 2*R*r7^3 - 2*R*r7^2*x7 + 2*R*r7*x7^2 + 2*R*r7*y7^2 + 2*R*x7^3 + 2*R*x7*y7^2 + 2*r1*r7^3 + 2*r1*r7^2*x7 - 2*r1*r7*x7^2 - 2*r1*r7*y7^2 - 2*r1*x7^3 - 2*r1*x7*y7^2 - r7^4 + 2*r7^2*x7^2 + 2*r7^2*y7^2 - x7^4 - 2*x7^2*y7^2 - y7^4)) - r1^2*r7^3 - r1^2*r7^2*x7 + r1^2*r7*x7^2 + r1^2*r7*y7^2 + r1^2*x7^3 + r1^2*x7*y7^2 + 2*r1*y7*sqrt(R*r1*(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 - 2*R^2*r1*r7 + 2*R^2*r1*x7 - 4*R^2*r7*x7 + 2*R*r1*r7^2 + 4*R*r1*r7*x7 + 2*R*r1*x7^2 + 2*R*r1*y7^2 - 2*R*r7^3 - 2*R*r7^2*x7 + 2*R*r7*x7^2 + 2*R*r7*y7^2 + 2*R*x7^3 + 2*R*x7*y7^2 + 2*r1*r7^3 + 2*r1*r7^2*x7 - 2*r1*r7*x7^2 - 2*r1*r7*y7^2 - 2*r1*x7^3 - 2*r1*x7*y7^2 - r7^4 + 2*r7^2*x7^2 + 2*r7^2*y7^2 - x7^4 - 2*x7^2*y7^2 - y7^4)))/(2*(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 + R^2*r1^2 - 4*R^2*r1*r7 + R^2*r7^2 - 2*R^2*r7*x7 + R^2*x7^2 + 2*R*r1^2*r7 + 2*R*r1^2*x7 - 2*R*r1*r7^2 + 2*R*r1*x7^2 + 4*R*r1*y7^2 + r1^2*r7^2 + 2*r1^2*r7*x7 + r1^2*x7^2)), (R^5 - 4*R^4*r1 + 3*R^4*r7 - R^4*x7 + 3*R^3*r1^2 - 8*R^3*r1*r7 + 2*R^3*r1*x7 + 3*R^3*r7^2 - 2*R^3*r7*x7 - R^3*x7^2 - R^3*y7^2 + 5*R^2*r1^2*r7 + 3*R^2*r1^2*x7 - 4*R^2*r1*r7^2 + 4*R^2*r1*y7^2 + R^2*r7^3 - R^2*r7^2*x7 - R^2*r7*x7^2 - R^2*r7*y7^2 + R^2*x7^3 + R^2*x7*y7^2 + R*r1^2*r7^2 + 2*R*r1^2*r7*x7 + R*r1^2*x7^2 - 3*R*r1^2*y7^2 - 2*R*r1*r7^2*x7 + 2*R*r1*x7^3 + 2*R*r1*x7*y7^2 + 2*R*y7*sqrt(R*r1*(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 - 2*R^2*r1*r7 + 2*R^2*r1*x7 - 4*R^2*r7*x7 + 2*R*r1*r7^2 + 4*R*r1*r7*x7 + 2*R*r1*x7^2 + 2*R*r1*y7^2 - 2*R*r7^3 - 2*R*r7^2*x7 + 2*R*r7*x7^2 + 2*R*r7*y7^2 + 2*R*x7^3 + 2*R*x7*y7^2 + 2*r1*r7^3 + 2*r1*r7^2*x7 - 2*r1*r7*x7^2 - 2*r1*r7*y7^2 - 2*r1*x7^3 - 2*r1*x7*y7^2 - r7^4 + 2*r7^2*x7^2 + 2*r7^2*y7^2 - x7^4 - 2*x7^2*y7^2 - y7^4)) - r1^2*r7^3 - r1^2*r7^2*x7 + r1^2*r7*x7^2 + r1^2*r7*y7^2 + r1^2*x7^3 + r1^2*x7*y7^2 + 2*r1*y7*sqrt(R*r1*(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 - 2*R^2*r1*r7 + 2*R^2*r1*x7 - 4*R^2*r7*x7 + 2*R*r1*r7^2 + 4*R*r1*r7*x7 + 2*R*r1*x7^2 + 2*R*r1*y7^2 - 2*R*r7^3 - 2*R*r7^2*x7 + 2*R*r7*x7^2 + 2*R*r7*y7^2 + 2*R*x7^3 + 2*R*x7*y7^2 + 2*r1*r7^3 + 2*r1*r7^2*x7 - 2*r1*r7*x7^2 - 2*r1*r7*y7^2 - 2*r1*x7^3 - 2*r1*x7*y7^2 - r7^4 + 2*r7^2*x7^2 + 2*r7^2*y7^2 - x7^4 - 2*x7^2*y7^2 - y7^4)))/(2*(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 + R^2*r1^2 - 4*R^2*r1*r7 + R^2*r7^2 - 2*R^2*r7*x7 + R^2*x7^2 + 2*R*r1^2*r7 + 2*R*r1^2*x7 - 2*R*r1*r7^2 + 2*R*r1*x7^2 + 4*R*r1*y7^2 + r1^2*r7^2 + 2*r1^2*r7*x7 + r1^2*x7^2)), 2*R*r1*y7*(R*r1 - R*r7 - R*x7 + r1*r7 + r1*x7 - r7^2 + x7^2 + y7^2)/(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 + R^2*r1^2 - 4*R^2*r1*r7 + R^2*r7^2 - 2*R^2*r7*x7 + R^2*x7^2 + 2*R*r1^2*r7 + 2*R*r1^2*x7 - 2*R*r1*r7^2 + 2*R*r1*x7^2 + 4*R*r1*y7^2 + r1^2*r7^2 + 2*r1^2*r7*x7 + r1^2*x7^2) - sqrt(-R*r1*(-R^2 - 2*R*r7 - r7^2 + x7^2 + y7^2)*(R^2 - 2*R*r1 - 2*R*x7 + 2*r1*r7 + 2*r1*x7 - r7^2 + x7^2 + y7^2))*(-R^2 + R*r1 - R*r7 + R*x7 + r1*r7 + r1*x7)/(R^4 - 2*R^3*r1 + 2*R^3*r7 - 2*R^3*x7 + R^2*r1^2 - 4*R^2*r1*r7 + R^2*r7^2 - 2*R^2*r7*x7 + R^2*x7^2 + 2*R*r1^2*r7 + 2*R*r1^2*x7 - 2*R*r1*r7^2 + 2*R*r1*x7^2 + 4*R*r1*y7^2 + r1^2*r7^2 + 2*r1^2*r7*x7 + r1^2*x7^2));

nextcircle3(32.99999999999988, -21.59999999999963, 46.20000000000031)

   (17.769230769230685, 27.138461538461794, 60.41538461538458)

nextcircle3(17.769230769230685, 27.138461538461794, 60.41538461538458)

   (9.624999999999968, 53.20000000000012, 51.97499999999994)

nextcircle3(9.624999999999968, 53.20000000000012, 51.97499999999994)

   (5.774999999999992, 65.52000000000011, 42.73499999999996)

nextcircle3(5.774999999999992, 65.52000000000011, 42.73499999999996)

   (3.786885245901612, 71.88196721311486, 35.59672131147536)

nextcircle3(3.786885245901612, 71.88196721311486, 35.59672131147536)

   (2.6551724137930997, 75.50344827586207, 30.268965517241313)

nextcircle3(2.6551724137930997, 75.50344827586207, 30.268965517241313)

   (1.9576271186440766, 77.73559322033896, 26.232203389830453)

nextcircle3(1.9576271186440766, 77.73559322033896, 26.232203389830453)

   (1.4999999999999876, 79.20000000000003, 23.099999999999948)

半径が 1.5 の円が出現した。これは，初円を 1 番目としたときの 12番目の円である。

4番目の円から12番目の円までの半径はすべて，実は分子が 231 の分数である。分母は 6, 7, 13, 24, 40, 61, 87, 118, 154 である。
分母は第 2 階差数列で，n 番目の項が 5//2*n^2 - 13//2*n + 10 の数列になっている。

for n = 1:9
   分母 = 5//2*n^2 - 13//2*n + 10
   @printf("n = %g;  %g 番目の円:  分母 = %g; 半径 = %g\n", n, n + 3, 分母, 231/分母)
end

   n = 1;   4 番目の円:  分母 =   6; 半径 = 38.5
   n = 2;   5 番目の円:  分母 =   7; 半径 = 33
   n = 3;   6 番目の円:  分母 =  13; 半径 = 17.7692
   n = 4;   7 番目の円:  分母 =  24; 半径 = 9.625
   n = 5;   8 番目の円:  分母 =  40; 半径 = 5.775
   n = 6;   9 番目の円:  分母 =  61; 半径 = 3.78689
   n = 7;  10 番目の円:  分母 =  87; 半径 = 2.65517
   n = 8;  11 番目の円:  分母 = 118; 半径 = 1.95763
   n = 9;  12 番目の円:  分母 = 154; 半径 = 1.5

なお，初円，二円，三円も分数で表される。

231/33, 231/19, 231/10

   (7.0, 12.157894736842104, 23.1)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2) = (168, 88, 14) ./ 2
   (x2, y2) = ((-R^2 + R*r1 + R*r2 + r1*r2)/(-R + r1), 2*sqrt(-R*r1*r2*(-R + r1 + r2))/(-R + r1))
   (r3, x3, y3) = (R*r1*(R^4*r1^2*r2 - R^4*r2^3 - 2*R^3*r1^3*r2 - 2*R^3*r1^2*r2^2 + 2*R^3*r1*r2^3 + R^2*r1^4*r2 + 2*R^2*r1^3*r2^2 - 2*R^2*r1^2*r2^3 + 2*R*r1^3*r2^3 - r1^4*r2^3 + 2*sqrt(r2^2*(R^6*r1^2 - 2*R^6*r1*r2 + R^6*r2^2 - 4*R^5*r1^3 + 6*R^5*r1^2*r2 - 2*R^5*r1*r2^2 + 6*R^4*r1^4 - 8*R^4*r1^3*r2 + 3*R^4*r1^2*r2^2 - 4*R^3*r1^5 + 8*R^3*r1^4*r2 - 4*R^3*r1^3*r2^2 + R^2*r1^6 - 6*R^2*r1^5*r2 + 3*R^2*r1^4*r2^2 + 2*R*r1^6*r2 - 2*R*r1^5*r2^2 + r1^6*r2^2))*sqrt(R*r1*r2*(R - r1 - r2)))/(R^5*r1^3 - 3*R^5*r1^2*r2 + 3*R^5*r1*r2^2 - R^5*r2^3 - 2*R^4*r1^4 + 7*R^4*r1^3*r2 - 4*R^4*r1^2*r2^2 - R^4*r1*r2^3 + R^3*r1^5 - 7*R^3*r1^4*r2 + 2*R^3*r1^3*r2^2 + 3*R^2*r1^5*r2 - 4*R^2*r1^4*r2^2 + 3*R*r1^5*r2^2 + R*r1^4*r2^3 + r1^5*r2^3), R*(-R^6*r1^3 + 3*R^6*r1^2*r2 - 3*R^6*r1*r2^2 + R^6*r2^3 + 3*R^5*r1^4 - 9*R^5*r1^3*r2 + 7*R^5*r1^2*r2^2 - R^5*r1*r2^3 - 3*R^4*r1^5 + 13*R^4*r1^4*r2 - 8*R^4*r1^3*r2^2 + R^3*r1^6 - 11*R^3*r1^5*r2 + 6*R^3*r1^4*r2^2 + 4*R^2*r1^6*r2 - 5*R^2*r1^5*r2^2 - R^2*r1^4*r2^3 + 3*R*r1^6*r2^2 + R*r1^5*r2^3 + 2*R*r1*sqrt(r2^2*(R^6*r1^2 - 2*R^6*r1*r2 + R^6*r2^2 - 4*R^5*r1^3 + 6*R^5*r1^2*r2 - 2*R^5*r1*r2^2 + 6*R^4*r1^4 - 8*R^4*r1^3*r2 + 3*R^4*r1^2*r2^2 - 4*R^3*r1^5 + 8*R^3*r1^4*r2 - 4*R^3*r1^3*r2^2 + R^2*r1^6 - 6*R^2*r1^5*r2 + 3*R^2*r1^4*r2^2 + 2*R*r1^6*r2 - 2*R*r1^5*r2^2 + r1^6*r2^2))*sqrt(R*r1*r2*(R - r1 - r2)) + 2*r1^2*sqrt(r2^2*(R^6*r1^2 - 2*R^6*r1*r2 + R^6*r2^2 - 4*R^5*r1^3 + 6*R^5*r1^2*r2 - 2*R^5*r1*r2^2 + 6*R^4*r1^4 - 8*R^4*r1^3*r2 + 3*R^4*r1^2*r2^2 - 4*R^3*r1^5 + 8*R^3*r1^4*r2 - 4*R^3*r1^3*r2^2 + R^2*r1^6 - 6*R^2*r1^5*r2 + 3*R^2*r1^4*r2^2 + 2*R*r1^6*r2 - 2*R*r1^5*r2^2 + r1^6*r2^2))*sqrt(R*r1*r2*(R - r1 - r2)))/(-R^6*r1^3 + 3*R^6*r1^2*r2 - 3*R^6*r1*r2^2 + R^6*r2^3 + 3*R^5*r1^4 - 10*R^5*r1^3*r2 + 7*R^5*r1^2*r2^2 - 3*R^4*r1^5 + 14*R^4*r1^4*r2 - 6*R^4*r1^3*r2^2 - R^4*r1^2*r2^3 + R^3*r1^6 - 10*R^3*r1^5*r2 + 6*R^3*r1^4*r2^2 + 3*R^2*r1^6*r2 - 7*R^2*r1^5*r2^2 - R^2*r1^4*r2^3 + 3*R*r1^6*r2^2 + r1^6*r2^3), -2*R*r1*sqrt(r2*(R*r1*(R - r1 - r2)*(-R^2*r1 + R^2*r2 + R*r1^2 - 2*R*r1*r2 + r1^2*r2)^2 + (-R^2*r1 + R^2*r2 + R*r1^2 - R*r1*r2 + r1^2*r2)*(R^4*r1^2 - 2*R^4*r1*r2 + R^4*r2^2 - 2*R^3*r1^3 + 6*R^3*r1^2*r2 + R^2*r1^4 - 6*R^2*r1^3*r2 - 2*R^2*r1^2*r2^2 + 2*R*r1^4*r2 + r1^4*r2^2)))/(R^4*r1^2 - 2*R^4*r1*r2 + R^4*r2^2 - 2*R^3*r1^3 + 6*R^3*r1^2*r2 + R^2*r1^4 - 6*R^2*r1^3*r2 - 2*R^2*r1^2*r2^2 + 2*R*r1^4*r2 + r1^4*r2^2) + 2*R*r1*sqrt(R*r1*r2*(R - r1 - r2))*(R*r1 - R*r2 + r1*r2)/(-R^3*r1^2 + 2*R^3*r1*r2 - R^3*r2^2 + R^2*r1^3 - 4*R^2*r1^2*r2 - R^2*r1*r2^2 + 2*R*r1^3*r2 + R*r1^2*r2^2 + r1^3*r2^2))
   (x2, y2) = (61.6, -46.2)
   (r3, x3, y3) = (12.157894736842104, 45.09473684210526, -55.92631578947368)
   (r4, x4, y4) = (23.099999999999962, 10.079999999999686, -60.05999999999997)
   (r5, x5, y5) = (38.50000000000007, -39.200000000000394, -23.100000000000193)
   (r6, x6, y6) = (32.99999999999988, -21.59999999999963, 46.20000000000031)
   (r7, x7, y7) = (17.769230769230685, 27.138461538461794, 60.41538461538458)
   (r8, x8, y8) = (9.624999999999968, 53.20000000000012, 51.97499999999994)
   (r9, x9, y9) = (5.774999999999992, 65.52000000000011, 42.73499999999996)
   (r10, x10, y10) = (3.786885245901612, 71.88196721311486, 35.59672131147536)
   (r11, x11, y11) = (2.6551724137930997, 75.50344827586207, 30.268965517241313)
   (r12, x12, y12) = (1.9576271186440766, 77.73559322033896, 26.232203389830453)
   (r13, x13, y13) = (1.4999999999999876, 79.20000000000003, 23.099999999999948)
   plot()
   circlef(0, 0, R)
   circlef(R - r1, 0, r1, 1)
   circlef(x2, y2, r2, 2)
   circlef(x3, y3, r3, 3)
   circlef(x4, y4, r4, 4)
   circlef(x5, y5, r5, 5)
   circlef(x6, y6, r6, 6)
   circlef(x7, y7, r7, 7)
   circlef(x8, y8, r8, 8)
   circlef(x9, y9, r9, 9)
   circlef(x10, y10, r10, 10)
   circlef(x11, y11, r11, 11)
   circlef(x12, y12, r12, 12)
   circlef(x13, y13, r13, 13)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R - r1, 0, "小円:r1,(R-r1,0)", :white, :center, delta=-delta/2)
       point(x2, y2, "初円:r2,(x2,y2)", :white, :right, delta=-delta/2)
       point(x3, y3, "二円:r3,(x3,y3)", :white, :right, delta=-delta/2)
       point(x13, y13, "終円:r13,(x13,y13) ", :white, :right, delta=-delta/2)
   end
end;