算額（その802）

算額（その802）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

菱形の中に大円，中円，小円が入っている。小円と中円の直径の積が 1 寸（平方寸），双股（小円と中円が菱形と接する 2 点間の距離）が 5 寸のとき，大円の直径はいかほどか。

菱形の対角線を 2a, 2b
大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (0, r1 + r2)
小円の半径と中心座標を r3, (r1 + r3, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive, r3::positive
c = sqrt(a^2 + b^2)
sinθ = b/c
cosθ = a/c
eq1 = 2r2*2r3 - 1
eq2 = c - (c - r1 - r2)*sinθ - (c -r1 - r3)*cosθ - 5//2
eq3 = r1 - a*sinθ
eq4 = r2 - (b - r1 - r2)*cosθ
eq5 = r3 - (a - r1 - r3)*sinθ;

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 Float64.(v), r.f_converged
end;

function H(u)
   (a, b, r1, r2, r3) = u
   return [
       4*r2*r3 - 1,  # eq1
       -a*(-r1 - r3 + sqrt(a^2 + b^2))/sqrt(a^2 + b^2) - b*(-r1 - r2 + sqrt(a^2 + b^2))/sqrt(a^2 + b^2) + sqrt(a^2 + b^2) - 2.5,  # eq2
       -a*b/sqrt(a^2 + b^2) + r1,  # eq3
       -a*(b - r1 - r2)/sqrt(a^2 + b^2) + r2,  # eq4
       -b*(a - r1 - r3)/sqrt(a^2 + b^2) + r3,  # eq5
   ]
end;

iniv = BigFloat[4, 5, 0.3, 3, 0.74]
res = nls(H, ini=iniv)

   ([5.0, 3.75, 3.0, 0.3333333333333333, 0.75], true)

大円の直径は 6 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r1, r2, r3) = res[1]
   @printf("大円の直径 = %g\n", 2r1)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:magenta, lw=0.5)
   circle(0, 0, r1)
   circle(0, r1 + r2, r2, :blue)
   circle(r1 + r3, 0, r3, :green)
   if more == true
       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;

算額（その801）

算額（その801）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

大円の周りに中円と小円が，大円と同時に隣の円とも外接するように配置されている。
大円と小円の直径がそれぞれ 21707 寸，5759 寸のとき，中円の直径はいかほどか。
なお，大円と小円の直径が上述の数値の場合には精要算法に示された図とはイメージがかなり異なる。小円の直径が 6800 寸程度の場合によく似た図（下図）になる。

大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (0, r1 + r2)
小円の半径と中心座標を r3, (x3, y3), (x4, y4), (x5, y5)
中円と小円の中心から大円の中心へ引いた2本の直線のなす角を θ
とおき，以下の連立方程式を解く。
eq3，eq4, eq5 は第二余弦定理を援用した。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1, r2, r3, x3, y3, x4, y4, x5, y5, θ
eq1 = x3^2 + y3^2 - (r1 + r3)^2
eq2 = x3^2 + (r1 + r2 - y3)^2 - (r2 + r3)^2
eq3 = (x4 - x3)^2 + (y3 - y4)^2 - 4r3^2
eq4 = (r1 + r2)^2 + (r1 + r3)^2 - 2(r1 + r2)*(r1 + r3)cos(θ) - (x3^2 + (r1 + r2 - y3)^2)
eq5 = 2(r1 + r3)^2 - 2(r1 + r3)^2*cos(PI/3 - θ) - ((x3 - x4)^2 + (y3 - y4)^2)
eq6 = (r1 + r2)^2 + (r1 + r3)^2 - 2(r1 + r2)*(r1 + r3)cos(PI/3) - (x4^2 + (r1 + r2 - y4)^2);

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r2, x3, y3, x4, y4, θ))

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 Float64.(v), r.f_converged
end;

function H(u)
   (r2, x3, y3, x4, y4, θ) = u
   return [
       x3^2 + y3^2 - (r1 + r3)^2,  # eq1
       x3^2 - (r2 + r3)^2 + (r1 + r2 - y3)^2,  # eq2
       -4*r3^2 + (-x3 + x4)^2 + (y3 - y4)^2,  # eq3
       -x3^2 + (r1 + r2)^2 + (r1 + r3)^2 - (r1 + r3)*(2*r1 + 2*r2)*cos(θ) - (r1 + r2 - y3)^2,  # eq4
       -2*(r1 + r3)^2*sin(θ + pi/6) + 2*(r1 + r3)^2 - (x3 - x4)^2 - (y3 - y4)^2,  # eq5
       -x4^2 - (r1/2 + r3/2)*(2*r1 + 2*r2) + (r1 + r2)^2 + (r1 + r3)^2 - (r1 + r2 - y4)^2,  # eq6
   ]
end;

(r1, r3) = (21707, 5759) .// 2
# (r1, r3) = (21707, 6800) .// 2
r1r3 = r1 + r3
iniv = BigFloat[r1r3, r1r3*cos(pi/3), r1r3*sin(pi/3), r1r3*cos(pi/6), r1r3*sin(pi/6), pi/6]
res = nls(H, ini=iniv)

   ([8893.500000000004, 8031.882766640407, 11139.306451612903, 11893.126870171644, 6866.499999999951, 0.624707483913989], true)

中円の直径は 17787 寸である。

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

r2 = 8893.5;  x3 = 8031.88;  y3 = 11139.3;  x4 = 11893.1;  y4 = 6866.5;  θ = 0.624707

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x3, y3, x4, y4, θ) = res[1]
   @printf("中円の直径 = %g\n", 2r2)
   @printf("r2 = %g;  x3 = %g;  y3 = %g;  x4 = %g;  y4 = %g;  θ = %g\n", r2, x3, y3, x4, y4, θ)
   plot()
   circle(0, 0, r1)
   rotate(0, r1 + r2, r2, :blue)
   #circle2(x3, y3, r3, :green)
   #circle2(x4, y4, r3, :green)
   rotate(x3, y3, r3, :green)
   rotate(x4, y4, r3, :green)
   x5 = (r1 + r3)*cos(θ - pi/6)
   y5 = (r1 + r3)*sin(θ - pi/6)
   rotate(x5, y5, r3, :green)
   if more == true
       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;

算額（その800）

算額（その800）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

底面，上面が菱形の角錐台の中に直径 6 寸の球がカッチリと入っている。上面の菱形の対角線の長いほうが 4 寸，短いほうが 3 寸のとき，底面の菱形の長いほうの直径はいかほどか。

上面の菱形の対角線 H，I (H > I)
底面の菱形の対角線 A，B (A > B) 
球の半径を r とおく。
z 軸の上方から下を見ると下図のように見える。

赤の菱形が上面，青の菱形が底面である。赤い円が球を表す。これでは球は角錐台に接しているとは思えないが，次の図を見ると納得できる。
球は平面 ABIH に接している。接点は図に示す直線 Oa 上にある。Oa は AB, HI と直交する。
Oa を含み x-y 平面に垂直な面で切断して横から見たのが下図である。

直線 ah が上図における Oa である。

r = 3, b = 1.2 のとき，以下の方程式を解き，a = 15/2, h = 50/7 である。

