算額（その1143）

算額（その1143）

二九 戸田市美女木 八幡社 天保4年(1833)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：直角三角形，正方形,繰り返し

直角三角形の中に甲，乙，丙，丁の正方形を 4 個容れる。甲と丁の正方形の一辺の長さが与えられたとき，丙の正方形の一辺の長さはいかほどか。

直角三角形の直角を挟む二辺の短い方を鈎，長い方を股とする。
甲，乙，丙，丁の正方形の一辺の長さを a, b, c, d とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms 鈎::positive, 股::positive,
     a::positive, b::positive, c::positive,　d::positive,
     e::positive, f::positive
eq1 = a/(股 - a) - 鈎/股
eq2 = b/(股 - a - b) - 鈎/股
eq3 = c/(股 - a - b - c) - 鈎/股
eq4 = d/(股 - a - b - c - d) - 鈎/股
res4 = solve([eq1, eq2, eq3, eq4], (鈎, 股, b, c))[1]

   (a^(4/3)/d^(1/3), a^(7/3)/(a^(4/3) - a*d^(1/3)), a^(2/3)*d^(1/3), a^(1/3)*d^(2/3))

res4[1](a => 4, d => 1).evalf() |> println
res4[2](a => 4, d => 1).evalf() |> println
res4[3](a => 4, d => 1).evalf() |> println
res4[4](a => 4, d => 1).evalf() |> println

   6.34960420787280
   10.8096575356773
   2.51984209978975
   1.58740105196820

正方形が増えるとどのような規則性があるのかやってみる。

eq5 = e/(股 - a - b - c - d - e) - 鈎/股
eq6 = f/(股 - a - b - c - d - e - f) - 鈎/股

# n = 5
res5 = solve([eq1, eq2, eq3, eq4, eq5], (鈎, 股, b, c, d))[1]

   (a^(5/4)/e^(1/4), a^(9/4)/(a^(5/4) - a*e^(1/4)), a^(3/4)*e^(1/4), sqrt(a)*sqrt(e), a^(1/4)*e^(3/4))

res5[1](a => 4, e => 1).evalf() |> println
res5[2](a => 4, e => 1).evalf() |> println
res5[3](a => 4, e => 1).evalf() |> println
res5[4](a => 4, e => 1).evalf() |> println
res5[5](a => 4, e => 1).evalf() |> println

   5.65685424949238
   13.6568542494924
   2.82842712474619
   2.00000000000000
   1.41421356237310

# n = 6
res6 = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鈎, 股, b, c, d, e))[1]

   (a^(6/5)/f^(1/5), a^(11/5)/(a^(6/5) - a*f^(1/5)), a^(4/5)*f^(1/5), a^(3/5)*f^(2/5), a^(2/5)*f^(3/5), a^(1/5)*f^(4/5))

res6[1](a => 4, f => 1).evalf() |> println
res6[2](a => 4, f => 1).evalf() |> println
res6[3](a => 4, f => 1).evalf() |> println
res6[4](a => 4, f => 1).evalf() |> println
res6[5](a => 4, f => 1).evalf() |> println
res6[6](a => 4, f => 1).evalf() |> println

   5.27803164309158
   16.5192518405067
   3.03143313302080
   2.29739670999407
   1.74110112659225
   1.31950791077289

正方形の数を n，最初と最後の正方形の一辺の長さを a，z とすれば，以下のようになる。
鈎 = a^(n/(n - 1))/z^(1/(n - 1))
股 = a^((2n - 1)/(n - 1)) - a*z^(1/(n - 1))
i 番目の正方形の一辺の長さは
a^((n - i)/(n - 1))*z^(i/(n - 1)), i = 1,...,n
である。

n = 6; a = 4; z = 1
println("鈎 = ", a^(n/(n - 1))/z^(1/(n - 1)))
println("股 = ", a^((2n - 1)/(n - 1))/(a^(n/(n - 1)) - a*z^(1/(n - 1))))
for i = 1:n
   println("i = $i; ", a^((n - i)/(n - 1))*z^(i/(n - 1)))
end

   鈎 = 5.278031643091577
   股 = 16.519251840506687
   i = 1; 4.0
   i = 2; 3.0314331330207964
   i = 3; 2.2973967099940698
   i = 4; 1.7411011265922482
   i = 5; 1.3195079107728942
   i = 6; 1.0

