算額（その980）
一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
正三角形の中に斜線 2 本，大円 1 個，小円 3 個を入れる。大円の直径が 500 寸のとき，小円の直径はいかほどか。
正三角形の一辺の長さを 2a
斜線と正三角形の底辺との交点座標を (b, 0)
大円の半径と中心座標を r1, (0, 2r1 + r1)
小円の半径と中心座標を r2, (0, r2), (x2, y2)
とおき，手動で連立方程式を解いていったが，最後の最後で，次数が高すぎて SymPy の手に負えなくなったので，その時点で r1 = 1 を代入し数値解を求めた。問の「大円の直径が 500 寸」に対応するには，単に得られた数値ベクトルを 250 倍すればよいだけではある。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive,
r1::positive, r2::positive, x2::positive, y2::positive
x2 = √Sym(3)*(y2 - 2r2 - r1)
eq1 = x2^2 + (y2 - 2r2 - r1)^2 - (r1 - r2)^2
eq2 = 2r1 - (√Sym(3)*a - 2r2 - r1)
eq3 = r2*sqrt(3a^2 + b^2)- b*(√Sym(3)*a - r2)
eq4 = dist2(0, √Sym(3)*a, b, 0, x2, y2, r2);
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, r2, x2, y2))
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
-(r1 - r2)^2 + 4*(-r1 - 2*r2 + y2)^2, # eq1
-sqrt(3)*a + 3*r1 + 2*r2, # eq2
-b*(sqrt(3)*a - r2) + r2*sqrt(3*a^2 + b^2), # eq3
(3*a^2*b^2 + 6*sqrt(3)*a^2*b*r1 + 12*sqrt(3)*a^2*b*r2 - 6*sqrt(3)*a^2*b*y2 + 9*a^2*r1^2 + 36*a^2*r1*r2 - 18*a^2*r1*y2 + 33*a^2*r2^2 - 36*a^2*r2*y2 + 9*a^2*y2^2 - 2*sqrt(3)*a*b^2*y2 - 6*a*b*r1*y2 - 12*a*b*r2*y2 + 6*a*b*y2^2 - b^2*r2^2 + b^2*y2^2)/(3*a^2 + b^2), # eq4
ans_y2 = solve(eq1, y2)[2]
ans_y2 |> println
3*r1/2 + 3*r2/2
ans_a = solve(eq2, a)[1]
ans_a |> println
sqrt(3)*(r1 + 2*r2/3)
eq3 = eq3(a => ans_a) |> simplify;
ans_b = solve(eq3, b)[1]
ans_b |> println
r2*sqrt(9*r1 + 6*r2)/(3*sqrt(r1))
ans_b |> simplify
ans_b |> println
r2*sqrt(9*r1 + 6*r2)/(3*sqrt(r1))
eq4 = eq4(y2 => ans_y2, a => ans_a, b => ans_b) |> simplify |> numerator |> simplify;
eq4 |> println
-54*r1^(7/2)*r2*sqrt(3*r1 + 2*r2) + 42*r1^(3/2)*r2^3*sqrt(3*r1 + 2*r2) + 12*sqrt(r1)*r2^4*sqrt(3*r1 + 2*r2) + 81*r1^5 - 54*r1^4*r2 - 180*r1^3*r2^2 - 72*r1^2*r2^3 - 9*r1*r2^4 - 6*r2^5
eq4 = eq4(r1 => 1)
eq4 |> println
-6*r2^5 + 12*r2^4*sqrt(2*r2 + 3) - 9*r2^4 + 42*r2^3*sqrt(2*r2 + 3) - 72*r2^3 - 180*r2^2 - 54*r2*sqrt(2*r2 + 3) - 54*r2 + 81
ans_r2 = solve(eq4, r2)[2]
ans_r2 |> println
-1 - 4/(3*(4 + 4*sqrt(93)/9)^(1/3)) + (4 + 4*sqrt(93)/9)^(1/3)
r1 = 1
r2 = -1 - 4/(3*(4 + 4*sqrt(93)/9)^(1/3)) + (4 + 4*sqrt(93)/9)^(1/3)
b = r2*sqrt(9r1 + 6r2)/(3*sqrt(r1))
a = √3*(r1 + 2*r2/3)
y2 = 3r1/2 + 3r2/2
x2 = √3*(y2 - 2r2 - r1)
(r1, r2, a, b, x2, y2) |> println
250 .*(r1, r2, a, b, x2, y2) |> println
(1, 0.36465560765603855, 2.153118834052318, 0.4065711730355355, 0.5502243839218582, 2.046983411484058)
(250, 91.16390191400964, 538.2797085130795, 101.64279325888388, 137.55609598046456, 511.7458528710145)
大円の直径が 500 寸のとき，小円の直径は 182.32780382801928 寸である。
前にも 1 問イレギュラーなのがあったが，このシリーズでは解は整数値に極めて近い解になるのがシバリであるが，この問も解がきれいな数値にならない。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1
r2 = -1 - 4/(3*(4 + 4*sqrt(93)/9)^(1/3)) + (4 + 4*sqrt(93)/9)^(1/3)
b = r2*sqrt(9r1 + 6r2)/(3*sqrt(r1))
a = √3*(r1 + 2*r2/3)
y2 = 3r1/2 + 3r2/2
x2 = √3*(y2 - 2r2 - r1)
(r1, r2, a, b, x2, y2) = 250 .*(r1, r2, a, b, x2, y2)
@printf("大円の直径 = %g; 小円の直径 = %g\n", 2r1, 2r2)
@printf("a = %g; b = %g; r2 = %g; x2 = %g; y2 = %g\n", a, b, r2, x2, y2)
plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
circle(0, 2r2 + r1, r1)
circle(0, r2, r2, :green)
circle2(x2, y2, r2, :green)
segment(b, 0, 0, √3a)
segment(-b, 0, 0, √3a)
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, 2r2 + r1, "大円:r1,(0,2r2+r1)", :red, :center, delta=-delta)
point(x2, y2, "小円:r2\n(x2,y2)", :black, :center, delta=-delta/2)
point(0, r2, "小円:r2\n(0,r2)", :green, :center, delta=-delta/2)
point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
point(b, 0, " b", :blue, :left, :bottom, delta=delta/2)
point(0, √3a, " √3a", :blue, :left, :bottom, delta=delta/2)
end
end;
