算額（その504）
高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8
算額（その164)
https://blog.goo.ne.jp/r-de-r/e/ad36c91a7122981d5a236e8addbe3cff
と基本的に同じであるが，乙円，丙円の半径を未知数として一般解を求める。
二等辺三角形内に甲円，乙円，丙円を入れる。乙円と丙円の直径が与えられたとき，丙円の直径を求めよ。
一般性を失わずに，丙円の直径を 1 とする。乙円と丙円の直径が与えられたとき，乙円の直径は「乙円の直径/丙円の直径」，丙円の直径は 1 として与える。結果は「丙円の直径/乙円の直径」倍になる。
二等辺三角形の底辺の長さを 2a，高さを b とする
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (0, 2r1 + r2)
丙円の半径と中心座標を r3, (x3, r3)
として，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a, b, r1, r2, r3, x3
r3 = 1
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = r3/(a - x3) - r1/a
eq3 = r2/(b - 2r1 - r2) - r1/(b - r1)
eq4 = r1/(b - r1) - a/sqrt(a^2 + b^2)
res = solve([eq1, eq2, eq3, eq4], (a, b, r1, x3))
4-element Vector{NTuple{4, Sym}}:
(-(1 - 2*sqrt(r2))^(3/2)/sqrt(r2), 2*(4*sqrt(r2) - 4*r2 - 1)/(2*sqrt(r2) + r2 - 1), 1 - 2*sqrt(r2), 2*sqrt(1 - 2*sqrt(r2)))
((1 - 2*sqrt(r2))^(3/2)/sqrt(r2), 2*(4*sqrt(r2) - 4*r2 - 1)/(2*sqrt(r2) + r2 - 1), 1 - 2*sqrt(r2), -2*sqrt(1 - 2*sqrt(r2)))
(-(2*sqrt(r2) + 1)^(3/2)/sqrt(r2), 2*(2*sqrt(r2) + 1)^2/(2*sqrt(r2) - r2 + 1), 2*sqrt(r2) + 1, -2*sqrt(2*sqrt(r2) + 1))
((2*sqrt(r2) + 1)^(3/2)/sqrt(r2), 2*(2*sqrt(r2) + 1)^2/(2*sqrt(r2) - r2 + 1), 2*sqrt(r2) + 1, 2*sqrt(2*sqrt(r2) + 1))
4 組の解が得られるが，4番目の解が適解である。
甲円の直径は， 4r3*(sqrt(r2/r3) + 1/2) である。
(r2, r3) = (1, 2)
t = r2/r3
u = 2sqrt(t) + 1
a = r3*(u^(3/2)/sqrt(t))
b = r3*(2u^2/(u - t))
r1 = r3*u
x3 = r3*(2sqrt(u))
@printf("r2 = %g; r3 = %g; a = %g; b = %g; r1 = %g; x3 = %g\n", r2, r3, a, b, r1, x3)
@printf("甲円の直径 = %g\n", 4r3*(sqrt(r2/r3) + 1/2))
r2 = 1; r3 = 2; a = 10.6098; b = 12.1793; r1 = 4.82843; x3 = 6.2151
甲円の直径 = 9.65685
乙円，丙円の半径を 1, 2 としたとき，甲円の直径は 9.65685 である。
using Plots
function draw(r2, r3, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
t = r2/r3
u = 2sqrt(t) + 1
a = r3*(u^(3/2)/sqrt(t))
b = r3*(2u^2/(u - t))
r1 = r3*u
x3 = r3*(2sqrt(u))
@printf("r2 = %g; r3 = %g; a = %g; b = %g; r1 = %g; x3 = %g\n", r2, r3, a, b, r1, x3)
@printf("甲円の直径 = %g\n", 4r3*(sqrt(r2/r3) + 1/2))
plot([a, 0, -a, 0], [0, b, 0, 0,], color=:black, lw=0.5)
circle(0, r1, r1)
circle(0, 2r1 + r2, r2, :blue)
circle(x3, r3, r3, :green)
circle(-x3, r3, 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(0, r1, " 甲円:r1,(0,r1)", :magenta, :left, :vcenter)
point(0, 2r1 + r2, " 乙円:r2,(0,2r1+r2)", :magenta, :left, :vcenter)
point(x3, r3, " 丙円:r3,(x3,r3)", :magenta, :center, :top, delta=-delta)
else
plot!(showaxis=false)
end
end;
r2 = 1; r3 = 4; a = 22.6274; b = 18.2857; r1 = 8; x3 = 11.3137
甲円の直径 = 16
