SymPy で算額の問題を解く

算額とは？

いろいろな算額を SymPy で解いています。

SymPy は Pyhon で使うのが一般的とは思いますが，Julia の Python で書いています。

Julia なので，Python の SymPy とはちょっと書き方にくせがありますが，かえって書きやすい面もあります。

ご意見，ご感想いただければうれしいです。

よろしければ，「いいね！」してください。

神壁算法 關龍藤田貞資門人 龜井隠岐守家士 堀田人助泉尹 天明八年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

直径 97 寸 5 分の大円の中に累円と，大円と甲円と累円に挟まれる円がある。
末円の直径が 1 分のとき，初円から末円まで何個の挟円があるか（つまり，末円は初円から数えて何番目か）。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (x1, y1); x1 = 0, y1 = r1
乙円の半径と中心座標を r2, (x2, y2); x2 = r0 - r2, y2 = 0
乙円以降 i 番目の円の半径と中心座標を ri, (xi, yi)
初円の半径と中心座標を R1, X1, Y1
初円以降 i 番目の挟円の半径と中心座標を Ri, (Xi, Yi)
とおき，順次パラメータを決定していく。

from sympy import *

var('r0, r1, r2, x2, y2, r3, x3, y3')

(r0, r1, r2, x2, y2, r3, x3, y3)

まず，乙円については，半径 r2 = r0/3 になる。

r1 = r0/2
x2 = r0 - r2
y2 = 0
eq1 = Eq(x2**2 + r1**2, (r1 + r2)**2)
res2 = solve(eq1, r2)
print(res2[0])

r0/3

次に，丙円については，半径，中心座標は (r0/6, 2*r0/3, r0/2) になる。

r2 = r0/3
x2 = r0 - r2
y2 = 0
eq2 = Eq(x3**2 + y3**2, (r0 - r3)**2)
eq3 = Eq(x3**2 + (r1 - y3)**2, (r1 + r3)**2)
eq4 = Eq((x2 - x3)**2 + (y2 - y3)**2, (r2 + r3)**2)
res3 = solve([eq2, eq3, eq4], (r3, x3, y3))
print(res3)

[(r0/6, 2*r0/3, r0/2), (r0/2, 0, -r0/2)]

更に，丁円については，半径，中心座標は (r0/11, 6r0/11, 8r0/11) になる。

var('r4, x4, y4')
r3 = r0/6
x3 = 2*r0/3
y3 = r0/2
eq5 = Eq(x4**2 + y4**2, (r0 - r4)**2)
eq6 = Eq(x4**2 + (r1 - y4)**2, (r1 + r4)**2)
eq7 = Eq((x3 - x4)**2 + (y3 - y4)**2, (r3 + r4)**2)
res4 = solve([eq5, eq6, eq7], (r4, x4, y4))
res4

[(r0/11, 6*r0/11, 8*r0/11), (r0/3, 2*r0/3, 0)]

var('r5, x5, y5')
r4 = r0/11
x4 = 6*r0/11
y4 = 8*r0/11
eq8 = Eq(x5**2 + y5**2, (r0 - r5)**2)
eq9 = Eq(x5**2 + (r1 - y5)**2, (r1 + r5)**2)
eq10 = Eq((x4 - x5)**2 + (y4 - y5)**2, (r4 + r5)**2)
res5 = solve([eq8, eq9, eq10], (r5, x5, y5))
res5

[(r0/18, 4*r0/9, 5*r0/6), (r0/6, 2*r0/3, r0/2)]

var('r6, x6, y6')
(r5, x5, y5) = res5[0]
eq11 = Eq(x6**2 + y6**2, (r0 - r6)**2)
eq12 = Eq(x6**2 + (r1 - y6)**2, (r1 + r6)**2)
eq13 = Eq((x5 - x6)**2 + (y5 - y6)**2, (r5 + r6)**2)
res6 = solve([eq11, eq12, eq13], (r6, x6, y6))
res6

[(r0/27, 10*r0/27, 8*r0/9), (r0/11, 6*r0/11, 8*r0/11)]

var('r7, x7, y7')
(r6, x6, y6) = res6[1]
eq14 = Eq(x7**2 + y7**2,(r0 - r7)**2)
eq15 = Eq(x7**2 + (r1 - y7)**2, (r1 + r7)**2)
eq16 = Eq((x6 - x7)**2 + (y6 - y7)**2, (r6 + r7)**2)
res7= solve([eq14, eq15, eq16], (r7, x7, y7))
res7

[(r0/18, 4*r0/9, 5*r0/6), (r0/6, 2*r0/3, r0/2)]

var('r8, x8, y8')
(r7, x7, y7) = res7[1]
eq17 = Eq(x8**2 + y8**2, (r0 - r8)**2)
eq18 = Eq(x8**2 + (r1 - y8)**2, (r1 + r8)**2)
eq19 = Eq((x7 - x8)**2 + (y7 - y8)**2, (r7 + r8)**2)
res8= solve([eq17, eq18, eq19], (r8, x8, y8))
res8

