算額（その1117）

算額（その1117）

十 群馬県甘楽郡妙義町妙義 妙義神社 寛政9年(1797)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

四 岩手県花巻市南笹間 東光寺 慶応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;

算額（その1116）

算額（その1116）

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

盤の上に大球を載せ，その周りに小球を数個互いに接するように並べる。大球の直径と小球の個数が与えられたとき，小球の直径を求めよ。

注：単に大球と書いているが，図にあるように，半球である

大球の半径と中心座標を R, (0, 0, 0)
小球の半径と中心座標を r, (x, 0, r), (x*cos(2π/n), x*sin(2π/n)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms n::positiveinteger, R::positive, r::positive,
     x::positive
eq1 = x^2 + r^2 - (R + r)^2
eq2 = x*sind(Sym(180)/n) - r;

(r, x) = solve([eq1, eq2], (r, x))[2];
r |> println
x |> println

   R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n)
   R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))

小球の半径は大球の半径の (sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n) 倍である。

n = 12
R = 5
r = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n)
x = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))
@printf("大球の直径が %g で，%i 個の小球を並べるとき，小球の直径は %g である。\n小球の中心座標は以下の通りである。\n", 2R, n, 2r)
for i = 1:n
   println((x*cosd(360*i/n), x*sind(360*i/n)))
end

   大球の直径が 10 で，12 個の小球を並べるとき，小球の直径は 3.34335 である。
   小球の中心座標は以下の通りである。
   (5.593527388798881, 3.229424543642579)
   (3.229424543642579, 5.593527388798881)
   (0.0, 6.458849087285158)
   (-3.229424543642579, 5.593527388798881)
   (-5.593527388798881, 3.229424543642579)
   (-6.458849087285158, 0.0)
   (-5.593527388798881, -3.229424543642579)
   (-3.229424543642579, -5.593527388798881)
   (0.0, -6.458849087285158)
   (3.229424543642579, -5.593527388798881)
   (5.593527388798881, -3.229424543642579)
   (6.458849087285158, 0.0)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = 12
   R = 5
   r = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))*sin(pi/n)
   x = R*(sqrt(sin(pi/n)^2 + 1) + sin(pi/n))
   plot(showaxis=false)
   circle(0, 0, R, :blue, beginangle=0, endangle=180)
   circle(x, r, r)
   segment(x, R, x, -3R, :green)
   segment(R, R, R, -3R, :blue)
   segment(x + r, R, x + r, -3R, :red)
   segment(x - r, R, x - r, -3R, :red)
   segment(0, 0, x, r)
   segment(x, 0, x, r, :green)
   segment(-R, 0, R + 2r, 0)
   base = -1.7R
   circle(0, base, R, :blue)
   circle(0, base, x, :green)
   for i = 0:n - 1
       (x0, y0) = (x*cosd(i*360/n), base + x*sind(i*360/n))
       circle(x0, y0, r)
   end
   i = 0.5
   (x0, y0) = (x*cosd(i*360/n), base + x*sind(i*360/n))
   segment(0, base, x + r, base)
   segment(0, base, x, y0)
   segment(x, base, x0, y0)
   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, 0, " R", :blue, :left, delta=-delta)
       point(x, r, " (x,r)", :red, :left, :vcenter)
       circle(0, base, 5delta, beginangle=0, endangle=360/2n)
       circle(0, base, 5.3delta, beginangle=0, endangle=360/2n)
       point(5delta, base + 1.5delta, "θ=π/n", mark=false)
       point(3delta, -20delta, "上から見た図", mark=false)
       point(3delta, 8delta, "横から見た図", mark=false)
       plot!(xlims=(-5delta, x + r + 5delta), ylims=(base - 8delta, R + 2delta))
   end
end;

算額（その1115）

算額（その1115）

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

円錐台の中に直径 4 寸の大球 2 個と直径 2 寸の小球 2 個が入っている。円錐台の高さを求めよ。
注：大球，小球は円錐台の側面にそれぞれ 2 点で接し，小球は円錐台の上面に接している。

2 個の小球が重なって見える方向から見た図

2 個の大球が重なって見える方向から見た図

円錐の高さ，底面の円の半径を h，a
大球の半径と中心座標を r1, (r1, 0, r1)
小球の半径と中心座標を r2, (0, r2, z2)
とおき，以下の連立方程式を解く。
円錐台の高さは，小球の中心の z 座標値に小球の半径を加えたもの z2 + r2 である。

include("julia-source.txt")

using SymPy
@syms a::positice, h::positive, r1::positive,
     r2::positive, z2::positive
eq1 = r1^2 + r2^2 + (z2 - r1)^2 - (r1 + r2)^2
eq2 = a + h - sqrt(a^2 + h^2) - 2r1
eq3 = dist2(0, h, a, 0, r2, z2, r2)
eq4 = r2*(h + a + sqrt(a^2 + h^2)) - a*h;

z2 は eq1 を解くことにより求めることができる。

ans_z2 = solve(eq1,z2)[2]
ans_z2 |> println
ans_z2 + r2 |> println

   sqrt(2)*sqrt(r1)*sqrt(r2) + r1
   sqrt(2)*sqrt(r1)*sqrt(r2) + r1 + r2

円錐台の高さは sqrt(2r1*r2) + r1 + r2 である。
大球，小球の半径が 4/2 寸, 2/2 寸のとき，円錐台の高さは 5 寸である。

ans_z2(r1 => 4//2, r2 => 2//2) |> println
(ans_z2 + r2)(r1 => 4//2, r2 => 2//2) |> println

   4
   5

この連立方程式は，変数名のままでは適切な解が得られない。r1, r2 に特定の数値を与え，a, h を求める。

eq12 = eq2(z2 => ans_z2)
eq13 = eq3(z2 => ans_z2)/a 
eq14 = eq4(z2 => ans_z2)

eq12 = eq12(r1 => 4//2, r2 => 2//2)
eq13 = eq13(r1 => 4//2, r2 => 2//2)
eq14 = eq14(r1 => 4//2, r2 => 2//2)

円錐台の底面の円の直径は 8 寸，円錐の高さは 6 寸である。

res = solve([eq12, eq13], (a, h))[1]

   (8, 6)

function draw(type, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, a, h, z2) = (4/2, 2/2, 8, 6, 4)
   b = (h - (z2 + r2))/h*a
   plot([a, 0, -a, a], [0, h, 0, 0], color=:green, lw=0.5)
   if type == 1
       circle2(r1, r1, r1, :blue)
       circle(0, z2, r2)
   else
       circle(0, r1, r1, :blue)
       circle2(r2, z2, r2)
   end
   segment(-b, z2 + r2, b, z2 + r2)
   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)
       if type == 1
           point(r1, r1, "大球:r1,(r1,r1)", :blue, :left, :bottom, delta=delta/2)
           point(0, z2, "小球:r2,(0, z2)", :red, :left, :bottom, delta=delta/2)
       else
           point(0, r1, "大球:r1,(0,r1)", :blue, :left, :bottom, delta=delta/2)
           point(r2, z2, "小球:r2,(r2, z2)", :red, :left, :bottom, delta=delta/2)
       end
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(0, h, "h", :green, :center, :bottom, delta=delta/2)
       point(0, z2 + r2, " z2+r2", :green, :center, :bottom, delta=delta/2)
   end
end;