function draw(n, a, z, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   鈎 = a^(n/(n - 1))/z^(1/(n - 1))
   股 = a^((2n - 1)/(n - 1))/(a^(n/(n - 1)) - a*z^(1/(n - 1)))
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   prev = 0
   for i = 1:n
       side = a^((n - i)/(n - 1))*z^(i/(n - 1))
       @printf("%d 個目の正方形の一辺の長さは %g\n", i, side)
       #plot!([prev, prev + loc[i], prev + loc[i]], [loc[i], loc[i], 0], color=:red, lw=0.5)
       plot!([prev, prev + side, prev + side], [side, side, 0], color=i, lw=0.5)
       prev += side
   end
   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(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta/2)
       point(股, 0, " 股", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その1142）

算額（その1142）

二五 武州妻沼 聖天宮 文政11年(1828)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：累円，外円

大円の中に等円 3 個を容れ，甲円から始まり，乙円，丙円，...，己円，... を容れる。
等円の直径が 3.5 寸のとき，己円の直径はいかほどか。

等円の半径と中心座標を r, (0, r), (0, 0), (0, -r)
大円の半径と中心座標を 3r, (0, 0)
中円の半径と中心座標を 2r, (0, r)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
:
とおき以下の連立方程式を解く。

まず，甲円について，r1, x1, y1 を求める。

include("julia-source.txt")

using SymPy

@syms r::positive, r1::positive, x1::positive,　y1::negative
r = 3.5/2
eq1 = x1^2 + y1^2 - (3r - r1)^2
eq2 = x1^2 + (r - y1)^2 - (2r + r1)^2
eq3 = x1^2 + (y1 + 2r)^2 - (r + r1)^2
res = solve([eq1, eq2, eq3], (r1, x1, y1))[1]

   (1.50000000000000, 3.00000000000000, -2.25000000000000)

大円，中円，等円，甲円の内接・外接関係は大円，中円，甲円，丙円と同じである。
丙円の半径，中心位置を求める連立方程式は先程の連立方程式の変数を変更して解く。

まずは計算を float 変数で行う方法について述べる。

@syms r::real, r1::real, x1::real, y1::real, r2::real, x2::real, y2::real
#(r1, x1, y1) = (arg_r1, arg_x1, arg_y1)
eq4 = x2^2 + y2^2 - (3r - r2)^2
eq5 = x2^2 + (r - y2)^2 - (2r + r2)^2
eq6 = (x2 - x1)^2 + (y2 - y1)^2 - (r1 + r2)^2
res = solve([eq4, eq5, eq6], (r2, x2, y2))[1]

   ((27*r^3 + 9*r^2*r1 - 63*r^2*y1 - 3*r*r1^2 - 6*r*r1*y1 + 27*r*x1^2 + 33*r*y1^2 - r1^3 + 5*r1^2*y1 + r1*x1^2 + r1*y1^2 - 5*x1^2*y1 - 2*sqrt(6)*x1*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 5*y1^3)/(2*(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2)), (72*r^2*x1 - 12*r*r1*x1 - 12*r*x1*y1 + 3*sqrt(6)*r*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 12*r1^2*x1 + sqrt(6)*r1*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) + 12*x1^3 + 12*x1*y1^2 - 5*sqrt(6)*y1*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4))/(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2), 5*sqrt(6)*x1*sqrt(-(-9*r^2 - 6*r*r1 - r1^2 + x1^2 + y1^2)*(-3*r^2 + 4*r*r1 - 2*r*y1 - r1^2 + x1^2 + y1^2))/(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2) + (-81*r^3 - 9*r^2*r1 + 135*r^2*y1 + 21*r*r1^2 - 30*r*r1*y1 + 9*r*x1^2 - 15*r*y1^2 + 5*r1^3 - 25*r1^2*y1 - 5*r1*x1^2 - 5*r1*y1^2 + 25*x1^2*y1 + 25*y1^3)/(2*(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2)))

これを関数として定義する。関数は，現在の円の半径と中心位置を与えれば次の円の半径と中心位置を返す。返された結果を関数に与えれば，更に次の円の半径と中心位置を返す。これを繰り返す。

この方法の利点は，関数内部の変数は普通の float 変数なので，計算時間が短い。
欠点は，関数定義が長い。

