算額（その194）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（89）
長野県中野市田上 田上観音堂 文化6年(1809)
問3. 大円の中に，甲，乙，丙，丁，戊，己の 5 個の円が外円に内接すると同時に隣同士の円は外接している。さらに，子，丑，寅，卯，辰の 5 個の円が甲円と乙，丙，丁，戊，己円と外接している。大円の径が 76 寸のとき，辰円の径はいかほどか。
大円の中心座標，半径を (0, R0, R0) ただし R0 = 2
甲円の中心座標，半径を (0, R1, R1) および (0, 3R1, R1) ただし R1 = 1
乙円の中心座標，半径を (X2, Y2, R2)
丙円の中心座標，半径を (X3, Y3, R3)
丁円の中心座標，半径を (X4, Y4, R4)
戊円の中心座標，半径を (X5, Y5, R5)
己円の中心座標，半径を (X6, Y6, R6)
子円の中心座標，半径を (x2, y2, r2)
丑円の中心座標，半径を (x3, y3, r3)
寅円の中心座標，半径を (x4, y4, r4)
卯円の中心座標，半径を (x5, y5, r5)
辰円の中心座標，半径を (x6, y6, r6)
とおき方程式を解く。なお，解として得られるのは大円の半径を 2 としたときの各円の半径なので，元の大円の大きさ，単位（寸）での直径はそれぞれを 38 倍する。大円の半径は 2 の倍数のほうが良いが，大きさは任意に指定できない。SymPy のせいであろうが，16 以上にすると平方根の計算においてエラーになる。
using SymPy
@syms R0::positive, R1::positive,
X2::positive, Y2::positive, R2::positive,
x2::positive, y2::positive, r2::positive;
R0 = 2 # 大円の半径
R1 = R0 // 2 # 甲円の半径
X2 = R0 - R2
Y2 = 2R1
eq1 = X2^2 + R1^2 - (R1 + R2)^2
R2 = solve(eq1, R2)[1] # 乙円の半径
$\frac{2}{3}$
R0 - R0//3, Y2, R2 # 乙円の中心座標と半径
(4//3, 2//1, 2/3)
eq2 = x2^2 + (y2 - R1)^2 - (R1 + r2)^2
eq3 = x2^2 + (3R1 - y2)^2 - (R1 + r2)^2
eq4 = (R0 - R2 - x2)^2 + (y2 - 2R1)^2 - (r2 + R2)^2;
(r2, x2, y2) = solve([eq2, eq3, eq4], (r2, x2, y2))[1] # 子円の中心座標と半径
(2/15, 8/15, 2)
乙円から始めて，丙円以降の中心座標と半径を求める漸化式を求める。
@syms Xn, Yn, Rn, Xnp1, Ynp1, Rnp1,
xn, yn, rn, xnp1, ynp1, rnp1;
eq5 = (Xn - Xnp1)^2 + (Yn - Ynp1)^2 - (Rn + Rnp1)^2 |> expand
eq6 = Xnp1^2 + (R1 - Ynp1)^2 - (R1 + Rnp1)^2 |> expand
eq7 = Xnp1^2 + (2R1 - Ynp1)^2 - (R0 - Rnp1)^2 |> expand;
res = solve([eq5, eq6, eq7], (Xnp1, Ynp1, Rnp1));
println(res[2])
出力省略
2番めの組が適切解である。これに基づいて，現在の円の中心座標と半径を与えて，次の円の中心座標と半径を求める関数を定義する。
NEXTCIRCLE(Xn, Yn, Rn) = return ((-4*Rn^2*Xn - 4*Rn*Xn + sqrt(2)*Rn*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2) + 4*Xn^3 + 4*Xn*Yn^2 - 12*Xn*Yn + 16*Xn + 3*sqrt(2)*Yn*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2) - 4*sqrt(2)*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2))/(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16), 3*(-Rn^3 - 3*Rn^2*Yn + 4*Rn^2 + Rn*Xn^2 + Rn*Yn^2 + 3*Xn^2*Yn + 4*Xn^2 - 2*sqrt(2)*Xn*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2) + 3*Yn^3 - 4*Yn^2)/(2*(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16)), -sqrt(2)*Xn*sqrt(-(Rn^2 - 2*Rn - Xn^2 - Yn^2 + 2*Yn)*(Rn^2 + 4*Rn - Xn^2 - Yn^2 + 4*Yn))/(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16) - (Rn^3 + 3*Rn^2*Yn - 4*Rn^2 - Rn*Xn^2 - Rn*Yn^2 - 3*Xn^2*Yn - 4*Xn^2 - 3*Yn^3 + 4*Yn^2)/(2*(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16)));
下側の甲円と乙円，丙円から始めて，丑円以降の中心座標と半径を求める漸化式を求める。
eq8 = xnp1^2 + (R1 -ynp1)^2 - (R1 + rnp1)^2 |> expand
eq9 = (Xn - xnp1)^2 + (Yn - ynp1)^2 - (Rn + rnp1)^2 |> expand
eq10 = (Xnp1 - xnp1)^2 + (ynp1 - Ynp1)^2 - (Rnp1 + rnp1)^2 |> expand
res2 = solve([eq8, eq9, eq10], (xnp1, ynp1, rnp1));
println(res2[1])
出力省略
1番めの組が適切解である。これに基づいて，丑円以降の円の中心座標と半径を求める関数を定義する。
function nextcircle(Xn, Yn, Rn, Xnp1, Ynp1, Rnp1)
term = sqrt(-Rn^4*Rnp1^2 + 2*Rn^4*Rnp1 + Rn^4*Xnp1^2 + Rn^4*Ynp1^2 - 2*Rn^4*Ynp1 + 2*Rn^3*Rnp1^3 - 2*Rn^3*Rnp1^2 - 2*Rn^3*Rnp1*Xnp1^2 - 2*Rn^3*Rnp1*Ynp1^2 + 4*Rn^3*Rnp1*Ynp1 - 4*Rn^3*Rnp1 - 2*Rn^3*Xnp1^2 - 2*Rn^3*Ynp1^2 + 4*Rn^3*Ynp1 - Rn^2*Rnp1^4 - 2*Rn^2*Rnp1^3 + 2*Rn^2*Rnp1^2*Xn^2 - 2*Rn^2*Rnp1^2*Xn*Xnp1 + 2*Rn^2*Rnp1^2*Xnp1^2 + 2*Rn^2*Rnp1^2*Yn^2 - 2*Rn^2*Rnp1^2*Yn*Ynp1 - 2*Rn^2*Rnp1^2*Yn + 2*Rn^2*Rnp1^2*Ynp1^2 - 2*Rn^2*Rnp1^2*Ynp1 + 8*Rn^2*Rnp1^2 - 4*Rn^2*Rnp1*Xn^2 + 4*Rn^2*Rnp1*Xn*Xnp1 + 2*Rn^2*Rnp1*Xnp1^2 - 4*Rn^2*Rnp1*Yn^2 + 4*Rn^2*Rnp1*Yn*Ynp1 + 4*Rn^2*Rnp1*Yn + 2*Rn^2*Rnp1*Ynp1^2 - 8*Rn^2*Rnp1*Ynp1 - 2*Rn^2*Xn^2*Xnp1^2 - 2*Rn^2*Xn^2*Ynp1^2 + 4*Rn^2*Xn^2*Ynp1 + 2*Rn^2*Xn*Xnp1^3 + 2*Rn^2*Xn*Xnp1*Ynp1^2 - 4*Rn^2*Xn*Xnp1*Ynp1 - Rn^2*Xnp1^4 - 2*Rn^2*Xnp1^2*Yn^2 + 2*Rn^2*Xnp1^2*Yn*Ynp1 + 2*Rn^2*Xnp1^2*Yn - 2*Rn^2*Xnp1^2*Ynp1^2 + 2*Rn^2*Xnp1^2*Ynp1 - 2*Rn^2*Yn^2*Ynp1^2 + 4*Rn^2*Yn^2*Ynp1 + 2*Rn^2*Yn*Ynp1^3 - 2*Rn^2*Yn*Ynp1^2 - 4*Rn^2*Yn*Ynp1 - Rn^2*Ynp1^4 + 2*Rn^2*Ynp1^3 + 2*Rn*Rnp1^4 - 2*Rn*Rnp1^3*Xn^2 - 2*Rn*Rnp1^3*Yn^2 + 4*Rn*Rnp1^3*Yn - 4*Rn*Rnp1^3 + 2*Rn*Rnp1^2*Xn^2 + 4*Rn*Rnp1^2*Xn*Xnp1 - 4*Rn*Rnp1^2*Xnp1^2 + 2*Rn*Rnp1^2*Yn^2 + 4*Rn*Rnp1^2*Yn*Ynp1 - 8*Rn*Rnp1^2*Yn - 4*Rn*Rnp1^2*Ynp1^2 + 4*Rn*Rnp1^2*Ynp1 + 2*Rn*Rnp1*Xn^2*Xnp1^2 + 2*Rn*Rnp1*Xn^2*Ynp1^2 - 4*Rn*Rnp1*Xn^2*Ynp1 + 4*Rn*Rnp1*Xn^2 - 8*Rn*Rnp1*Xn*Xnp1 + 2*Rn*Rnp1*Xnp1^2*Yn^2 - 4*Rn*Rnp1*Xnp1^2*Yn + 4*Rn*Rnp1*Xnp1^2 + 2*Rn*Rnp1*Yn^2*Ynp1^2 - 4*Rn*Rnp1*Yn^2*Ynp1 + 4*Rn*Rnp1*Yn^2 - 4*Rn*Rnp1*Yn*Ynp1^2 + 4*Rn*Rnp1*Ynp1^2 + 2*Rn*Xn^2*Xnp1^2 + 2*Rn*Xn^2*Ynp1^2 - 4*Rn*Xn^2*Ynp1 - 4*Rn*Xn*Xnp1^3 - 4*Rn*Xn*Xnp1*Ynp1^2 + 8*Rn*Xn*Xnp1*Ynp1 + 2*Rn*Xnp1^4 + 2*Rn*Xnp1^2*Yn^2 - 4*Rn*Xnp1^2*Yn*Ynp1 + 4*Rn*Xnp1^2*Ynp1^2 - 4*Rn*Xnp1^2*Ynp1 + 2*Rn*Yn^2*Ynp1^2 - 4*Rn*Yn^2*Ynp1 - 4*Rn*Yn*Ynp1^3 + 8*Rn*Yn*Ynp1^2 + 2*Rn*Ynp1^4 - 4*Rn*Ynp1^3 + Rnp1^4*Xn^2 + Rnp1^4*Yn^2 - 2*Rnp1^4*Yn - 2*Rnp1^3*Xn^2 - 2*Rnp1^3*Yn^2 + 4*Rnp1^3*Yn - Rnp1^2*Xn^4 + 2*Rnp1^2*Xn^3*Xnp1 - 2*Rnp1^2*Xn^2*Xnp1^2 - 2*Rnp1^2*Xn^2*Yn^2 + 2*Rnp1^2*Xn^2*Yn*Ynp1 + 2*Rnp1^2*Xn^2*Yn - 2*Rnp1^2*Xn^2*Ynp1^2 + 2*Rnp1^2*Xn^2*Ynp1 + 2*Rnp1^2*Xn*Xnp1*Yn^2 - 4*Rnp1^2*Xn*Xnp1*Yn - 2*Rnp1^2*Xnp1^2*Yn^2 + 4*Rnp1^2*Xnp1^2*Yn - Rnp1^2*Yn^4 + 2*Rnp1^2*Yn^3*Ynp1 + 2*Rnp1^2*Yn^3 - 2*Rnp1^2*Yn^2*Ynp1^2 - 2*Rnp1^2*Yn^2*Ynp1 + 4*Rnp1^2*Yn*Ynp1^2 - 4*Rnp1^2*Yn*Ynp1 + 2*Rnp1*Xn^4 - 4*Rnp1*Xn^3*Xnp1 + 2*Rnp1*Xn^2*Xnp1^2 + 4*Rnp1*Xn^2*Yn^2 - 4*Rnp1*Xn^2*Yn*Ynp1 - 4*Rnp1*Xn^2*Yn + 2*Rnp1*Xn^2*Ynp1^2 - 4*Rnp1*Xn*Xnp1*Yn^2 + 8*Rnp1*Xn*Xnp1*Yn + 2*Rnp1*Xnp1^2*Yn^2 - 4*Rnp1*Xnp1^2*Yn + 2*Rnp1*Yn^4 - 4*Rnp1*Yn^3*Ynp1 - 4*Rnp1*Yn^3 + 2*Rnp1*Yn^2*Ynp1^2 + 8*Rnp1*Yn^2*Ynp1 - 4*Rnp1*Yn*Ynp1^2 + Xn^4*Xnp1^2 + Xn^4*Ynp1^2 - 2*Xn^4*Ynp1 - 2*Xn^3*Xnp1^3 - 2*Xn^3*Xnp1*Ynp1^2 + 4*Xn^3*Xnp1*Ynp1 + Xn^2*Xnp1^4 + 2*Xn^2*Xnp1^2*Yn^2 - 2*Xn^2*Xnp1^2*Yn*Ynp1 - 2*Xn^2*Xnp1^2*Yn + 2*Xn^2*Xnp1^2*Ynp1^2 - 2*Xn^2*Xnp1^2*Ynp1 + 2*Xn^2*Yn^2*Ynp1^2 - 4*Xn^2*Yn^2*Ynp1 - 2*Xn^2*Yn*Ynp1^3 + 2*Xn^2*Yn*Ynp1^2 + 4*Xn^2*Yn*Ynp1 + Xn^2*Ynp1^4 - 2*Xn^2*Ynp1^3 - 2*Xn*Xnp1^3*Yn^2 + 4*Xn*Xnp1^3*Yn - 2*Xn*Xnp1*Yn^2*Ynp1^2 + 4*Xn*Xnp1*Yn^2*Ynp1 + 4*Xn*Xnp1*Yn*Ynp1^2 - 8*Xn*Xnp1*Yn*Ynp1 + Xnp1^4*Yn^2 - 2*Xnp1^4*Yn + Xnp1^2*Yn^4 - 2*Xnp1^2*Yn^3*Ynp1 - 2*Xnp1^2*Yn^3 + 2*Xnp1^2*Yn^2*Ynp1^2 + 2*Xnp1^2*Yn^2*Ynp1 - 4*Xnp1^2*Yn*Ynp1^2 + 4*Xnp1^2*Yn*Ynp1 + Yn^4*Ynp1^2 - 2*Yn^4*Ynp1 - 2*Yn^3*Ynp1^3 + 2*Yn^3*Ynp1^2 + 4*Yn^3*Ynp1 + Yn^2*Ynp1^4 + 2*Yn^2*Ynp1^3 - 8*Yn^2*Ynp1^2 - 2*Yn*Ynp1^4 + 4*Yn*Ynp1^3)
return ((Rn^3*Rnp1*Xnp1 - Rn^3*Xnp1 - Rn^2*Rnp1^2*Xn - Rn^2*Rnp1^2*Xnp1 + 2*Rn^2*Rnp1*Xn - Rn^2*Rnp1*Xnp1 + Rn^2*Xn*Ynp1^2 - 2*Rn^2*Xn*Ynp1 + Rn^2*Xnp1^3 - Rn^2*Xnp1*Yn*Ynp1 + Rn^2*Xnp1*Yn + Rn^2*Xnp1*Ynp1^2 - Rn^2*Xnp1*Ynp1 + 2*Rn^2*Xnp1 + Rn*Rnp1^3*Xn - Rn*Rnp1^2*Xn + 2*Rn*Rnp1^2*Xnp1 - Rn*Rnp1*Xn^2*Xnp1 - Rn*Rnp1*Xn*Xnp1^2 - Rn*Rnp1*Xn*Ynp1^2 + 2*Rn*Rnp1*Xn*Ynp1 - 2*Rn*Rnp1*Xn - Rn*Rnp1*Xnp1*Yn^2 + 2*Rn*Rnp1*Xnp1*Yn - 2*Rn*Rnp1*Xnp1 + Rn*Xn^2*Xnp1 + Rn*Xn*Xnp1^2 - Rn*Xn*Ynp1^2 + 2*Rn*Xn*Ynp1 - 2*Rn*Xnp1^3 + Rn*Xnp1*Yn^2 + 2*Rn*Xnp1*Yn*Ynp1 - 4*Rn*Xnp1*Yn - 2*Rn*Xnp1*Ynp1^2 + 2*Rn*Xnp1*Ynp1 - Rn*Ynp1*term + Rn*term - Rnp1^3*Xn + Rnp1^2*Xn^3 + Rnp1^2*Xn*Yn^2 - Rnp1^2*Xn*Yn*Ynp1 - Rnp1^2*Xn*Yn + Rnp1^2*Xn*Ynp1 + 2*Rnp1^2*Xn + Rnp1^2*Xnp1*Yn^2 - 2*Rnp1^2*Xnp1*Yn - 2*Rnp1*Xn^3 + Rnp1*Xn^2*Xnp1 + Rnp1*Xn*Xnp1^2 - 2*Rnp1*Xn*Yn^2 + 2*Rnp1*Xn*Yn*Ynp1 + 2*Rnp1*Xn*Yn + Rnp1*Xn*Ynp1^2 - 4*Rnp1*Xn*Ynp1 - Rnp1*Xnp1*Yn^2 + 2*Rnp1*Xnp1*Yn + Rnp1*Yn*term - Rnp1*term - Xn^3*Ynp1^2 + 2*Xn^3*Ynp1 + Xn^2*Xnp1*Yn*Ynp1 - Xn^2*Xnp1*Yn - Xn^2*Xnp1*Ynp1 + Xn*Xnp1^2*Yn*Ynp1 - Xn*Xnp1^2*Yn - Xn*Xnp1^2*Ynp1 - Xn*Yn^2*Ynp1^2 + 2*Xn*Yn^2*Ynp1 + Xn*Yn*Ynp1^3 - Xn*Yn*Ynp1^2 - 2*Xn*Yn*Ynp1 - Xn*Ynp1^3 + 2*Xn*Ynp1^2 - Xnp1^3*Yn^2 + 2*Xnp1^3*Yn + Xnp1*Yn^3*Ynp1 - Xnp1*Yn^3 - Xnp1*Yn^2*Ynp1^2 - Xnp1*Yn^2*Ynp1 + 2*Xnp1*Yn^2 + 2*Xnp1*Yn*Ynp1^2 - 2*Xnp1*Yn*Ynp1 - Yn*term + Ynp1*term)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)), (Rn^3*Rnp1*Ynp1 - Rn^3*Rnp1 - Rn^3*Ynp1 + Rn^3 - Rn^2*Rnp1^2*Yn - Rn^2*Rnp1^2*Ynp1 + 2*Rn^2*Rnp1^2 + 2*Rn^2*Rnp1*Yn - Rn^2*Rnp1*Ynp1 - Rn^2*Rnp1 - Rn^2*Xn*Xnp1*Ynp1 + Rn^2*Xn*Xnp1 + Rn^2*Xnp1^2*Yn + Rn^2*Xnp1^2*Ynp1 - Rn^2*Yn + Rn^2*Ynp1^3 - Rn^2*Ynp1^2 + Rn^2*Ynp1 + Rn*Rnp1^3*Yn - Rn*Rnp1^3 - Rn*Rnp1^2*Yn + 2*Rn*Rnp1^2*Ynp1 - Rn*Rnp1^2 - Rn*Rnp1*Xn^2*Ynp1 + Rn*Rnp1*Xn^2 - 4*Rn*Rnp1*Xn*Xnp1 - Rn*Rnp1*Xnp1^2*Yn + Rn*Rnp1*Xnp1^2 - Rn*Rnp1*Yn^2*Ynp1 + Rn*Rnp1*Yn^2 - Rn*Rnp1*Yn*Ynp1^2 + Rn*Rnp1*Ynp1^2 + Rn*Xn^2*Ynp1 - Rn*Xn^2 + 2*Rn*Xn*Xnp1*Ynp1 + 2*Rn*Xn*Xnp1 - Rn*Xnp1^2*Yn - 2*Rn*Xnp1^2*Ynp1 - Rn*Xnp1^2 + Rn*Xnp1*term + Rn*Yn^2*Ynp1 - Rn*Yn^2 + Rn*Yn*Ynp1^2 - 2*Rn*Ynp1^3 + Rn*Ynp1^2 - Rnp1^3*Yn + Rnp1^3 + Rnp1^2*Xn^2*Yn + Rnp1^2*Xn^2*Ynp1 - Rnp1^2*Xn*Xnp1*Yn + Rnp1^2*Xn*Xnp1 + Rnp1^2*Yn^3 - Rnp1^2*Yn^2 + Rnp1^2*Yn - Rnp1^2*Ynp1 - 2*Rnp1*Xn^2*Yn - Rnp1*Xn^2*Ynp1 - Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1*Yn + 2*Rnp1*Xn*Xnp1 - Rnp1*Xn*term + Rnp1*Xnp1^2*Yn - Rnp1*Xnp1^2 - 2*Rnp1*Yn^3 + Rnp1*Yn^2*Ynp1 + Rnp1*Yn^2 + Rnp1*Yn*Ynp1^2 - Rnp1*Ynp1^2 + Xn^3*Xnp1*Ynp1 - Xn^3*Xnp1 - Xn^2*Xnp1^2*Yn - Xn^2*Xnp1^2*Ynp1 + 2*Xn^2*Xnp1^2 + Xn^2*Yn - Xn^2*Ynp1^3 + Xn^2*Ynp1^2 + Xn^2*Ynp1 + Xn*Xnp1^3*Yn - Xn*Xnp1^3 + Xn*Xnp1*Yn^2*Ynp1 - Xn*Xnp1*Yn^2 + Xn*Xnp1*Yn*Ynp1^2 - 2*Xn*Xnp1*Yn - Xn*Xnp1*Ynp1^2 - 2*Xn*Xnp1*Ynp1 + Xn*term - Xnp1^2*Yn^3 + Xnp1^2*Yn^2 + Xnp1^2*Yn + Xnp1^2*Ynp1 - Xnp1*term + Yn^3 - Yn^2*Ynp1 - Yn*Ynp1^2 + Ynp1^3)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)), sqrt(-(Rn^2 - 2*Rn - Xn^2 - Yn^2 + 2*Yn)*(Rnp1^2 - 2*Rnp1 - Xnp1^2 - Ynp1^2 + 2*Ynp1)*(Rn^2 - 2*Rn*Rnp1 + Rnp1^2 - Xn^2 + 2*Xn*Xnp1 - Xnp1^2 - Yn^2 + 2*Yn*Ynp1 - Ynp1^2))*(Xn*Ynp1 - Xn - Xnp1*Yn + Xnp1)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)) - (Rn^3*Xnp1^2 + Rn^3*Ynp1^2 - 2*Rn^3*Ynp1 + Rn^3 - Rn^2*Rnp1*Xn*Xnp1 - Rn^2*Rnp1*Yn*Ynp1 + Rn^2*Rnp1*Yn + Rn^2*Rnp1*Ynp1 - Rn^2*Rnp1 + Rn^2*Xn*Xnp1 - Rn^2*Xnp1^2 + Rn^2*Yn*Ynp1 - Rn^2*Yn - Rn^2*Ynp1^2 + Rn^2*Ynp1 - Rn*Rnp1^2*Xn*Xnp1 - Rn*Rnp1^2*Yn*Ynp1 + Rn*Rnp1^2*Yn + Rn*Rnp1^2*Ynp1 - Rn*Rnp1^2 - Rn*Xn^2*Xnp1^2 - Rn*Xn^2*Ynp1^2 + 2*Rn*Xn^2*Ynp1 - Rn*Xn^2 + Rn*Xn*Xnp1^3 + Rn*Xn*Xnp1*Ynp1^2 - 2*Rn*Xn*Xnp1*Ynp1 + 2*Rn*Xn*Xnp1 - Rn*Xnp1^2*Yn^2 + Rn*Xnp1^2*Yn*Ynp1 + Rn*Xnp1^2*Yn - Rn*Xnp1^2*Ynp1 - Rn*Xnp1^2 - Rn*Yn^2*Ynp1^2 + 2*Rn*Yn^2*Ynp1 - Rn*Yn^2 + Rn*Yn*Ynp1^3 - Rn*Yn*Ynp1^2 - Rn*Ynp1^3 + Rn*Ynp1^2 + Rnp1^3*Xn^2 + Rnp1^3*Yn^2 - 2*Rnp1^3*Yn + Rnp1^3 - Rnp1^2*Xn^2 + Rnp1^2*Xn*Xnp1 - Rnp1^2*Yn^2 + Rnp1^2*Yn*Ynp1 + Rnp1^2*Yn - Rnp1^2*Ynp1 + Rnp1*Xn^3*Xnp1 - Rnp1*Xn^2*Xnp1^2 + Rnp1*Xn^2*Yn*Ynp1 - Rnp1*Xn^2*Yn - Rnp1*Xn^2*Ynp1^2 + Rnp1*Xn^2*Ynp1 - Rnp1*Xn^2 + Rnp1*Xn*Xnp1*Yn^2 - 2*Rnp1*Xn*Xnp1*Yn + 2*Rnp1*Xn*Xnp1 - Rnp1*Xnp1^2*Yn^2 + 2*Rnp1*Xnp1^2*Yn - Rnp1*Xnp1^2 + Rnp1*Yn^3*Ynp1 - Rnp1*Yn^3 - Rnp1*Yn^2*Ynp1^2 - Rnp1*Yn^2*Ynp1 + Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1^2 - Rnp1*Ynp1^2 - Xn^3*Xnp1 + 2*Xn^2*Xnp1^2 - Xn^2*Yn*Ynp1 + Xn^2*Yn + Xn^2*Ynp1 - Xn*Xnp1^3 - Xn*Xnp1*Yn^2 + 4*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - Xn*Xnp1*Ynp1^2 - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn*Ynp1 + Xnp1^2*Yn + Xnp1^2*Ynp1 - Yn^3*Ynp1 + Yn^3 + 2*Yn^2*Ynp1^2 - Yn^2*Ynp1 - Yn*Ynp1^3 - Yn*Ynp1^2 + Ynp1^3)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)));
end
乙 (1.33333, 2.00000, 0.66667)
丙 (1.33333, 1.00000, 0.33333)
丁 (1.09091, 0.54545, 0.18182)
戊 (0.88889, 0.33333, 0.11111)
己 (0.74074, 0.22222, 0.07407)
子 (0.53333, 2.00000, 0.13333) 元の単位で 5.06667 寸
丑 (1.04348, 1.30435, 0.08696) 元の単位で 3.30435 寸
寅 (1.02564, 0.76923, 0.05128) 元の単位で 1.94872 寸
卯 (0.88889, 0.47619, 0.03175) 元の単位で 1.20635 寸
辰 (0.75789, 0.31579, 0.02105) 元の単位で 0.80000 寸
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 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")
R0 = 2 # 大円の半径
R1 = R0 // 2 # 甲円の半径
X2, Y2, R2 = (4//3, 2//1, 2//3)
XnA = zeros(6)
YnA = zeros(6)
RnA = zeros(6)
plot()
circle(0, R0, R0, :black)
circle(0, R1, R1, :blue)
circle(0, 3R1, R1, :green)
circle(X2, Y2, R2, :green)
circle(-X2, Y2, R2, :green)
(XnA[1], YnA[1], RnA[1]) = (Xn, Yn, Rn) = (X2, Y2, R2)
large = Char["乙丙丁戊己 "...]
for i = 1:5
@printf("%s (%.5f, %.5f, %.5f)\n", large[i], Xn, Yn, Rn)
circle(Xn, Yn, Rn, :green)
circle(-Xn, Yn, Rn, :green)
more && point(Xn, Yn, " " * large[i], mark=false)
i == 5 && break
(XnA[i+1], YnA[i+1], RnA[i+1]) = (Xn, Yn, Rn) = NEXTCIRCLE(Xn, Yn, Rn)
end
(xn, yn, rn) = x2.evalf(), y2.evalf(), r2.evalf()
circle(xn, yn, rn, :blue, fill=true)
circle(-xn, yn, rn, :blue, fill=true)
more && point(xn, yn, "子 ", :blue, :right, mark=false)
small = Char["子丑寅卯辰"...]
@printf("%s (%.5f, %.5f, %.5f) 元の単位で %.5f 寸\n", small[1], xn, yn, rn, 2rn * 76 / 4)
for i = 1:4
(xn, yn, rn) = nextcircle(XnA[i], YnA[i], RnA[i], XnA[i+1], YnA[i+1], RnA[i+1])
@printf("%s (%.5f, %.5f, %.5f) 元の単位で %.5f 寸\n", small[i+1], xn, yn, rn, 2rn * 76 / 4)
circle(xn, yn, rn, :blue, fill=true)
circle(-xn, yn, rn, :blue, fill=true)
more && i < 5 && point(xn, yn, small[i+1] * " ", :blue, :right, :bottom, mark=false)
end
if more == true
point(0, 0, "O ", :green, :right, :bottom)
point(0, R1, "R1 ", :green, :right, :bottom)
point(0, R0, "R0=2R1", :green, :center, :bottom)
point(0, R0+R1, "3R1 ", :green, :right, :bottom)
point(0, 2R0, "2R0 ", :green, :right, :bottom)
point(0, 3R1, " 甲", :green, mark=false)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