底面と上面の相似比は a/b = 25/4 である。
したがって，上面の菱形の長い方の対角線が 4 寸のとき，底面上面の菱形の長い方の対角線は 25 寸である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r, a, b, h
r = 3
b = 12//10
eq1 = (2r)^2 + (a - b)^2 - (a + b)^2
eq2 = (h - 2r)/b - h/a
(a, h) = solve([eq1, eq2], (a, h))[1]

   (15/2, 50/7)

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 3
   b = 1.2
   a = 15/2
   h = 50/7
   @printf("相似比 = %g/%g = %g\n", a, b, a/b)
   plot([0, a, 0, 0], [0, 0, h, 0], color=:black, lw=0.5)
   segment(0, 2r, b, 2r)
   circle(0, r, r)
   if more == true
       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, r, " r", :red, :left, :vcenter)
       point(0, h, " h", :black, :left, :bottom, delta=delta)
       point(b, 2r, "(b,2r)", :black, :left, :bottom, delta=delta)
       point(0, 2r, " 2r", :black, :left, :bottom, delta=delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
   end
   savefig("/Users/aoki/Downloads/fig0.png");
   A = 25/2
   B = A*3/4
   H = 4/2
   I = H*3/4
   x0 = B*4/5*3/5
   y0 = B*4/5*4/5
   plot([A, 0, -A, 0, A], [0, B, 0, -B, 0], color=:blue, lw=0.5)
   plot!([H, 0, -H, 0, H], [0, I, 0, -I, 0], color=:red, lw=0.5)
   circle(0, 0, r, :red)
   circle(0, 0, 7.5, :gray80)
   segment(0, 0, x0, y0)
   hline!([0], color=:gray80, lw=0.5)
   vline!([0], color=:gray80, lw=0.5)
   point(A, 0, " A", :blue, :left, :vcenter)
   point(0, B, "B", :blue, :center, :bottom, delta=2delta)
   point(H, 0, " H", :red, :left, :vcenter)
   point(0, I, "I", :red, :center, :bottom, delta=2delta)
   point(x0, y0, " a", :blue, :left, :bottom, delta=2delta)
   point(0, 0, "O", :black, :right, delta=-delta)
end;

算額（その799）

算額（その799）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

等脚台形の中に斜線を引き，区画された領域に甲円と乙円を入れる。乙円の直径が 15 寸，上頭（上底）の長さが 12 寸，斜線の長さが 20 寸のとき，甲円の直径はいかほどか。

等脚台形の高さを h，下頭（下底），上頭の長さをそれぞれ 2a, 2b
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (0, h - r2)
とおいて以下の連立方程式を解く。

なお，斜の長さについては算法助術の公式 39 を援用した。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a::positive, b::positive, h::positive,
     r1::negative, r2::positive, 斜::positive, d
eq1 = h^2 + (a - b)^2 - (斜 + a + b)^2
eq2 = dist(a, 0, b, h, 0, r1) - r1^2
eq2 = numerator(apart(eq2, d))
eq3 = dist(a, 0, b, h, 0, h - r2) - r2^2
eq3 = numerator(apart(eq3, d))
res = solve([eq1, eq2, eq3], (r1, a, h));
res[2]

   (-b*斜/(2*r2) + r2 + r2*斜/(2*b), -斜/2 + r2^2/b + r2^2*斜/(2*b^2), r2*(2*b + 斜)/b)

2 組の解が得られるが，2 番目のものが適解である。

b = 12//2
斜 = Sym(20)
r2 = 15//2
(-b*斜/(2*r2) + r2 + r2*斜/(2*b), -斜/2 + r2^2/b + r2^2*斜/(2*b^2), r2*(2*b + 斜)/b)

   (12, 15, 40)

甲円の半径の一般解は r2*(1 + 斜/2b) - b*斜/2r2 である

r2*(1 + 斜/2b) - b*斜/2r2 |> println

   12

図を描くために斜と等脚辺の交点座標 (x1, y1), (x2, y2) を求める。
一般解は求めにくいようなので，問の条件下で交点座標を求める。

@syms x1::positive, y1::positive, x2::negative, y2::positive, d
(a, h, r1) = (15, 40, 12)
eq4 = dist(x1, y1, x2, y2, 0, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = dist(x1, y1, x2, y2, 0, h - r2) - r2^2
eq5 = numerator(apart(eq5, d))
eq6 = y1*(b - a) - h*(x1 - a)
eq7 = y2*(a - b) - h*(x2 + a)
res2 = solve([eq4, eq5, eq6, eq7], (x1, y1, x2, y2));
res2[2]

   (9*sqrt(10)/41 + 390/41, 1000/41 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41, 40*sqrt(10)/41 + 1000/41)

2 組の解が得られるが，2 番目のものが適解である。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 12//2
   斜 = Sym(20)
   r2 = 15//2
   (a, h, r1) = (15, 40, 12)
   @printf("甲円の直径 = %g;  下頭 = %g;  高さ = %g\n", 2r1, 2a, h)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, h - r2, r2, :blue)
   if more == true
       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)
       (x1, y1, x2, y2) = (9*sqrt(10)/41 + 390/41, 1000/41 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41, 40*sqrt(10)/41 + 1000/41)
       segment(x1, y1, x2, y2, :green)
       point(b, h, "(b,h)", :black, :right, :bottom, delta=delta/2)
       point(a, 0, "(a,0) ", :black, :right, :bottom, delta=delta/2)
       point(x1, y1, "(x1,y1)", :black, :right, delta=-delta)
       point(x2, y2, "(x2,y2) ", :black, :left, delta=-2delta)
       point(0, r1, " 甲円:r1,(0,r1)", :red, :left, :vcenter)
       point(0, h - r2, " 乙円:r2\n (0,h-r2)", :blue, :left, :vcenter)
   end
end;

算額（その798）

算額（その798）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

半円の中に離矢を隔てて大円と小円が入っている。両円の直径の和が 4 寸，離矢が 2 寸 4 分のとき，外円（半円）の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (x2, r2)
2 円の直径の和を α，離矢を β
とおいて以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r1::positive, x1::positive, r2::positive,
     x2::negative, R::positive,
     α::positive, β::positive
eq1 = 2r1 + 2r2 - α
eq2 = x1^2 + r1^2 - (R - r1)^2
eq3 = x2^2 + r2^2 - (R - r2)^2
eq4 = (R^2 - (x1 - r1)^2) - β^2
eq5 = (x2 + r2) - (x1 - r1)
res = solve([eq1, eq2, eq3, eq4, eq5], (R, r1, x1, r2, x2));

2 通りの解が得られるが，大円と小円の位置関係だけが違うものである。
問の位置関係からいえば，2 番目のものが適解である。
この場合，外円の直径は α^2/(4(2β - α)) となる。
2 円の直径の和 α が 4 寸，離矢 β が 2.4 寸ならば，外円の直径は 5 寸である。

α = 4
β = 2.4
α^2/(4(2β - α))

   5.000000000000001

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

