算額（その611）

算額（その611）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/10

等脚台形内に，上円，中円，下円，甲円，乙円と斜線を入れる。中円の直径が 4 寸，子（中円の直径から甲円と乙円の直径の和を引いた長さ）が 1 寸のとき，台形の高さはいかほどか。

台形の下底，上底，高さを a, b, h とする。また，斜線と斜辺の交点座標を (x3, (a - x3)*h/(a - b)), (x4, (a - x4)*h/(a - b)) とする。
上円の半径と中心座標を r1, (0, h - r1)
中円の半径と中心座標を r2, (0, 2r3 + r2)
下円の半径と中心座標を r3, (0, r3)
甲円の半径と中心座標を r4, (0, h - 2r1 - r4)
乙円の半径と中心座標を r5, (0, 2r3 + r5)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive, x3::positive, x4::positive,
     r1::positive, r2::positive, r3::positive, r4::positive, r5::positive, 子::positive
eq1 = distance(a, 0, b, h, 0, h - r1) - r1^2
eq2 = distance(a, 0, b, h, 0, 2r3 + r2) - r2^2
eq3 = distance(a, 0, b, h, 0, r3) - r3^2
eq4 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, h - r1) - r1^2
eq5 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, r3) - r3^2
eq6 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, h - 2r1 - r4) - r4^2
eq7 = distance(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b), 0, 2r3 + r5) - r5^2
eq8 = 2r2 - 2r4 - 2r5 - 子
eq9 = 2r1 + 2r2 + 2r3 - h;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, b, h, x3, x4, r1, r3, r4, r5) = u
   return [
       h^2*(a*r1 + b*h - b*r1)^2/(a^2 - 2*a*b + b^2 + h^2)^2 - r1^2 + (h - h*(a^2 - a*b + h^2 - h*r1)/(a^2 - 2*a*b + b^2 + h^2) - r1)^2,  # eq1
       h^2*(a*h - a*r2 - 2*a*r3 + b*r2 + 2*b*r3)^2/(a^2 - 2*a*b + b^2 + h^2)^2 - r2^2 + (-h*(a^2 - a*b + h*r2 + 2*h*r3)/(a^2 - 2*a*b + b^2 + h^2) + r2 + 2*r3)^2,  # eq2
       h^2*(a*h - a*r3 + b*r3)^2/(a^2 - 2*a*b + b^2 + h^2)^2 - r3^2 + (-h*(a^2 - a*b + h*r3)/(a^2 - 2*a*b + b^2 + h^2) + r3)^2,  # eq3
       h^2*(a*r1*x3^2 - a*r1*x4^2 + b*h*x3^2 - b*h*x4^2 - b*r1*x3^2 + b*r1*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r1^2 + (h - h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2 - h*r1*x3^2 + 2*h*r1*x3*x4 - h*r1*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) - r1)^2,  # eq4
       h^2*(a*h*x3^2 - a*h*x4^2 - a*r3*x3^2 + a*r3*x4^2 + b*r3*x3^2 - b*r3*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r3^2 + (-h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + h*r3*x3^2 - 2*h*r3*x3*x4 + h*r3*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) + r3)^2,  # eq5
       h^2*(2*a*r1*x3^2 - 2*a*r1*x4^2 + a*r4*x3^2 - a*r4*x4^2 + b*h*x3^2 - b*h*x4^2 - 2*b*r1*x3^2 + 2*b*r1*x4^2 - b*r4*x3^2 + b*r4*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r4^2 + (h - h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2 - 2*h*r1*x3^2 + 4*h*r1*x3*x4 - 2*h*r1*x4^2 - h*r4*x3^2 + 2*h*r4*x3*x4 - h*r4*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) - 2*r1 - r4)^2,  # eq6
       h^2*(a*h*x3^2 - a*h*x4^2 - 2*a*r3*x3^2 + 2*a*r3*x4^2 - a*r5*x3^2 + a*r5*x4^2 + 2*b*r3*x3^2 - 2*b*r3*x4^2 + b*r5*x3^2 - b*r5*x4^2 - 2*h*x3^2*x4 + 2*h*x3*x4^2)^2/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2)^2 - r5^2 + (-h*(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - a*b*x3^2 - 2*a*b*x3*x4 - a*b*x4^2 - 2*a*x3^2*x4 - 2*a*x3*x4^2 + 2*b*x3^2*x4 + 2*b*x3*x4^2 + 2*h*r3*x3^2 - 4*h*r3*x3*x4 + 2*h*r3*x4^2 + h*r5*x3^2 - 2*h*r5*x3*x4 + h*r5*x4^2)/(a^2*x3^2 + 2*a^2*x3*x4 + a^2*x4^2 - 2*a*b*x3^2 - 4*a*b*x3*x4 - 2*a*b*x4^2 + b^2*x3^2 + 2*b^2*x3*x4 + b^2*x4^2 + h^2*x3^2 - 2*h^2*x3*x4 + h^2*x4^2) + 2*r3 + r5)^2,  # eq7
       2*r2 - 2*r4 - 2*r5 - 子,  # eq8
       -h + 2*r1 + 2*r2 + 2*r3,  # eq9
   ]