r = 3.5/2
nextcircle2(r1, x1, y1) = (((27*r^3 + 9*r^2*r1 - 63*r^2*y1 - 3*r*r1^2 - 6*r*r1*y1 + 27*r*x1^2 + 33*r*y1^2 - r1^3 + 5*r1^2*y1 + r1*x1^2 + r1*y1^2 - 5*x1^2*y1 - 2*sqrt(6)*x1*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 5*y1^3)/(2*(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2)), (72*r^2*x1 - 12*r*r1*x1 - 12*r*x1*y1 + 3*sqrt(6)*r*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) - 12*r1^2*x1 + sqrt(6)*r1*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4) + 12*x1^3 + 12*x1*y1^2 - 5*sqrt(6)*y1*sqrt(-27*r^4 + 18*r^3*r1 - 18*r^3*y1 + 12*r^2*r1^2 - 12*r^2*r1*y1 + 12*r^2*x1^2 + 12*r^2*y1^2 - 2*r*r1^3 - 2*r*r1^2*y1 + 2*r*r1*x1^2 + 2*r*r1*y1^2 + 2*r*x1^2*y1 + 2*r*y1^3 - r1^4 + 2*r1^2*x1^2 + 2*r1^2*y1^2 - x1^4 - 2*x1^2*y1^2 - y1^4))/(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2), 5*sqrt(6)*x1*sqrt(-(-9*r^2 - 6*r*r1 - r1^2 + x1^2 + y1^2)*(-3*r^2 + 4*r*r1 - 2*r*y1 - r1^2 + x1^2 + y1^2))/(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2) + (-81*r^3 - 9*r^2*r1 + 135*r^2*y1 + 21*r*r1^2 - 30*r*r1*y1 + 9*r*x1^2 - 15*r*y1^2 + 5*r1^3 - 25*r1^2*y1 - 5*r1*x1^2 - 5*r1*y1^2 + 25*x1^2*y1 + 25*y1^3)/(2*(9*r^2 + 6*r*r1 - 30*r*y1 + r1^2 - 10*r1*y1 + 24*x1^2 + 25*y1^2))));

nextcircle2(1.5, 3.0, -2.25)  # 乙円を求める

   (1.05, 4.2, 2.220446049250313e-16)

nextcircle2(1.05, 4.2, 0)  # 丙円を求める

   (0.7, 4.2, 1.75)

nextcircle2(0.7, 4.2, 1.75)  # 丁円を求める

   (0.47727272727272724, 3.818181818181818, 2.8636363636363633)

nextcircle2(0.47727272727272724, 3.818181818181818, 2.8636363636363633)  # 戊円を求める

   (0.3387096774193549, 3.3870967741935485, 3.5564516129032255)

nextcircle2(0.3387096774193549, 3.3870967741935485, 3.5564516129032255)  # 己円を求める

   (0.25000000000000006, 2.9999999999999996, 3.9999999999999982)

次の方法は，計算を SymPy 変数で行う方法である。毎回連立方程式を解くので，計算時間がかかる場合がある。

function nextcircle(arg_r1, arg_x1, arg_y1)
   @syms r1::real, x1::real, y1::real, r2::real, x2::real, y2::real
   (r1, x1, y1) = (arg_r1, arg_x1, arg_y1)
   eq4 = x2^2 + y2^2 - (3r - r2)^2
   eq5 = x2^2 + (r - y2)^2 - (2r + r2)^2
   eq6 = (x2 - x1)^2 + (y2 - y1)^2 - (r1 + r2)^2
   res = solve([eq4, eq5, eq6], (r2, x2, y2))[1]
end;

nextcircle(1.5, 3.0, -2.25)  # 乙円を求める

   (1.05000000000000, 4.20000000000000, 0.0)

nextcircle(1.05000000000000, 4.20000000000000, 0.0)  # 丙円を求める

   (0.700000000000000, 4.20000000000000, 1.75000000000000)

nextcircle(0.700000000000000, 4.20000000000000, 1.75000000000000)  # 丁円を求める

   (0.477272727272727, 3.81818181818182, 2.86363636363636)

nextcircle(0.477272727272727, 3.81818181818182, 2.86363636363636)  # 戊円を求める

   (0.338709677419356, 3.38709677419355, 3.55645161290322)