[(r0/11, 6*r0/11, 8*r0/11), (r0/3, 2*r0/3, 0)]

累円のパラメータ，半径，中心座標を表す数列には規則性があり，以下のような一般項になる。
甲円の半径 r1，中心座標 (x1, y1) は r0/2, (0, r1)
乙円の半径 r2，中心座標 (x2, y2) は r0/3, (r0 - r2, 0)
丙円以降は i ≧ 3 として，ri, (xi, yi) は以下のようになる。
ri = r0/(i2 - 2i + 3)
xi = ri*(2i - 2)
yi = ri*(i2 - 2i)

次に，互いに接する，外円と 2 個の累円に挟まれる円のパラメータを計算する。

var('R1, X1')
X1 = r0 - 2*r2 - R1
eq20 = Eq(X1**2 + r1**2, (r1 + R1)**2)
res11 = solve(eq20, R1)
print(res11[0])

r0/15

R1 = r0/15
X1 = r0 - 2*r2 - R1
(R1, X1)

(r0/15, 4*r0/15)

var('R2, X2, Y2')
(x1, y1, x2, y2) = (0, r1, r0 - r2, 0)
(R1, X1) = (r0/15, 4*r0/15)
eq21 = Eq( X2**2       + (Y2 - y1)**2, (R2 + r1)**2)
eq22 = Eq((X2 - x2)**2 + (Y2 - y2)**2, (R2 + r2)**2)
eq23 = Eq((X2 - x3)**2 + (Y2 - y3)**2, (R2 + r3)**2)
res12 = solve([eq21, eq22, eq23], (R2, X2, Y2))
res12

[(-r0, 0, 0), (r0/23, 12*r0/23, 8*r0/23)]

var('R3, X3, Y3')
(R2, X2, Y2) = res12[1]
eq24 = Eq( X3**2       + (Y3 - y1)**2, (R3 + r1)**2)
eq25 = Eq((X3 - x3)**2 + (Y3 - y3)**2, (R3 + r3)**2)
eq26 = Eq((X3 - x4)**2 + (Y3 - y4)**2, (R3 + r4)**2)
res13 = solve([eq24, eq25, eq26], (R3, X3, Y3))
res13

[(-r0, 0, 0), (r0/39, 20*r0/39, 8*r0/13)]

var('R4, X4, Y4')
(R3, X3, Y3) = res13[1]
eq27 = Eq( X4**2       + (Y4 - y1)**2, (R4 + r1)**2)
eq28 = Eq((X4 - x4)**2 + (Y4 - y4)**2, (R4 + r4)**2)
eq29 = Eq((X4 - x5)**2 + (Y4 - y5)**2, (R4 + r5)**2)
res14 = solve([eq27, eq28, eq29], (R4, X4, Y4))
res14

[(-r0, 0, 0), (r0/63, 4*r0/9, 16*r0/21)]

var('R5, X5, Y5')
(R4, X4, Y4) = res14[1]
eq30 = Eq( X5**2       + (Y5 - y1)**2, (R5 + r1)**2)
eq31 = Eq((X5 - x5)**2 + (Y5 - y5)**2, (R5 + r5)**2)
eq32 = Eq((X5 - x6)**2 + (Y5 - y6)**2, (R5 + r6)**2)
res14 = solve([eq30, eq31, eq32], (R5, X5, Y5))
res14

[(-r0, 0, 0), (r0/63, 4*r0/9, 16*r0/21)]

挟円のパラメータは
初円は
R1 = r0/15, X1 = 4*r0/15, Y1 = 0
二円以降 i ≧ 2 においては
半径 Ri = r0/(4i2 - 4i + 15)
x座標 Xi = Ri * (8i - 4)
y座標 Yi = Ri * (4i2 - 4i)

問は，「末円の径が一分であるとき，それは何番目の挟円か」ということである。
以下の方程式を解けば，16番目の挟円であることがわかる。

var('i')
r0 = 975/20  # 直径 = 97.5
eq = Eq(r0/(4*i**2 - 4*i + 15), 1/20)  # 直径 = 0.1 寸
print(solve(eq, i)[1])

16.0000000000000

一般的な解を求めるには，末円の半径を Rn とすれば

var('i, r0, Rn')
eq = r0/(4*i**2 - 4*i + 15) - Rn
res = solve(eq, i)[1]
print(res)

(Rn + sqrt(Rn*(-14*Rn + r0)))/(2*Rn)

SymPy ではこれ以上簡単にできないが，手計算で，(1 + sqrt(r0/Rn - 14))/2 となる。

「外円の半径を末円の半径で割り，14を引いたものを開平し1を加えて2で割る」となり，「術」と一致する。
なお，当たり前であるが，術は，半径を直径と読み替えても成り立つ。

ちなみに，20番目の挟円の半径は 0.031759 ほどであり，数値を丸めるとたしかに20番目である事がわかる。

(sqrt(97.5/(2*0.031759) - 14) + 1)/2

$\displaystyle 19.999986880355$