end;

r2 = 4//2
子 = 1
iniv = BigFloat[8.5, 0.5, 16.0, 4.5, 0.9, 0.8, 5.2, 0.2, 1.3]
res = nls(H, ini=iniv)

   (BigFloat[8.472135954999579392818347337462552470881319131313934445299371062594802345966918, 0.472135954999579392818347337462552470881231767596567280477615683449053318730822, 16.00000000000000000000000000000000000000005685400937294517590029875627505751969, 4.47213595499957939281834733746255247088125870170645007332837022256876548756112, 0.8944271909999158785636694674925104941762430878938323936182682049502660273589708, 0.7639320225002103035908263312687237645593763638467272007538116532706423285987476, 5.236067977499789696409173668731276235440652063157959271834138496107495199876993, 0.1909830056250525758977065828171809411398420533761430006087923002379663556053446, 1.309016994374947424102293417182819058860157946623856999391207699762033644413767], true)

中円の直径が 4 寸，子が 1 寸のとき，等脚台形の高さは 16 寸である。

その他のパラメータは以下の通りである。
r2 = 2;  子 = 1;  a = 8.47214;  b = 0.472136;  h = 16;  x3 = 4.47214;  x4 = 0.894427;  r1 = 0.763932;  r3 = 5.23607;  r4 = 0.190983;  r5 1.30902

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, h, x3, x4, r1, r3, r4, r5) = res[1]
   @printf("高さ = %g;  r2 = %g;  子 = %g;  a = %g;  b = %g;  h = %g;  x3 = %g;  x4 = %g;  r1 = %g;  r3 = %g;  r4 = %g;  r5 %g\n",
       h, r2, 子, a, b, h, x3, x4, r1, r3, r4, r5)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:gray60, lw=0.5)
   circle(0, h - r1, r1, :blue)
   circle(0, 2r3 + r2, r2, :magenta)
   circle(0, r3, r3, :green)
   circle(0, h - 2r1 - r4, r4, :brown)
   circle(0, 2r3 + r5, r5, :red)
   segment(x3, (a - x3)*h/(a - b), -x4, (a - x4)*h/(a - b))
   segment(-x3, (a - x3)*h/(a - b), x4, (a - x4)*h/(a - b))
   segment(0, 2r3 + 2r5, 0, h - 2r1 - 2r4, lw=3)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, h - r1, "上円:r1,(0,h-r1)    ", :blue, :right, :vcenter)
       point(0, h - 2r1 - r4, "     甲円:r4,(0,h-2r1-r4)", :brown, :left, :vcenter)
       point(0, 2r3 + r2, "          中円:r2,(0,2r3+r2)", :magenta, :left, :vcenter)
       point(0, 2r3 + r5, "          乙円:r5,(0,2r3+r5)", :red, :left, :vcenter)
       point(0, r3, " 下円:r3,(0,r3)", :green, :left, :vcenter)
       point(x3, (a - x3)*h/(a - b), " x3", :black, :left, :vcenter)
       point(x4, (a - x4)*h/(a - b), " x4", :black, :left, :vcenter)
       point(a, 0, "a", :black, :left, :bottom, delta=delta/2)
       point(b, h, "(b,h)", :black, :left, :bottom, delta=delta/2)
       point(0, 2r3 + 1.5r2, "子  ", :black, :right, :vcenter, mark=false)
   end
end;

算額（その610）

算額（その610）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/10

外円内に弦を引き，等円を 4 個配置する。
外円の直径が与えられたとき等円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r), (x, y); y = 3r - R
とおき，連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive, y::negative
y = 2r - R + r
eq1 = x^2 + y^2 - (R - r)^2
eq2 = x^2 + (R - r - y)^2 - 4r^2
res = solve([eq1, eq2], (r, x));

res[1][1] |> simplify |> println
res[1][2] |> simplify |> println

   R*(3 - sqrt(5))/2
   R*sqrt(-22 + 10*sqrt(5))

等円の半径は外円の半径の (3 - √5)/2 倍である。
例えば，外円の直径が 987 のとき，等円の直径は 377.0004531038537 である。