R = 2.5;  r1 = 1.2;  x1 = 0.5;  r2 = 0.8;  x2 = -1.5

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (α, β) = (4, 2.4)
   (R, r1, x1, r2, x2) = (5//2, 1.2, 0.5, 0.8, -1.5000000000000002)
   @printf("2円の直径の和が %g，離矢が %g のとき，外円の直径は %g 寸である。\n", α, β, 2R)
   @printf("R = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g\n", R, r1, x1, r2, x2)
   plot()
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, r1, r1, :blue)
   circle(x2, r2, r2, :green)
   segment(x1 - r1, 0, x1 - r1, β, :magenta)
   if more == true
       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(x1, r1, "大円:r1,(x1,r1)", :blue, :center, delta=-delta)
       point(x2, r2, "小円:r2,(x2,r2)", :green, :center, delta=-delta)
       point(x1 - r1, β, "(x1-r1,β) ", :magenta, :right, :bottom, delta=delta)
       point(x1 - r1, 3β/4, "離矢 ", :magenta, :right, :vcenter, mark=false)
       point(R, 0, " R", :red, :left, :bottom, delta=delta)
   end
end;

算額（その797）

算額（その797）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

正方形の中に2本の斜線を引き，できた区画の中に大円，中円，小円を入れる。大円，中円，小円の直径および2本の斜線の長さと正方形の一辺の長さの和が 55 寸のとき，正方形の一辺の長さはいかほどか。

大円の半径と中心座標を r1, (r1, a - r1)
中円の半径と中心座標を r2, (a - r2, y2)
小円の半径と中心座標を r3, (r3, r3)
とおいて以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r1::positive, r2::positive, y2::positive,
     r3::positive, a::positive, b::positive,
     c::positive, d::positive
eq1 = dist(b, 0, a, a, r1, a - r1) - r1^2
eq2 = dist(b, 0, a, a, a - r2, y2) - r2^2
eq3 = dist(b, 0, a, a, r3, r3) - r3^2
eq4 = dist(c,0, 0, d, r1, a - r1) - r1^2
eq5 = dist(c,0, 0, d, a - r2, y2) - r2^2
eq7 = c + d - sqrt(c^2 + d^2) - 2r3
eq8 = a + (a - b) - sqrt(a^2 + (a - b)^2) - 2r2
eq9 = 2r1 + 2r2 + 2r3 + sqrt(a^2 + (a - b)^2) + sqrt(c^2 + d^2) + a - 55;

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 Float64.(v), r.f_converged
end;

function H(u)
   (r1, r2, y2, r3, a, b, c, d) = u
   return [
       -r1^2 + (a - a*(a*(a - r1) + (a - b)*(-b + r1))/(a^2 + (a - b)^2) - r1)^2 + (-b + r1 - (a - b)*(a*(a - r1) + (a - b)*(-b + r1))/(a^2 + (a - b)^2))^2,  # eq1
       -r2^2 + (-a*(a*y2 + (a - b)*(a - b - r2))/(a^2 + (a - b)^2) + y2)^2 + (a - b - r2 - (a - b)*(a*y2 + (a - b)*(a - b - r2))/(a^2 + (a - b)^2))^2,  # eq2
       -r3^2 + (-a*(a*r3 + (a - b)*(-b + r3))/(a^2 + (a - b)^2) + r3)^2 + (-b + r3 - (a - b)*(a*r3 + (a - b)*(-b + r3))/(a^2 + (a - b)^2))^2,  # eq3
       -r1^2 + (a - d*(-c*(-c + r1) + d*(a - r1))/(c^2 + d^2) - r1)^2 + (-c + c*(-c*(-c + r1) + d*(a - r1))/(c^2 + d^2) + r1)^2,  # eq4
       -r2^2 + (-d*(-c*(a - c - r2) + d*y2)/(c^2 + d^2) + y2)^2 + (a - c + c*(-c*(a - c - r2) + d*y2)/(c^2 + d^2) - r2)^2,  # eq5
       c + d - 2*r3 - sqrt(c^2 + d^2),  # eq7
       2*a - b - 2*r2 - sqrt(a^2 + (a - b)^2),  # eq8
       a + 2*r1 + 2*r2 + 2*r3 + sqrt(a^2 + (a - b)^2) + sqrt(c^2 + d^2) - 55,  # eq9
   ]
end;

iniv = BigFloat[44, 30, 30, 21, 124, 35, 85, 62]./10
res = nls(H, ini=iniv)

   ([4.0, 3.0, 3.0, 2.0, 12.0, 3.0, 8.0, 6.0], true)

正方形の一辺の長さは 12 寸である。

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