nextcircle(0.338709677419356, 3.38709677419355, 3.55645161290322) # 己円を求める

   (0.250000000000001, 3.00000000000000, 4.00000000000000)

1.　デカルトの円定理を用いる方法

円の直径（半径）だけを求めるならば，デカルトの円定理を使って，累円の半径を求めてゆくことができる。
半径 r1, prev, nxt の 3 円が内接する外円の半径を R とすると，以下の式が成り立つ。

using SymPy
@syms R, r1, prev, nxt
nxtr = -1/(1/r1 + 1/prev + 1/nxt - 2sqrt(1/(r1*prev) + 1/(prev*nxt) + 1/(nxt*r1))) - R;

これを解いて nxt を求める式を得る。

result = solve(nxtr, nxt)[1]
result |> println

   (R*prev*r1*(R*prev + R*r1 - prev*r1) - 2*sqrt(-R^3*prev^3*r1^3*(-R + prev + r1)))/(R^2*prev^2 - 2*R^2*prev*r1 + R^2*r1^2 + 2*R*prev^2*r1 + 2*R*prev*r1^2 + prev^2*r1^2)

これを関数にする。この関数は，円の半径を入力すれば，その円に外接する次の円の半径を返す。

res2(prev, R, r1) = (R*prev*r1*(R*prev + R*r1 - prev*r1) - 2*sqrt(-R^3*prev^3*r1^3*(-R + prev + r1)))/(R^2*prev^2 - 2*R^2*prev*r1 + R^2*r1^2 + 2*R*prev^2*r1 + 2*R*prev*r1^2 + prev^2*r1^2);

甲円が求まった段階で，次々に関数を適用してゆく。

using SymPy

@syms r::positive, r1::positive, x1::positive,　y1::negative
r = 3.5/2
eq1 = x1^2 + y1^2 - (3r - r1)^2
eq2 = x1^2 + (r - y1)^2 - (2r + r1)^2
eq3 = x1^2 + (y1 + 2r)^2 - (r + r1)^2
res = solve([eq1, eq2, eq3], (r1, x1, y1))[1]

   (1.50000000000000, 3.00000000000000, -2.25000000000000)

r1 = 1.50000000000000  # 甲円

   1.5

nxt = res2(r1, 3r, 2r)  # 乙円

   1.05

nxt = res2(1.05, 3r, 2r)  # 丙円

   0.6999999999999997

nxt = res2(0.6999999999999997, 3r, 2r)  # 丁円

   0.47727272727272696

nxt = res2(0.47727272727272696, 3r, 2r)  # 戊円

   0.33870967741935454

nxt = res2(0.33870967741935454, 3r, 2r)  # 己円

   0.24999999999999983

nxt = res2(0.24999999999999983, 3r, 2r)  # 庚円

   0.19090909090909078

nxt = res2(0.19090909090909078, 3r, 2r)  # 辛円

   0.14999999999999988

nxt = res2(0.14999999999999988, 3r, 2r)  # 壬円

   0.12068965517241372

nxt = res2(0.12068965517241372, 3r, 2r)  # 癸円

   0.09905660377358488

各円の半径を分数表示すると，分子はどれでも 21，分母は二階等差数列である。

rationalize(1.5) |> println
rationalize(1.05) |> println
rationalize(0.699999999999999, tol=1e-8) |> println
rationalize(0.47727272727272, tol=1e-8) |> println
rationalize(0.33870967741935454, tol=1e-8) |> println
rationalize(0.25, tol=1e-8) |> println
rationalize(0.19090909090909078, tol=1e-8) |> println
rationalize(0.14999999999999988, tol=1e-8) |> println
rationalize(0.12068965517241372, tol=1e-8) |> println
rationalize(0.09905660377358488, tol=1e-8) |> println

   3//2
   21//20
   7//10
   21//44
   21//62
   1//4
   21//110
   3//20
   7//58
   21//212

