精要算法
和算研究所，佐藤健一：和算百科，p.200, 丸善出版株式会社，2017/10/10.
キーワード：円4個，三角形
#Julia #SymPy #算額 #和算 #数学
三角形の中に，全円，大円，中円，小円が内接している。大円の直径が 9 寸，中円の直径が 4 寸，小円の直径が 1 寸のとき，全円の直径はいくらか。
以下の論文の著者からコメントを貰った。
全円，大円，中円，小円の半径を \(R, r_1, r_2, r_3\) とすると，\(R = \sqrt{r_1 r_2} + \sqrt{r_2 r_3} + \sqrt{r_3 r_1}\) であることが，以下の論文で示されている。
HIROSHI OKUMURA:THREE CIRCLE CHAINS ARISING FROM THREE LINES, Global Journal of Advanced Research
on Classical and Modern Geometries
ISSN: 2284-5569, Vol.9, (2020), Issue 1, pp.1-7
以下は，小生の解。
図のように記号を定める。
大円，中円，小円の半径を \(r_1, r_2, r_3\)
大円のある方の三角形の頂点から全円と大円との接点までの距離を \(b, b_1\)
中円のある方の三角形の頂点から全円と大円との接点までの距離を \(a, a_1\)
小円のある方の三角形の頂点から全円と大円との接点までの距離を \(c, c_1\)
小円の中心座標を \( (x_3, y_3)\)
三角形の頂点を \( (0, 0), (a + b, 0), (x, y)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive,
a::positive, b::positive, c::positive,
a1::positive, b1::positive, c1::positive,
r1::positive, r2::positive, r3::positive
a1 = a - 2sqrt(R*r2)
b1 = b - 2sqrt(R*r1)
c1 = c - 2sqrt(R*r3)
eq1 = R*a1 - r2*a # R/a - r2/a1
eq2 = R*b1 - r1*b # R/b - r1/b1
eq3 = R*c1 - r3*c # R/c - r3/c1
eq4 = R^2*(a + b + c) - a*b*c; # 算法助術の公式36
\(r_1, r_2, r_3\) が与えられたとき，\(a, b, c\) は eq1, eq2, eq3 を解くことで求められる（\(R\) の関数になる）。
ans_a = solve(eq1, a)[1]
@show(ans_a)
ans_a = 2*R^(3/2)*sqrt(r2)/(R - r2)
\(\displaystyle \frac{2 R^{\frac{3}{2}} \sqrt{r_{2}}}{R - r_{2}}\)
ans_b = solve(eq2, b)[1]
@show(ans_b)
ans_b = 2*R^(3/2)*sqrt(r1)/(R - r1)
\(\displaystyle \frac{2 R^{\frac{3}{2}} \sqrt{r_{1}}}{R - r_{1}}\)
ans_c = solve(eq3, c)[1]
@show(ans_c)
ans_c = 2*R^(3/2)*sqrt(r3)/(R - r3)
\(\displaystyle \frac{2 R^{\frac{3}{2}} \sqrt{r_{3}}}{R - r_{3}}\)
\(a, b, c\) を eq4 に代入し，\(R\) を求める。
eq40 = eq4(a => ans_a, b => ans_b, c => ans_c)
\(\displaystyle - \frac{8 R^{\frac{9}{2}} \sqrt{r_{1}} \sqrt{r_{2}} \sqrt{r_{3}}}{\left(R - r_{1}\right) \left(R - r_{2}\right) \left(R - r_{3}\right)} + R^{2} \left(\frac{2 R^{\frac{3}{2}} \sqrt{r_{1}}}{R - r_{1}} + \frac{2 R^{\frac{3}{2}} \sqrt{r_{2}}}{R - r_{2}} + \frac{2 R^{\frac{3}{2}} \sqrt{r_{3}}}{R - r_{3}}\right)\)
ans_R = solve(eq40, R)[2]
\(\displaystyle \frac{4 \sqrt{r_{1}} \sqrt{r_{2}} \sqrt{r_{3}} + \sqrt{r_{1}} r_{2} + \sqrt{r_{1}} r_{3} + r_{1} \sqrt{r_{2}} + r_{1} \sqrt{r_{3}} + \sqrt{r_{2}} r_{3} + r_{2} \sqrt{r_{3}} + \sqrt{2 r_{1}^{\frac{3}{2}} r_{2}^{\frac{3}{2}} + 6 r_{1}^{\frac{3}{2}} \sqrt{r_{2}} r_{3} + 6 r_{1}^{\frac{3}{2}} r_{2} \sqrt{r_{3}} + 2 r_{1}^{\frac{3}{2}} r_{3}^{\frac{3}{2}} + 6 \sqrt{r_{1}} r_{2}^{\frac{3}{2}} r_{3} + 2 \sqrt{r_{1}} \sqrt{r_{2}} r_{3}^{2} + 2 \sqrt{r_{1}} r_{2}^{2} \sqrt{r_{3}} + 6 \sqrt{r_{1}} r_{2} r_{3}^{\frac{3}{2}} + 2 r_{1}^{2} \sqrt{r_{2}} \sqrt{r_{3}} + r_{1}^{2} r_{2} + r_{1}^{2} r_{3} + 6 r_{1} r_{2}^{\frac{3}{2}} \sqrt{r_{3}} + 6 r_{1} \sqrt{r_{2}} r_{3}^{\frac{3}{2}} + r_{1} r_{2}^{2} + 10 r_{1} r_{2} r_{3} + r_{1} r_{3}^{2} + 2 r_{2}^{\frac{3}{2}} r_{3}^{\frac{3}{2}} + r_{2}^{2} r_{3} + r_{2} r_{3}^{2}}}{2 \sqrt{r_{1}} + 2 \sqrt{r_{2}} + 2 \sqrt{r_{3}}}\)
「問」では「大円の直径が 9 寸，中円の直径が 4 寸，小円の直径が 1 寸」なので，半径を代入すれば，全円の半径が 11/2 寸，直径は 11 寸となり，「答」に一致する。
ans_R(r1 => 9//2, r2 => 4//2, r3 => 1//2) |> println
11/2
答えがあっているので，安心して式 \(ans_R\) を簡約化することに専心することができる。
まず式中の \(\sqrt{r_1}, \sqrt{r_2}, \sqrt{r_3}\) を \(s, t, u\) とする。
@syms s, t, u
ans_R2 = ans_R(√r1 => s, √r2 => t, √r3 => u)
\(\displaystyle \frac{s^{2} t + s^{2} u + s t^{2} + 4 s t u + s u^{2} + t^{2} u + t u^{2} + \sqrt{s^{4} t^{2} + 2 s^{4} t u + s^{4} u^{2} + 2 s^{3} t^{3} + 6 s^{3} t^{2} u + 6 s^{3} t u^{2} + 2 s^{3} u^{3} + s^{2} t^{4} + 6 s^{2} t^{3} u + 10 s^{2} t^{2} u^{2} + 6 s^{2} t u^{3} + s^{2} u^{4} + 2 s t^{4} u + 6 s t^{3} u^{2} + 6 s t^{2} u^{3} + 2 s t u^{4} + t^{4} u^{2} + 2 t^{3} u^{3} + t^{2} u^{4}}}{2 s + 2 t + 2 u}\)
\(\LaTeX\) 式は長くなり，本ブログでははみ出すので，テキストで表示しておく。
ans_R2 |> println
(s^2*t + s^2*u + s*t^2 + 4*s*t*u + s*u^2 + t^2*u + t*u^2 + sqrt(s^4*t^2 + 2*s^4*t*u + s^4*u^2 + 2*s^3*t^3 + 6*s^3*t^2*u + 6*s^3*t*u^2 + 2*s^3*u^3 + s^2*t^4 + 6*s^2*t^3*u + 10*s^2*t^2*u^2 + 6*s^2*t*u^3 + s^2*u^4 + 2*s*t^4*u + 6*s*t^3*u^2 + 6*s*t^2*u^3 + 2*s*t*u^4 + t^4*u^2 + 2*t^3*u^3 + t^2*u^4))/(2*s + 2*t + 2*u)
分子のルートの中身を \(v\) とおく。
v = s^4*t^2 + 2*s^4*t*u + s^4*u^2 + 2*s^3*t^3 + 6*s^3*t^2*u + 6*s^3*t*u^2 + 2*s^3*u^3 + s^2*t^4 + 6*s^2*t^3*u + 10*s^2*t^2*u^2 + 6*s^2*t*u^3 + s^2*u^4 + 2*s*t^4*u + 6*s*t^3*u^2 + 6*s*t^2*u^3 + 2*s*t*u^4 + t^4*u^2 + 2*t^3*u^3 + t^2*u^4
\(\displaystyle s^{4} t^{2} + 2 s^{4} t u + s^{4} u^{2} + 2 s^{3} t^{3} + 6 s^{3} t^{2} u + 6 s^{3} t u^{2} + 2 s^{3} u^{3} + s^{2} t^{4} + 6 s^{2} t^{3} u + 10 s^{2} t^{2} u^{2} + 6 s^{2} t u^{3} + s^{2} u^{4} + 2 s t^{4} u + 6 s t^{3} u^{2} + 6 s t^{2} u^{3} + 2 s t u^{4} + t^{4} u^{2} + 2 t^{3} u^{3} + t^{2} u^{4}\)
\(v\) は因子分解すると以下のように簡単になり，\(s + t > 0, s + u > 0, t + u > 0\) なので，平方根の部分は \( ( (s + t)(s + u)(t + u)\) になる。
v |> factor
\(\displaystyle \left(s + t\right)^{2} \left(s + u\right)^{2} \left(t + u\right)^{2}\)
以上により，分子は以下のようになる。
x = s^2*t + s^2*u + s*t^2 + 4*s*t*u + s*u^2 + t^2*u + t*u^2 + (s + t)*(s + u)*(t + u)
\(\displaystyle s^{2} t + s^{2} u + s t^{2} + 4 s t u + s u^{2} + t^{2} u + t u^{2} + \left(s + t\right) \left(s + u\right) \left(t + u\right)\)
これを因数分解すると更に簡単な式になる。
x |> factor
\(\displaystyle 2 \left(s + t + u\right) \left(s t + s u + t u\right)\)
分母が \(2(s + t + u)\) なので，結局，全円の半径は \(st + su + tu\)である。
これを \(s, t, u\) を元の \(\sqrt{r_1}, \sqrt{r_2}, \sqrt{r_3}\) に戻すと以下のようになる。
ans_R2 = s*t + s*u + t*u
\(\displaystyle s t + s u + t u\)
ans_R3 = ans_R2(s => √r1, t => √r2, u => √r3)
@show(ans_R3)
ans_R3 = sqrt(r1)*sqrt(r2) + sqrt(r1)*sqrt(r3) + sqrt(r2)*sqrt(r3)
\(\displaystyle \sqrt{r_{1}} \sqrt{r_{2}} + \sqrt{r_{1}} \sqrt{r_{3}} + \sqrt{r_{2}} \sqrt{r_{3}}\)
全円の半径は \(\sqrt{r_1 r_2} + \sqrt{r_2 r_3} + \sqrt{r_3 r_1}\) である。
正しく数式変換できたかどうか，「大円の直径が 9 寸，中円の直径が 4 寸，小円の直径が 1 寸」なので，半径を代入」し，全円の半径が 11/2 寸になることを確認できる。
ans_R3(r1 => 9//2, r2 => 4//2, r3 => 1//2) |> println
11/2
「術」は以下のようになっている。
@syms 大径, 中径, 小径
寄位 = sqrt(中径*小径)
全径 = sqrt( (2*2 + 中径 + 小径)*大径) + 寄位
\(\displaystyle \sqrt{中径\cdot 小径} + \sqrt{大径 \left(中径 + 小径 + 4\right)}\)
全径は \(\sqrt{中径\cdot 小径} + \sqrt{大径\cdot(中径 + 小径 + 4)}\) となる。
大径，中径，小径をさまざまに変えて全径を得るために関数にしてみる。
function 全径f(大径, 中径, 小径)
寄位 = sqrt(中径*小径)
sqrt( (2*2 + 中径 + 小径)*大径) + 寄位
end;
「大円の直径が 9 寸，中円の直径が 4 寸，小円の直径が 1 寸」のとき，全円の直径は 11 寸となり，「答」に一致する。
全径f(9, 4, 1)
11.0
得られる三角形は高さのわりに底辺が長いもので，図に描くと見栄えが悪い（原因は大円の直径と中円，小円の直径が違いすぎる）。
なので，それぞれの直径が 9, 8, 7 寸のときの全円の直径を求めてみる。ページの前の方にある図はこの場合のものである。
全円の直径は 20.560011604169905 寸ということになる。
全径f(9, 8, 7)
20.560011604169905
しかし，上で求めた式 \(ans_R\) によれば（半径を与えて，直径を得ることに注意），23.9058500809802 寸になる。
ans_R3(r1 => 9//2, r2 => 8//2, r3 => 7//2).evalf() * 2 |> println
23.9058500809802
解析解は大円，中円，小円がどのような値でも正しい結果が与えられるはずである。
ここで，「術」の式を見てみると不自然な点がある。「大径, 中径, 小径は大きさの順ではあるが，図形としては 3 個の円は固定されるものではない，つまり，結果の式は大径, 中径, 小径について交代式でなければならない」ということである。\(ans_R\) はこの条件を満たしている。
また，術式の中に「2×2」という部分があることも不可思議である。なぜ「4」としないのか。
術式は，大径 = 9, 中径 = 4, 小径 = 1 を使って，全径 = 11 という結果が出るように，適当に細工をしたのではないかと疑われる。
今まで，佐藤健一を含め，だれもこの術式がおかしいことに気づかなかったのであろうか。
なお，問題文の通りの図は平べったい三角形になるが，それは，「問と答の中に出てくる数は整数」という縛りがあるため，虱潰しで探索した結果 9, 4, 1 11 という結果になったのであろう。
コンピュータを使って虱潰しで候補を探索すると面白い（？）ことがわかる。全円の半径の式に \(\sqrt{r_1}, \sqrt{r_2}, \sqrt{r_3}\) があることからもわかるが，それぞれの直径は平方数（および定数倍）である必要があるということである。
大径 = 100, 中径 = 81, 小径 = 64 を使って，全径 = 242 などを得るには根気が続かなかったのであろう。
全径f(100, 81, 64)
194.06555615733703 これは間違い！！
ans_R3(r1 => 100//2, r2 => 81//2, r3 => 64//2).evalf() * 2 |> println
242.000000000000
描画関数プログラムのソースを見る
function get_x3_y3(a, b, x, y, r3)
@syms x3, y3
eq1 = dist2(0, 0, x, y, x3, y3, r3)
eq2 = dist2(a + b, 0, x, y, x3, y3, r3)
res = solve([eq1, eq2], (x3, y3))
return (real(res[2][1]), real(res[2][2]))
end;
function draw(r1, r2, r3, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = sqrt(r1*r2) + sqrt(r1*r3) + sqrt(r2*r3)
string = @sprintf("r1 = %g\nr2 = %g\nr3 = %g\n R = %g", r1, r2, r3, R)
println(string)
a = 2R^(3/2)*√r2/(R - r2)
b = 2R^(3/2)*√r1/(R - r1)
c = 2R^(3/2)*√r3/(R - r3)
a1 = a - 2sqrt(R*r2)
b1 = b - 2sqrt(R*r1)
c1 = c - 2sqrt(R*r3)
plot()
circle(b1, r1, r1)
circle(b, R, R, :blue)
circle(b + a - a1, r2, r2, :green)
θ1 = atand(R, b)
θ2 = atand(R, a)
(x, y) = Float64.(intersection(0, 0, b, b*tand(2θ1), b + a, 0, b, -tand(2θ2)*(b - (b + a))))
plot!([0, b + a, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
(x3, y3) = Float64.(get_x3_y3(a, b, x, y, r3))
circle(x3, y3, r3, :magenta)
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
#point(x, y)
point(b, R, "大円:R,(b,R)", :blue, :center, :bottom, delta=delta/2)
point(b1, r1, "大円:r1\n(b1,r1)", :red, :center, :bottom, delta=delta/2)
point(b + a - a1, r2, "中円:r2\n(b+a-a1,r2)", :green, :center, :bottom, delta=delta/2)
point(x3, y3, "小円:r3,(x3,y3)", :magenta, :left, :vcenter, deltax=10delta)
point( (a + b)/10, y3, string, :black, :left, mark=false)
end
end;
draw(9/2, 8/2, 7/2, true)
draw(100, 81, 64, true)
draw(9/2, 4/2, 1/2, false)
以下のアイコンをクリックして応援してください