r1 = 4;  r2 = 3;  y2 = 3;  r3 = 2;  a = 12;  b = 3;  c = 8;  d = 6

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y2, r3, a, b, c, d) = res[1]
   @printf("正方形の一辺の長さは %g 寸である。\n", a)
   @printf("r1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  a = %g;  b = %g;  c = %g;  d = %g\n", r1, r2, y2, r3, a, b, c, d)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(r1, a - r1, r1)
   circle(a - r2, y2, r2, :blue)
   circle(r3, r3, r3, :green)
   segment(c, 0, 0, d, :magenta)
   segment(b, 0, a, a, :magenta)
   if more == true
       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(r1, a - r1, "大円:r1,(r1,a-r1)", :red, :center, delta=-delta)
       point(a - r2, y2, "中円:r2,(a-r2,y2)", :blue, :center, delta=-delta)
       point(r3, r3, "小円:r3,(r3,r3)", :green, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(b, 0, "   b", :magenta, :left, :bottom, delta=delta/2)
       point(c, 0, "c   ", :magenta, :right, :bottom, delta=delta/2)
       point(0, d, "d ", :magenta, :right, :vcenter)
       plot!(xlims=(-0.5, 12.5))
   end
end;

算額（その796）

算額（その796）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直角三角形内に界斜（斜線）によって隔てられた領域に全円と小円が入っている。全円は三篇と界斜に接し，小円は二辺と界斜に接している。全円，小円の直径がそれぞれ 8 寸，2 寸のとき，界斜の長さはいかほどか。

全円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (x2, r1)
斜線と直角三角形の辺との交点座標を (c, 0), (d, e)
とおいて以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a::positive, b::positive, c::positive,
     d::positive, e::positive, r1::positive,
     r2::positive, x2::positive
#x2 = b - b*r2/r1 + r2
eq1 = (a + b - 2r1) - sqrt(a^2 + b^2)
eq4 = dist(0, a, b, 0, r1, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = r2/(b - x2) - r1/(b - r1)
#eq5 = r2*(b - r1) - r1*(b - x2);
#eq5 = dist(0, a, b, 0, x2, r2) - r2^2
#eq5 = (sqrt(a^2 + b^2) + b + a)*r1 - a*b

eq2 = dist(c, 0, d, e, r1, r1) - r1^2
eq2 = numerator(apart(eq2, d))/e
eq3 = dist(c, 0, d, e, x2, r2) - r2^2
eq3 = numerator(apart(eq3, d))/e
eq6 = e*b - a*(b - d);

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 Float64.(v), r.f_converged
end;

function H(u)
   (a, b, c, d, e, x2) = u
   return [
       a + b - 2*r1 - sqrt(a^2 + b^2),  # eq1
       c^2*e - 2*c^2*r1 + 2*c*d*r1 - 2*c*e*r1 + 2*c*r1^2 - 2*d*r1^2,  # eq2
       c^2*e - 2*c^2*r2 + 2*c*d*r2 - 2*c*e*x2 + 2*c*r2*x2 - 2*d*r2*x2 - e*r2^2 + e*x2^2,  # eq3
       a^2*b^2 - 2*a^2*b*r1 - 2*a*b^2*r1 + 2*a*b*r1^2,  # eq4
       -r1/(b - r1) + r2/(b - x2),  # eq5
       -a*(b - d) + b*e,  # eq6
   ]
end;

(r1, r2) = (8, 2) .// 2
iniv = BigFloat[12, 16, 4.6, 13.2, 2.1, 13]
res = nls(H, ini=iniv)

   ([12.0, 16.0, 4.468871125850725, 13.22490309931942, 2.0813226755104353, 13.0], true)

c = 4.468871125850725, d = 13.22490309931942, e = 2.0813226755104353 となり，斜線の長さは 9 寸である。

c = 4.468871125850725
d = 13.22490309931942
e = 2.0813226755104353
sqrt((d - c)^2 + e^2)

   9.0

その他のパラメータは以下のとおりである。
a = 12;  b = 16;  c = 4.46887;  d = 13.2249;  e = 2.08132;  x2 = 13
なお，x2 = b - b*r2/r1 + r2 の関係がある。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (8, 2) .// 2
   (a, b, c, d, e, x2) = res[1]
   # x2 = b - b*r2/r1 + r2
   l = sqrt((d - c)^2 + e^2)
   @printf("全円の直径が %g，小円の直径が %g のとき，界斜の長さは %g である。\n", 2r1, 2r2, l)
   @printf("a = %g;  b = %g;  c = %g;  d = %g;  e = %g;  x2 = %g\n", a, b, c, d, e, x2)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :blue)
   segment(c, 0, d, e, :green)
   if more == true
       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(r1, r1, "全円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "全円:r2,(x2,r2) ", :black, :right, delta=-delta)
       point(0, a, " a", :black, :left, :bottom)
       point(b, 0, "b", :black, :left, :bottom, delta=delta)
       point(c, 0, "(c,0)", :green, :center, :bottom, delta=2delta)
       point(d, e, "(d,e)", :green, :center, :bottom, delta=2delta)
   end
end;

算額（その795）

算額（その795）

宮崎智也，吉岡倖佑，原本博史，安部利之：大西佐兵衛『雑題』の現代解について ~第8巻 第1問から第4問~，愛媛大学教育学部附属 科学教育研究センター紀要，Vol.2，pp.17-22，2023/03/31.
https://ed.ehime-u.ac.jp/CRESE/wp-content/uploads/2023/03/Vol2No3.pdf

円弧内に等円と小円を入れる。等円の直径と元の長さを用いて，分数の形で小円の直径を求めよ。

弦の長さを 2x0，弦と y 軸の交点座標を (0, y0)
円弧がその一部である円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r1, (r1, y0 + r1)
小円の半径と中心座標を r2, (x2, y0 + r2)
とおき，以下の連立方程式を解く。

たとえば，等円の直径が 18 寸，弦の長さが 54 寸のとき，小円の直径は 8 寸である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, r2::positive,
     x2::positive, x0::positive, y0::positive
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = r1^2 + (y0 + r1)^2 - (R - r1)^2
eq3 = x2^2 + (y0 + r2)^2 - (R - r2)^2
eq4 = x0^2 + y0^2 - R^2;
res = solve([eq1, eq2, eq3, eq4], (R, r2, x2, y0));

小円の半径は，r1, x0 を使って以下のように表すことができる。

r2 = res[1][2]
r2|> println

   r1*(-r1 + x0)^2*(r1 + x0)^2/(3*r1^2 + x0^2)^2

等円の直径が 18 寸，弦の長さが 54 寸のとき，小円の直径は 8 寸である。

r2(r1 => 18//2, x0 => 54//2) |> println

   4

その他のパラメータは以下のとおりである。
r1 = 9;  x0 = 27;  R = 28.125;  r2 = 4;  x2 = 21;  y0 = 7.875

SymPy の簡約化は中途半端なので，手動で簡約化する。

r22 = r1*(r1^2 - x0^2)^2 / (3r1^2 + x0^2)^2;
r2(r1 => 18//2, x0 => 54//2) |> println
r22(r1 => 18//2, x0 => 54//2) |> println

   4
   4

直径 d1 = 2r1, 弦の長さ x = 2x0 を使って等円の「直径」を求めるときは，次のようになる。

@syms d1, x
r23 = d1*(d1^2 - x^2)^2/(3*d1^2 + x^2)^2
r23(d1 => 18, x => 54) |> println

   8

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x0) = (18, 54) .// 2
   (R, r2, x2, y0) = (225/8, 4, 21, 63/8)
   @printf("等円の直径が %g，弦の長さが %g のとき，小円の直径は %g\n", 2r1, 2x0, 2r2)
   @printf("r1 = %g;  x0 = %g;  R = %g;  r2 = %g;  x2 = %g;  y0 = %g\n", r1, x0, R, r2, x2, y0)
   θ = atand(y0, x0)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180-θ, n=500)
   circle2(r1, y0 + r1, r1, :green)
   circle2(x2, y0 + r2, r2, :blue)
   segment(-x0, y0, x0, y0)
   segment(0, 0, x0, y0)
   if more == true
       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, R, " R", :black, :left, :bottom, delta=delta)
       point(0, y0, " y0", :black, :left, delta=-delta)
       point(x0, y0, "(x0,y0)", :black, :right, delta=-4delta)
       point(r1, y0 + r1, "等円:r1,(r1,y0+r1)", :green, :center, delta=-delta)
       point(x2, y0 + r2, "小円:r2,(x2,y0+r2)", :black, :right, delta=-delta)
   end
end;

算額（その794）

算額（その794）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

福島県郡山市日和田町宮下 天満宮 明治14年(1881)
http://www.wasan.jp/fukusima/miyasitatenma.html

後者の写真が不鮮明で長い間お蔵入りになっていたが，「精要算法」にあった。問題文と答えも同じである。ただし，後者の図では乙円と戊円は接していないが，精要算法にあるように，乙円と戊円も接している。

外円内に 7 個の円が入っている。外円，甲円の直径がそれぞれ 890 寸，267 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (0, r4 - R)
戊円の半径と中心座標を r5, (0, r5 + 2r4 - R)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, x2::positive, y2::positive, x3::positive, y3::negative, r1::positive, r2::positive, r3::positive, r4::positive, r5::positive;
# (R, r1) = (890,267) .// 2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq3 = x3^2 + (r5 + 2r4 - R - y3)^2 - (r3 + r5)^2
eq4 = x3^2 + (y3 - r4 + R)^2 - (r3 + r4)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = x3^2 + y3^2 - (R - r3)^2
eq7 = x2^2 + (r5 + 2r4 - R - y2)^2 - (r2 + r5)^2
eq8 = r1 + r4 + r5 - R;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8],
    (x2, y2, x3, y3, r2, r3, r4, r5));
