算額（その1545）
埼玉県川口市三ツ和 氷川神社 享和4年(1804)
山口正義：やまぶき 第 39 号
https://yamabukiwasan.sakura.ne.jp/ymbk39.pdf
キーワード：円5個，直角三角形
#Julia, #SymPy, #算額, #和算
直角三角形の中に，甲，乙，丙，丁，戊の 5 円を容れる。直角三角形の股（底辺）と勾配が与えられたとき，各辺の長さと各円の直径を求めよ。
図形としては算額（その206）と同じである。何を与えて何を求めるかが違うだけである。
注：勾配といえば Δy/Δx であるが，当時はその 1/10 の値を使っていたのだろうか，また，勾配の有効桁が 4 桁しかないので，本来は整数値を意図していたはずのものに端数が出てしまう。そこで，鈎は問題文の「高倍（勾配）四寸弐分八厘五毛」を使わずに，3 とする。
鈎，股，弦をそのまま，「鈎」，「股」，「弦」
甲円の半径と中心座標を r1, (股 - r1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
戊円の半径と中心座標を r5, (x5, r5)
とおき，順次解く。
長いので，先に答えを書いておく。
鈎 = 3; 股 = 7; 弦 = 7.61577
直径： 甲円 = 2.38423; 乙円 = 1.58596; 丙円 = 1.05496; 丁円 = 0.701746; 戊円 = 0.466793
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms 鈎, 股, 弦, r1, r2, x2, r3, x3, r4, x4, r5, x5
弦 = sqrt(鈎^2 + 股^2);
# r1:甲円の半径
eq1 = (鈎 + 股 - 弦) - 2r1
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
ans_r1(鈎 => 3, 股 => 7).evalf() |> println
1.19211344706805
@syms x1, x2
r1 = ans_r1 # 1.19211344706805
x1 = 股 - r1
eq2 = (x1 - x2) - 2sqrt(r1*r2)
eq3 = r1/x1 - r2/x2;
res = solve([eq2, eq3], r2, x2)[1];
# r2:乙円の半径
ans_r2 = res[1]
ans_r2 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))^2
ans_r2(鈎 => 3, 股 => 7).evalf() |> println
0.792979085478491
# x2:乙円の中心の x 座標
ans_x2 = res[2]
ans_x2 |> println
(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
ans_x2(鈎 => 3, 股 => 7).evalf() |> println
3.86333413034970
累円の半径は，公比 p = r2/r1 = 0.665187602261171 の 等比数列をなす。
# p:公比
p = res[1]/ans_r1
p |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^2)
p(鈎 => 3, 股 => 7).evalf() |> println
0.665187602261170
# r3:丙円の半径
r3 = ans_r2*p
r3 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))^2*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^2/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^4)
r3(鈎 => 3, 股 => 7).evalf() |> println
0.527479856512693
# r4:丁円の半径
r4 = ans_r2*p^2
r4 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^2*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)
r4(鈎 => 3, 股 => 7).evalf() |> println
0.350873060994744
# r5:戊円の半径
r5 = ans_r2*p^3
r5 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))^4*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^4/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^8)
r5(鈎 => 3, 股 => 7).evalf() |> println
0.233396410141131
互いに接する円の半径が既知ならば，中心間の水平距離は d12 = 2sqrt(r1*r2) であるが，算額の場合は，距離も等比級数に従う。
# d23:乙円と丙円の中心の水平距離
d23 = 2sqrt(ans_r2*r3)
d23 |> println
2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6))
d23(鈎 => 3, 股 => 7).evalf() |> println
1.29349216344864
# d34:丙円と丁円の中心の水平距離
d34 = 2sqrt(r3*r4)
d34 |> println
2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^5*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^5/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^10))
d34(鈎 => 3, 股 => 7).evalf() |> println
0.860414950748014
# d45:丁円と戊円の中心の水平距離
d45 = 2sqrt(r4*r5)
d45 |> println
2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^7*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^7/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^5*(股 - 鈎 + sqrt(股^2 + 鈎^2))^14))
d45(鈎 => 3, 股 => 7).evalf() |> println
0.572337358037734
# x3:丙円の中心の x 座標
x2 = res[2]
x3 = x2 - d23
x3 |> println
-2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)) + (3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
x3(鈎 => 3, 股 => 7).evalf() |> println
2.56984196690106
# x4:丁円の中心の x 座標
x4 = x3 - d34
x4 |> println
-2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^5*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^5/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^10)) - 2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)) + (3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
x4(鈎 => 3, 股 => 7).evalf() |> println
1.70942701615304
# x5:戊円の中心の x 座標
x5 = x4 - d45
x5 |> println
-2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^7*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^7/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^5*(股 - 鈎 + sqrt(股^2 + 鈎^2))^14)) - 2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^5*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^5/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^10)) - 2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)) + (3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
x5(鈎 => 3, 股 => 7).evalf() |> println
1.13708965811531
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(鈎, 股) = (3, 7)
弦 = sqrt(鈎^2 + 股^2)
r1 = 1.19211344706805
(r2, x2) = (0.792979085478493, 3.86333413034969)
(r3, r4, r5) = (0.527479856512695, 0.350873060994746, 0.233396410141133)
(x3, x4, x5) = (2.56984196690104, 1.70942701615302, 1.13708965811529)
@printf("鈎 = %g; 股 = %g; 弦 = %g\n", 鈎, 股, 弦)
@printf("直径： 甲円 = %g; 乙円 = %g; 丙円 = %g; 丁円 = %g; 戊円 = %g\n", 2r1, 2r2, 2r3, 2r4, 2r5)
plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:red, lw=0.5)
circle(股 - r1, r1, r1, 1)
circle(x2, r2, r2, 2)
circle(x3, r3, r3, 3)
circle(x4, r4, r4, 4)
circle(x5, r5, r5, 5)
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)
end
end;
draw(true)