987*(3 - √5)/2

   377.0004531038537

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

   等円の直径 = 754.001;  x = 592.759;  y = 144.001;  R = 987

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 987
   (r, x) = R .* ((3 - sqrt(5))/2, sqrt(-22 + 10*sqrt(5)))
   y = 2r - R + r
   @printf("等円の直径 = %g;  x = %g;  y = %g;  R = %g\n", 2r, x, y, R)
   plot()
   circle(0, 0, R)
   circle(0, R - r, r, :blue)
   circle(0, r - R, r, :blue)
   circle(x, y, r, :blue)
   circle(-x, y, r, :blue)
   x0 = sqrt(R^2 - (2r - R)^2)
   segment(-x0, 2r - R, x0, 2r - R, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - r, " R-r", :blue, :left, :vcenter)
       point(0, r - R, " r-R", :blue, :left, :vcenter)
       point(0, 2r - R, " 2r-R", :blue, :left, :bottom, delta=delta/2)
       point(x, y, "等円:r,(x,y)", :blue, :center, :top, delta=-delta)
   end
end;

算額（その609）

算額（その609）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/10

直角三角形内に，甲円，乙円，丙円，丁円を入れる。甲円と乙円の直径がわかっているとき丙円の直径を求めよ。

直角三角形の直角を挟む二辺を鈎，股とする（鈎 < 股）
甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (y2, y2)
丙円の半径と中心座標を r3, (r3, y3)
丁円の半径と中心座標を r4, (r4, y4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, r1::positive, y1, r2::negative, y2::positive, r3::positive, y3::negative, r4::positive, y4::positive
(x2, x3, x4) = (y2, y3, y4)
eq1 = (x2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq2 = (r1 - r4)^2 + (y4 - r1)^2 - (r1 + r4)^2
eq3 = (x2 - r3)^2 + (y3 - y2)^2 - (r2 + r3)^2
eq4 = (x2 - r4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq5 = (r3 - r4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq6 = distance(股, 0, 0, 鈎, r3, y3) - r3^2
eq7 = distance(股, 0, 0, 鈎, x2, y2) - r2^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (鈎, 股, y2, r3, y3, r4, y4) = u
   return [
       2*(-r1 + y2)^2 - (r1 + r2)^2,  # eq1
       (-r1 + y4)^2 + (r1 - r4)^2 - (r1 + r4)^2,  # eq2
       -(r2 + r3)^2 + (-r3 + y2)^2 + (-y2 + y3)^2,  # eq3
       -(r2 + r4)^2 + (-r4 + y2)^2 + (y2 - y4)^2,  # eq4
       (r3 - r4)^2 - (r3 + r4)^2 + (y3 - y4)^2,  # eq5
       -r3^2 + (r3 - 股*(r3*股 - y3*鈎 + 鈎^2)/(股^2 + 鈎^2))^2 + (y3 - 鈎*(-r3*股 + y3*鈎 + 股^2)/(股^2 + 鈎^2))^2,  # eq6
       -r2^2 + (y2 - 股*(y2*股 - y2*鈎 + 鈎^2)/(股^2 + 鈎^2))^2 + (y2 - 鈎*(-y2*股 + y2*鈎 + 股^2)/(股^2 + 鈎^2))^2,  # eq7
   ]
end;

(r1, r2) = (9, 12)
iniv = BigFloat[75.6, 62.1, 25.6, 11.3, 44.0, 6.8, 26.5]
res = nls(H, ini=iniv)

   (BigFloat[54.28853657836170682638366694115126793584060670779405752143209928308113720353747, 80.96484995197442957569524070254618022769191543098055398929484432041730484231544, 23.84924240491749801241773160420182982498155469145795476835513724890269102385213, 8.705582264797570515028006208872174856682679273246349361204917648446301735851336, 37.96981939330018192322196742329539101348115699766049412475605009936624145822072, 5.925450840795221324838051774317623074624056457781674516833436921472992954425592, 23.60534937167296822394913519268628446794816052457474295677168184443295330276947], true)

   鈎 = 54.2885;  股 = 80.9648;  r1 = 9;  r2 = 12;  y2 = 23.8492;  r3 = 8.70558;  y3 = 37.9698;  r4 = 5.92545;  y4 = 23.6053

甲円と乙円の直径が 18 寸，24 寸のとき丙円の直径は 17.4112 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, y2, r3, y3, r4, y4) = res[1]
   (x2, x3, x4) = (y2, y3, y4)
   @printf("丙円の直径 = %g\n", 2r3)
   @printf("鈎 = %g;  股 = %g;  r1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  y3 = %g;  r4 = %g;  y4 = %g\n",
       鈎, 股, r1, r2, y2, r3, y3, r4, y4)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:brown, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(r3, y3, r3, :green)
   circle(x3, r3, r3, :green)
   circle(r4, y4, r4, :magenta)
   circle(x4, r4, r4, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r1, r1, "甲円:r1\n(r1,r1)", :red, :center, delta=-delta)
       point(y2, y2, "乙円:r2\n(y2,y2)", :blue, :center, delta=-delta)
       point(r3, y3, "丙円:r3\n(r3,y3)", :green, :center, delta=-delta)
       point(r4, y4, "丁円:r4\n(r4,y4)", :magenta, :center, :vcenter)
       point(股, 0, "股", :brown, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :brown, :left, :bottom)
   end
end;

算額（その608）

算額（その608）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/11

半円内に2本の斜線を引き，区画された領域に甲円，乙円，丙円を入れる。乙円，丙円の直径が与えられたとき，甲円の直径を求めよ。

算額（その511）より一段階複雑になっている。

半円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
2 本の斜線が半円の円周上で交差する座標を (x, y)
とおき以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, x::negative, y::positive, r1::positive, x1::negative, r2::positive, x2::positive, y2::positive, r3::positive, x3::negative, y3::positive

y = sqrt(R^2 - x^2)
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = distance(R, 0, x, y, x1, r1) - r1^2
eq4 = y/(x + R) * y3/x3 + 1
eq5 = y/(x - R) * y2/x2 + 1
eq6 = (sqrt((x + R)^2 + y^2) + sqrt((R - x)^2 + y^2) + 2R)*r1 - 2R*y
eq7 = distance(R, 0, x, y, x2, y2) - r2^2
eq8 = distance(-R, 0, x, y, x3, y3) - r3^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (R, x, r1, x1, x2, y2, x3, y3) = u
   return [
       x2^2 + y2^2 - (R - r2)^2,  # eq1
       x3^2 + y3^2 - (R - r3)^2,  # eq2
       -r1^2 + (r1 - (R*r1 + R*sqrt(R^2 - x^2) + r1*x - x1*sqrt(R^2 - x^2))/(2*R))^2 + (x1 - (R^2 + R*x + R*x1 - r1*sqrt(R^2 - x^2) - x*x1)/(2*R))^2,  # eq3
       1 + y3*sqrt(R^2 - x^2)/(x3*(R + x)),  # eq4
       1 + y2*sqrt(R^2 - x^2)/(x2*(-R + x)),  # eq5
       -2*R*sqrt(R^2 - x^2) + r1*(2*R + sqrt(R^2 - x^2 + (R - x)^2) + sqrt(R^2 - x^2 + (R + x)^2)),  # eq6
       -r2^2 + (x2 - (R^2 + R*x + R*x2 - x*x2 - y2*sqrt(R^2 - x^2))/(2*R))^2 + (y2 - (R*y2 + R*sqrt(R^2 - x^2) + x*y2 - x2*sqrt(R^2 - x^2))/(2*R))^2,  # eq7
       -r3^2 + (x3 - (-R^2 + R*x + R*x3 + x*x3 + y3*sqrt(R^2 - x^2))/(2*R))^2 + (y3 - (R*y3 + R*sqrt(R^2 - x^2) - x*y3 + x3*sqrt(R^2 - x^2))/(2*R))^2,  # eq8
   ]
end;
(r2, r3) = (10, 5)

iniv = BigFloat[50, -10, 20, -7, 23, 33, -35, 27]
res = nls(H, ini=iniv)

   (BigFloat[50.00000000000000000000000000000000000000000000000000000000000000000000000630095, -14.00000000000000000000000000000000000000000000000000000000000000000000000571507, 20.00000000000000000000000000000000000000000000000000000000000000000000000238607, -10.00000000000000000000000000000000000000000000000000000000000000000000000404504, 24.00000000000000000000000000000000000000000000000000000000000000000000000337377, 32.00000000000000000000000000000000000000000000000000000000000000000000000631642, -36.00000000000000000000000000000000000000000000000000000000000000000000000619483, 27.00000000000000000000000000000000000000000000000000000000000000000000000251126], true)

乙円，丙円の直径が 10, 5 のとき，甲円の直径は 40 である。

その他のパラメータは以下の通り。
R = 50;  x = -14;  y = 48;  r1 = 20;  x1 = -10;  r2 = 10;  x2 = 24;  y2 = 32;  r3 = 5;  x3 = -36;  y3 = 27

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (10, 5)
   (R, x, r1, x1, x2, y2, x3, y3) = res[1]
   y = sqrt(R^2 - x^2)
   @printf("甲円の直径 = %g\n", 2r1)
   @printf("R = %g;  x = %g;  y = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n",
       R, x, y, r1, x1, r2, x2, y2, r3, x3, y3)
   plot([R, x, -R, R], [0, y, 0, 0], color=:brown, lw=0.5)
   circle(0, 0, R, :black, beginangle=0, endangle=180)
   circle(x1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x3, y3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, r1, "甲円:r1 \n(x1,r1) ", :red, :right, :vcenter)
       point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :center, :top, delta=-delta)
       point(x3, y3, " 丙円:r3,(x3,y3)", :green, :left, :bottom, delta=delta)
       point(0, R, " R", :black, :left, :bottom, delta=delta)
   end
end;