res[2] |> println

   (4*sqrt(3)*R*r1*(R - r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -R*(-3*R^2 + 6*R*r1 + r1^2)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*sqrt(3)*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R^2 + 2*R*r1 + 11*r1^2)/(R^2 + 2*R*r1 + 9*r1^2), -4*R*r1*(-R + r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R + r1)/(R + 3*r1), (3*R*r1 - 3*r1^2)/(R + 3*r1))

2 組の解が得られるが，2 番目のものが適解である。

外円，甲円の直径がそれぞれ 890 寸，267 寸のとき，乙円の直径は 8*R*r1*(R - r1)/(3*R^2 - 2*R*r1 + 3*r1^2) = 280 寸である。

2res[2][5] |> println
2res[2][5](R => 890//2, r1 => 267//2) |> println

   -8*R*r1*(-R + r1)/(3*R^2 - 2*R*r1 + 3*r1^2)
   280

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (890,267) ./ 2
   (x2, y2, x3, y3, r2, r3, r4, r5) = (4*sqrt(3)*R*r1*(R - r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -R*(-3*R^2 + 6*R*r1 + r1^2)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*sqrt(3)*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R^2 + 2*R*r1 + 11*r1^2)/(R^2 + 2*R*r1 + 9*r1^2), -4*R*r1*(-R + r1)/(3*R^2 - 2*R*r1 + 3*r1^2), -4*R*r1*(-R + r1)/(R^2 + 2*R*r1 + 9*r1^2), -R*(-R + r1)/(R + 3*r1), (3*R*r1 - 3*r1^2)/(R + 3*r1))
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r1, r1)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
   circle(0, r4 - R, r4, :brown)
   circle(0, r5 + 2r4 - R, r5, :blue)
   if more == true
       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, R, "外円:R,(0,0)", :red, :center, :bottom, delta=delta/2)
       point(0, R-r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :green, :center, delta=-delta)
       point(x3, y3, "丙円:r3,(x3,y3)", :magenta, :center, delta=-delta)
       point(0, r4 - R, "丁円:r4,(0,r4-R)", :brown, :center, delta=-delta)
       point(0, r5 + 2r4 - R, "戊円:r5,(0,r5+2r4-R)", :blue, :center, delta=-delta)
   end
end;

 

算額（その793）

算額（その793）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直線上に互いに接する 3 個の円（甲円，丙円，丁円）が載っており，その 3 個の円のいずれとも接するように乙円が載っている。甲円，丙円，丁円の直径がわかっているときに，乙円の直径を求めよ。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
なお，SymPy の能力では一般解は求まらないので，与えられた条件における解を求める。なお，後半で r2 の一般解を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     x2::positive, y2::positive, r3::positive,
     r4::positive, x4::positive