甲円を 1 番目の円，乙円を 2 番目の円と数えるとき，n 番目の円の半径は 21/(12 + 2n^2) である。
己円は 6 番目の円であるから，その半径は　21/(12 + 2*6^2) = 0.25 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   string = Char["甲乙丙丁戊己庚辛壬癸"...]
   plot(showaxis=false)
   circlef(0, 0, 3r)
   circlef(0, r, 2r, :blue)
   circlef(0, 2r, r, :magenta)
   circlef(0, 0, r, :magenta)
   circlef(0, -2r, r, :magenta)
   (r1, x1, y1) = (1.5, 3.0, -2.25)
   for i = 1:10
       @printf("%2d: %.15f, %.15f, % .15f, %3d, %.15f\n", i, r1, x1, y1, 12 + 2i^2, 21/(12 + 2i^2))
       circlef(x1, y1, r1, i)
       circlef(-x1, y1, r1, i)
       i < 7 && point(x1, y1, string[i], :white, :center, :vcenter, mark=false)
       (r1, x1, y1) = nextcircle2(r1, x1, y1)
   end
   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(0, -2r, "等円", :white, :center, :vcenter, mark=false)
   end
end;

算額（その1141）

算額（その1141）

二五 武州妻沼 聖天宮 文政11年(1828)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，外円，弦

外円の中に水平な弦を引き，大円，中円，小円を容れる。大円，中円の直径がそれぞれ 10 寸，2.5 寸のとき，小円の直径はいかほどか。

弦と y 軸の交点座標を (0, y); y = 2r2 - R
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, y + r1)
中円の半径と中心座標を r2, (x2, y + r2), (0, y - r2)
小円の半径と中心座標を r3, (0, y + r3)
とおき以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms R::positive, y::positive, r1::positive, x1::positive,
     r2::positive, x2::negative, r3::positive, x3::positive
y = 2r2 - R
eq1 = x1^2 + (y + r1)^2 - (R - r1)^2
eq2 = x2^2 + (y + r2)^2 - (R - r2)^2
eq3 = x3^2 + (y + r3)^2 - (R - r3)^2
eq4 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2;
(R, x1, x2, r3, x3) = solve([eq1, eq2, eq3, eq4, eq5], (R, x1, x2, r3, x3))[1];

R  |> println
x1 |> factor |> println
x2 |> factor |> println
r3 |> println
x3 |> factor |> println

   (2*r1 + r2)^2/(4*r1)
   r2^(3/2)/sqrt(r1)
   -sqrt(r2)*(2*r1 - r2)/sqrt(r1)
   r1^2*r2/(r1 + r2)^2
   sqrt(r2)*(2*r1^2 + r1*r2 + r2^2)/(sqrt(r1)*(r1 + r2))

小円の半径 r3 は，大円，中円の半径 r1, r2 により r1^2*r2/(r1 + r2)^2 倍である。
大円，中円の直径がそれぞれ 10 寸，2.5 寸のとき，小円の直径は 1.6 寸である。

その他のパラメータは以下の通りである。

   R = 6.32812; r1 = 5; x1 = 0.625;  r2 = 1.25;  x2 = -4.375;  r3 = 0.8;  x3 = 4.625

function draw(r1, r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R  = (2r1 + r2)^2/4r1
   x1 = r2^(3/2)/√r1
   x2 = √r2*(r2 - 2r1)/√r1
   r3 = r1^2*r2/(r1 + r2)^2
   x3 = √r2*(2r1^2 + r1*r2 + r2^2)/(√r1*(r1 + r2))
   @printf("大円，中円の直径が %g，%g のとき，小円の直径は %g である。\n", 2r1, 2r2, 2r3)
   @printf("R = %g; r1 = %g; x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", R, r1, x1, r2, x2, r3, x3)
   y = 2r2 - R
   plot()
   circle(0, 0, R)
   circle(x1, y + r1, r1, :magenta)
   circle(0, y - r2, r2, :blue)
   circle(x2, y + r2, r2, :blue)
   circle(x3, y + r3, r3, :green)
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :orange)
   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(0, 0, "", :red)
       point(0, R, "R", :red, :center, :bottom, delta=delta)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, y, "y", :black, :center, :bottom, delta=delta)
       point(x1, y + r1, "大円:r1,(x1,y+r1)", :magenta, :center, :bottom, delta=delta/2)
       point(x2, y + r2, "中円:r2\n(x2,y+r2)", :blue, :center, :bottom, delta=delta/2)
       point(0, y - r2, "中円:r2\n(0,y-r2)", :blue, :center, delta=-delta/2)
       point(x3, y + r3, "小円:r3,(x3,y+r3)", :black, :right, :vcenter, deltax=-delta)
   end
end;