(r1, r3, r4) = (100, 64, 48) .// 2
eq1 = (x1 - x2)^2 + (y2 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2 |> expand
eq3 = x2^2 + (y2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = (x4 - x2)^2 + (y2 - r4)^2 - (r2 + r4)^2 |> expand
eq5 = x4^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand;
res = solve([eq1, eq2, eq3, eq4, eq5], (x1, r2, x2, y2, x4))
(x1, r2, x2, y2, x4) = res[1]
@printf("乙円の直径 = %g\n", 2res[1][2])
println("(x1, r2, x2, y2, x4) = ",res[1])
@printf("x1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x4 = %g\n", x1, r2, x2, y2, x4)

   乙円の直径 = 72.9
   (x1, r2, x2, y2, x4) = (72*sqrt(3), 729/20, 26*sqrt(3), 1671/20, 32*sqrt(3))
   x1 = 124.708;  r2 = 36.45;  x2 = 45.0333;  y2 = 83.55;  x4 = 55.4256

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3, r4) = (100, 64, 48) .// 2
   (x1, r2, x2, y2, x4) = (72*sqrt(3), 729/20, 26*sqrt(3), 1671/20, 32*sqrt(3))
   plot()
   circle(x1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(0, r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta)
       point(0, r3, "丙円:r3,(0,r3)", :green, :center, delta=-delta)
       point(x4, r4, "丁円:r4\n(x4,r4)", :magenta, :center, delta=-delta)
   end
end;

1.　算法助術の公式47

直線上に載って互いに接している円の半径をそれぞれ r1，r2 とする。その 2 つの円に載っている円の半径を r3 とする。
上に載っている円のてっぺんから直線までの距離 h を求める。

公式集は，直径についての式になっているので，これを半径についての式に書き直す。
それぞれの円の直径を d1 = 2r1, d2 = 2r2, d3 = 2r3 とする。

d1*h + d2*h - d1*d2 = 2sqrt(d1*d2*d3*h)
r1*h + r2*h - 2r1*r2 = 2sqrt(2r1*r2*r3*h)

using SymPy
@syms r1, r2, r3, r4, h, d1, d2, d3, d4
eq47d = d1*h + d2*h - d1*d2 - 2sqrt(d1*d2*d3*h)
eq47r = 2r1*h + 2r2*h - 2r1*2r2 - 2sqrt(2r1*2r2*2r3*h)
eq47r = r1*h + r2*h - 2r1*r2 - 2sqrt(2r1*r2*r3*h)

2.　応用

乙円のてっぺんから直線までの距離を h とする。
甲円，乙円，丙円，丁円の半径をそれぞれ r1, r2, r3, r4 とする。
まず，甲円，丁円，乙円に対して公式47を適用する。
ただし，そのまま使って連立方程式を解くには SymPy の能力的に無理があるので，辺々を 2 乗して使う。

using SymPy
@syms r1, r2, r3, r4, h

eq11 = (r1*h + r4*h - 2r1*r4)^2 - 4(2r1*r4*r2*h);

同様に，丙円，丁円，乙円に対して公式47を適用する。

eq12 = (r3*h + r4*h - 2r3*r4)^2 - 4(2r3*r4*r2*h);

筆算で解く場合にはまず h について解き，それをどちらかの式に代入して r2 を求めるのだろうが，SymPy で連立方程式を解くと，自動的にやってくれる。

res = solve([eq11, eq12], (h, r2));
res[2]

   (4*r1*r3*(r4^2*(r1 + r3 + 2*r4)/(2*(r1*r3 - r4^2)) + r4 + r4^2*sqrt(r1*r3)*(2*r1*r3 + r1*r4 + r3*r4)/(2*r1*r3*(r1*r3 - r4^2)))/(2*r1*r3 + r1*r4 + r3*r4), r4^2*(r1 + r3 + 2*r4)/(4*(r1*r3 - r4^2)) + r4^2*sqrt(r1*r3)*(2*r1*r3 + r1*r4 + r3*r4)/(4*r1*r3*(r1*r3 - r4^2)))

2 番目のものが適解である。

甲円の直径が 100 寸, 丙円 の直径が 64 寸, 丁円の直径が 48 寸のとき, h =120，r2 = 36.45 である。直径は 72.9 である。

res[2][1](r1 => 100/2, r3 => 64/2, r4 => 48/2) |> println
2res[2][2](r1 => 100/2, r3 => 64/2, r4 => 48/2) |> println

   120.000000000000
   72.9000000000000

r2 の式は，長いままで，SymPy では簡約化できない。

res[2][2] |> println

   r4^2*(r1 + r3 + 2*r4)/(4*(r1*r3 - r4^2)) + r4^2*sqrt(r1*r3)*(2*r1*r3 + r1*r4 + r3*r4)/(4*r1*r3*(r1*r3 - r4^2))

術文は，「天 = sqrt(甲*丙), 地 = 4(天 - 丁)として，乙 = ((甲 + 丙)/天 + 2)丁^2/地

@syms d1::positive, d3::positive, d4::positive
h = sqrt(d1*d3)d4/(sqrt(d1*d3) - d4);
d2 = ((d1 + d3)/sqrt(d1*d3) + 2)*d4^2/4(sqrt(d1*d3)- d4)
d2(d1 => 100, d3 => 64, d4 => 48).evalf() |> println

   72.9000000000000

 

算額（その792）

算額（その792）

文化元年甲子九月 藤田貞資門人 北越村上藩 杉村勘助道定
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

圓壔（えんとう；円柱）から弧環（トーラス？）で欠き取った鼓状の物体内に大球 2 個と，数個の小球が入っている。小球は大球と外側面および両隣の小球に接する（小球間の隙間はない）。円柱の底面（円）の直径が 35 寸 9 分のとき，円柱の高さはいかほどか。

円柱の高さ，底面の半径と中心座標を 4r2, r1, (0, 0, 2r2)
大球の半径と中心座標を r2, (0, 0, r2)
小球の半径と中心座標を r4, (x4, 0, 0)
とおき，以下の連立方程式を解くが，小球の個数を試行錯誤で決める必要がある（?）。
n = 10 のときに，条件を満たす解が得られた。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     x3::positive, r4::positive, x4::positive,
     h::positive, n::integer
h = 2r2
n = 10
eq1 = x3^2 + r2^2 - (r2 + r3)^2
eq2 = (x3 - r1)^2 + h^2 - r3^2
eq3 = x4^2 + r2^2 - (r2 + r4)^2
eq4 = x4 + r4 + r3 - x3
eq5 = x4*sind(Sym(180)//n) - r4
res = solve((eq1, eq2, eq3, eq4, eq5), (r2, r3, x3, r4, x4));

2 組の解が得られるが，2 番目のものが適解である。

res[2][1] |> println
res[2][4] |> println
res[2][5] |> println

   r1*(29*sqrt(5 + 5*sqrt(5)) + 65*sqrt(1 + sqrt(5)) + 125 + 56*sqrt(5))/(2*(40*sqrt(1 + sqrt(5)) + 18*sqrt(5)*sqrt(1 + sqrt(5)) + 105 + 47*sqrt(5)))
   -3*sqrt(5)*r1*sqrt(5 + 5*sqrt(5))/5 - sqrt(5)*r1*sqrt(1 + sqrt(5)) - r1/4 + 3*sqrt(5)*r1/20 + 9*r1*sqrt(1 + sqrt(5))/4 + 27*r1*sqrt(5 + 5*sqrt(5))/20
   -sqrt(5)*r1/10 + r1*sqrt(1 + sqrt(5))/(2*sqrt(5) + 5) + r1/2

円柱の底面（円）の直径が 35 寸 9 分のとき，円柱の高さは 49.0003 寸，小球の個数は 10 個で，小球の半径は 2.58655 寸，小球の中心を通る円の半径は 8.37025 寸である。

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

r2 = 12.2501;  r3 = 46.4137;  x3 = 57.3705;  r4 = 2.58655;  x4 = 8.37025

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 35.9/2
   u = √5
   t = sqrt(1 + u)
   s = u*t
   (r2, r3, x3, r4, x4) = (
       359(29s + 65t + 125 + 56u)/(40(18s + 40t + 47u + 105)),
       359(30254s + 67650t + 58301u + 130365)/(200(987s + 2207t + 2584u + 5778)),
       359(987s + 2207t + 4253 + 1902u)/(10(521s + 1165t + 1364u + 3050)),
       (2513s -5385t + 1077u - 1795)/400,
       -359u/200 + 359t/(20(2u + 5)) + 359/40)
   h = 2r2
   n = 10
   @printf("円柱の高さ = %g;  小球の個数 = %g;  小球の半径 = %g\n", 4r2, n, r4)
   @printf("r2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 = %g\n", r2, r3, x3, r4, x4)
   θ = atand(h, x3 - r1)
   plot()
   circle22(0, r2, r2, :blue)
   circle(x3, 0, r3, beginangle=180 - θ, endangle=180 + θ)
   circle(x4, 0, r4, :green)
   segment(-r1, h, r1, h, :black)
   segment(r1, h, x3, 0)
   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, r2, " 大球:r2\n (0, r2)", :blue, :left, :vcenter)
       point(x4, 0, "   小球:r4,(x4, 0)", :green, :left, :vcenter)
       point(r1, h, "底面:r1,(0,2r2)", :black, :center, :bottom, delta=delta)
       point(x3, 0, "弧環:r3,(x3,0)", :red, :right, delta=-delta/2)
   end

   function rotate2(ox, oy, r, color=:red; angle=120, beginangle=0, endangle=360, by=0.5, n=0)
       for deg in 0:angle:360-1
           (ox2, oy2) = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [ox; oy]
           @printf("(x-%g)^2+(y-%g)^2+z^2=%g^2\n", ox2, oy2, r)
           circle(ox2, oy2, r, color; beginangle, endangle, by, n)
           beginangle += angle
           endangle += angle
       end
   end;
    pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   circle(0, 0, r1, :black, lw=0.1)
   circle(0, 0, r2, :blue, lw=0.1)
   circle(0, 0, x4, :orange)
   rotate2(x4, 0, r4, :green, angle=2*18)
   point(x4, 0, "   小球:r4\n   (x4, 0)", :green, :left, :vcenter)
   point(0, r1, "圓壔", :black, :center, :vcenter, mark=false)
   point(0, r2, "大球", :blue, :center, :vcenter, mark=false)
end;

算額（その791）

算額（その791）

寛政十一年己未四月 丸山良玄門人 北越中宿邑 米持杢左衛門富房
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

右鈎，左鈎，容円の直径が 9 寸，5 寸，3 寸のとき，雙股はいかほどか。

右鈎，左鈎，雙股 を R, L, a
容円の半径と中心座標を r, (x, r)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms x::positive, a::positive, d,

     r::positive, L::positive, R::positive
eq1 = dist(0, L, a, 0, x, r) - r^2
eq2 = dist(a, R, 0, 0, x, r) - r^2;
eq1 = numerator(apart(eq1, d))
eq2 = numerator(apart(eq2, d))
eq1 |> println
eq2 |> println

   L^2*a^2 - 2*L^2*a*x - L^2*r^2 + L^2*x^2 - 2*L*a^2*r + 2*L*a*r*x
   -R^2*r^2 + R^2*x^2 - 2*R*a*r*x

2 組の解が得られるが，最初のものが適解である。

右鈎，左鈎，容円の直径が 9 寸，5 寸，3 寸のとき，雙股は 112 寸である。

res = solve([eq1, eq2], (a, x));
res[1][1] |> println
res[1][1](L => 5, R => 9, r => 3//2) |> println

   2*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2))*(-L*R + L*r + R*r)/(-L + 2*r)
   12

容円の中心座標は (9/2, 3/2) である。

res[1][2] |> println
res[1][2](L => 5, R => 9, r => 3//2) |> println

   R*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2))
   9/2

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (L, R, r) = (5, 9, 3/2)
   (a, x) = (2*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2))*(-L*R + L*r + R*r)/(-L + 2*r), R*r*sqrt((L - 2*r)/(L*R^2 - 4*L*R*r + 4*L*r^2 - 2*R^2*r + 4*R*r^2)))
   @printf("右鈎 = %g;  左鈎 = %g;  容円直径 = %g\n", R, L, 2r)
   @printf("雙股 = %g;  x = %g\n", a, x)
   plot([0, 0, a, 0, a, a, 0], [0, L, 0, 0, R, 0], color=:blue, lw=0.5)
   circle(x, r, r)
   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(a, R, "(a,R)", :blue, :right, :bottom, delta=delta/2)
       point(0, L, " R", :blue, :left, :bottom, delta=delta/2)
       point(x, r, "容円:r,(x,r)", :red, :center, delta=-delta/2)
       point(0, L/2, " 左鈎", :blue, :left, :vcenter, mark=false)
       point(a, R/2, "右鈎 ", :blue, :right, :vcenter, mark=false)
       point(a/2, 0, "雙股", :blue, :center, :bottom, delta=delta/2, mark=false)
   end
end;

算額（その790）

算額（その790）

寛政八年丙辰十一月 丸山因平良玄門人 参州苅屋 林政右衛門盛保 
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

十字線を隔てて乾円，坤円，巽円，艮円の4円と，中央に容円を置く。艮円，坤円，巽円の直径がそれぞれ 15 寸，10 寸，6 寸のとき，乾円の直径を求めよ。

容円の半径と中心座標を r0, (x0, y0)
坤円の半径と中心座標を r1, (r1, r1)
乾円の半径と中心座標を r2, (r2, -r2)
巽円の半径と中心座標を r3, (-r3, r3)
艮円の半径と中心座標を r4, (-r4, -r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r0::positive, x0::negative, y0::positive,
     r1::positive, r2::positive, r3::positive, r4::positive
eq1 = (r1 - x0)^2 + (r1 - y0)^2 - (r0 + r1)^2
eq2 = (r2 - x0)^2 + (-r2 - y0)^2 - (r0 + r2)^2
eq3 = (-r3 - x0)^2 + (r3 - y0)^2 - (r0 + r3)^2
eq4 = (-r4 - x0)^2 + (-r4 - y0)^2 - (r0 + r4)^2
res = solve([eq1, eq2, eq3, eq4], (r2, r0, x0, y0));

4 番目の組が適解であるが，SymPy では簡約化できない長い式になる。

res[4][1] |> println

   sqrt(8*r1^2*r3*r4^2*(r1*r3 - r1*r4 + r3^2 + r3*r4)*(r1^2*r3 + r1^2*r4 + r1*r3^2 + 6*r1*r3*r4 - r1*r3*sqrt(r1^2 + 6*r1*r4 + r4^2) + r1*r4^2 - r1*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + r3^2*r4 - r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + r3*r4^2 - r3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2)) + (r1^3*r3^2 - 2*r1^3*r3*r4 + r1^3*r4^2 + r1^2*r3^3 + 6*r1^2*r3^2*r4 - r1^2*r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) - 10*r1^2*r3*r4^2 + r1^2*r4^3 + r1^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + 4*r1*r3^3*r4 - r1*r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) + 6*r1*r3^2*r4^2 - r1*r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - 2*r1*r3*r4^3 + r3^3*r4^2 - r3^3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + r3^2*r4^3 - r3^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2))^2)/(4*r1*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4)) + (-r1^3*r3^2 + 2*r1^3*r3*r4 - r1^3*r4^2 - r1^2*r3^3 - 6*r1^2*r3^2*r4 + r1^2*r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + 10*r1^2*r3*r4^2 - r1^2*r4^3 - r1^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2) - 4*r1*r3^3*r4 + r1*r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) - 6*r1*r3^2*r4^2 + r1*r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + 2*r1*r3*r4^3 - r3^3*r4^2 + r3^3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^2*r4^3 + r3^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2))/(4*r1*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4))

res[4][2] |> println

   (-r1^3*r3^2 + r1^3*r3*r4 - 2*r1^3*r4^2 - r1^2*r3^3 - 5*r1^2*r3^2*r4 + r1^2*r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + 2*r1^2*r3*r4^2 + r1^2*r3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - 2*r1^2*r4^3 - 2*r1*r3^3*r4 + r1*r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) - 5*r1*r3^2*r4^2 + 2*r1*r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + r1*r3*r4^3 + r1*r3*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^3*r4^2 + r3^3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^2*r4^3 + r3^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2))/(8*r1*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4))

res[4][3] |> println

   (3*r1^3*r3 - r1^3*r4 + 9*r1^2*r3*r4 - r1^2*r3*sqrt(r1^2 + 6*r1*r4 + r4^2) + r1^2*r4^2 - r1^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - 3*r1*r3^3 - 7*r1*r3^2*r4 - r3^3*r4 + r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^2*r4^2 + r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2))/(8*r1*(r1*r3 - r1*r4 + r3^2 + r3*r4))

res[4][4] |> println
res[4][4](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println

   -(r1 + r3)*(r3 - r4)*(r3 + r4)*sqrt(r1^2 + 6*r1*r4 + r4^2)/(8*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4)) + (r3 - r4)*(r1^2*r3 + r1^2*r4 + r1*r3^2 + 8*r1*r3*r4 - r1*r4^2 + 3*r3^2*r4 + 3*r3*r4^2)/(8*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4))

res[4][1](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
res[4][2](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
res[4][3](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
res[4][4](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println

   10.0000000000000
   2.70833333333333
   0.666666666666667
   -1.37500000000000

まえもって r1, r3, r4 に定数を代入しておいてから解くと数値解が求まる。

@syms r0::positive, x0::negative, y0::positive,
     r1::positive, r2::positive, r3::positive, r4::positive
@syms r0, x0, y0, r1, r2, r3, r4
(r1, r3, r4) = (10, 6, 15) .// 2
eq1 = (r1 - x0)^2 + (r1 - y0)^2 - (r0 + r1)^2
eq2 = (r2 - x0)^2 + (-r2 - y0)^2 - (r0 + r2)^2
eq3 = (-r3 - x0)^2 + (r3 - y0)^2 - (r0 + r3)^2
eq4 = (-r4 - x0)^2 + (-r4 - y0)^2 - (r0 + r4)^2
res = solve([eq1, eq2, eq3, eq4], (r2, r0, x0, y0))

   4-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (-18, -305/8, 17, -207/8)
    (55/2, -305/8, 17, -207/8)
    (-1/2, 65/24, 2/3, -11/8)
    (10, 65/24, 2/3, -11/8)

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3, r4) = (10, 6, 15) .// 2
   (r2, r0, x0, y0) = [10, 65/24, 2/3, -11/8]
   plot()
   circle(r1, r1, r1)
   circle(r2, -r2, r2, :blue)
   circle(-r3, r3, r3, :green)
   circle(-r4, -r4, r4, :magenta)
   circle(x0, y0, r0, :orange)
   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(x0, y0, "容円:r0,(x0,y0)", :black, :left, delta=-delta/2)
       point(r1, r1, "坤円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(r2, -r2, "乾円:r2,(r2,-r2)", :blue, :center, delta=-delta/2)
       point(-r3, r3, "巽円:r3\n(-r3,r3)", :green, :center, :bottom, delta=delta/2)
       point(-r4, -r4, "艮円:r4\n(-r4,-r4)", :magenta, :center, :bottom, delta=delta/2)
   end
end;

 

算額（その789）

算額（その789）

藤田貞資門人 東都 八木林平質 文化三年丙寅正月
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

十四 群馬県佐波郡玉村町樋越 神明宮 文化3年(1806)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

キーワード：円1個，円弧1個，直角三角形
#Julia, #SymPy, #算額, #和算

直角三角形の中に容円と 2 つの頂点を通る弧が入っている。容円は直角三角形の 2辺と接し，弧とも接している。直角を挟む二辺の短い方（鈎）と長い方（股）の長さがそれぞれ 9 寸，12 寸のとき，容円が最小であるときの直径を求めよ。

少し考えると，容円は大きくなる一方である。最小というのは0ではないが，暗黙のうちに，「弧と容円が図のような配置のとき」という条件がある。つまり，「弧は斜辺と 2 点で交わってはいけない」ということであり，斜辺と1点で交わるときは弧が薄くなっていけば容円は大きくなり続ける。弧が一番厚いのは，斜辺が弧の接線のときである。

鈎，股を a, b
弧がその一部である円の半径と中心座標を r1, (b/2, y1)
容円の半径と中心座標を r2, (r2, y2)
とおく。
斜辺が弧の接線であるということは，(b/2)/y1 = a/b ということである。
以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a, b, r1::positive, y1::negative, r2::positive, y2::positive, d
(a, b) = (Sym(9), Sym(12))
y1 = -b^2/2a
eq1 = (b/2 - r2)^2 + (y2 - y1)^2 - (r1 + r2)^2
eq2 = b^2/4 + y1^2 - r1^2
eq3 = (sqrt(a^2 + b^2) +a)*r2 + y2*b - a*b
res = solve([eq1, eq2, eq3], (r1, r2, y2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (10, 5/2, 4)

鈎，股の長さがそれぞれ 9 寸，12 寸のとき条件を満たすときの容円の直径の最小値は 5 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (9, 12)
   y1 = -b^2/2a
   (r1, r2, y2) = (10, 5/2, 4)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
   circle(r2, y2, r2)
   θ = atand(-y1, b/2)
   circle(b/2, y1, r1, :blue, lw=0.2)
   circle(b/2, y1, r1, :blue, beginangle=θ, endangle=180 - θ, lw=1)
   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(b/2, y1, "(b/2,y1)", :blue, :center, delta=-delta/2)
       point(r2, y2, " 容円:r2,(r2,y2)", :red, :left, :vcenter)
       point(b, 0, " b", :black, :left, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
       plot!(ylims=(1.15y1, 10))
   end
end;

算額（その788）

算額（その788）

享和3年癸亥五月 丸山良玄門人 豫州松山 大西佐兵衛義全
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

平田浩一：愛媛の和算と算額，日本数学会秋季総合分科会市民講演会，2013/9/23.
https://www.mathsoc.jp/publication/tushin/1804/1804hirata.pdf

弧環減球に中球 2 個と小球 1 個を除いた体積を求めよ。

弧環減球は円を両側から 2 円でくり抜き，それを回転させて得られる回転体である。図の太い黒で示した図形を y 軸を中心として回転する。上の図は，それを三次元表示したものである。

曲線が切り替わる点 (x0, y0) を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r1, r2, r3, x0, y0
(r2, r3) = (3, 2) .// 2
R = 2r2 + r3
eq1 = (r1 + r3)^2 + (r3 + r2)^2 - (r1 + r2)^2
eq2 = (x0 - r1 - r3)^2 + y0^2 - r1^2
eq3 = x0^2 + y0^2 - R^2
res = solve((eq1, eq2, eq3), (r1, x0, y0))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (5, 9/4, -5*sqrt(7)/4)
    (5, 9/4, 5*sqrt(7)/4)

[0, 5√7/4] で回転体の体積 res1 を求める。

@syms x
res1= integrate(PI*(-sqrt(25 - x^2) + 6)^2, (x, 0, res[2][3]))
res1 |> println
res1.evalf() |> println

   pi*(-150*asin(sqrt(7)/4) + 8365*sqrt(7)/192)
   21.5487629153420

[5√7/4, 4]で回転体の体積 res2 を求める。

res2 = integrate(PI*sqrt(16 - x^2)^2, (x, res[2][3], 4))
res2 |> println
res2.evalf() |> println

   -2965*sqrt(7)*pi/192 + 128*pi/3
   5.68345793167528

3 個の球の体積の和 res3 を求める。

res3 = 4PI/3*(2*(3//2)^3 + (2//2)^3)
res3 |> println
res3.evalf() |> println

   31*pi/3
   32.4631240870945

求める立体の体積 res を求める。

res = 2res1 + 2res2 - res3 
res |> println
res |> simplify |> println
res.evalf() |> println

   -2965*sqrt(7)*pi/96 + 2*pi*(-150*asin(sqrt(7)/4) + 8365*sqrt(7)/192) + 75*pi
   75*pi*(-16*asin(sqrt(7)/4) + 4 + 3*sqrt(7))/4
   22.0013176069401

求める体積は，75π(4 + 3√7 - 16asin(√7/4))/4 = 22.001317606940137 である。

75π*(4 + 3√7 - 16asin(√7/4))/4

   22.001317606940137

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (3, 2) .// 2
   r1 = 5
   R = 2r2 + r3
   println("R = $R")
   (r1, x0, y0) = (5, 9/4, 5*sqrt(7)/4)
   θ1 = atand(y0, x0)
   θ2 = atand(y0, r1 + r3 - x0)
   plot()
   circle(0, 0, R, :palevioletred1)
   circle(0, 0, R, :black, beginangle=θ1, endangle=90, lw=2)
   circle(0, 0, R, :black, beginangle=270, endangle=360-θ1, lw=2)
   circle22(0, r3 + r2, r2, :magenta)
   circle(0, 0, r3, :orange)
   circle(r1 + r3, 0, r1, :blue)
   circle(r1 + r3, 0, r1, :black, beginangle=180-θ2, endangle=180+θ2, lw=2)
   segment(0, r3 + 2r2, 0, -r3 - 2r2, :black, lw=2)
   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(r1 + r3, 0, "r1+r3")
       point(0, r3, " r3", :black, :left, :bottom, delta=delta)
       point(0, r3+ r2, " r3+r2", :magenta, :left, :vcenter)
       point(x0, y0, " (x0,y0)", :black, :left, :vcenter, mark=false)
   end
end;