算額（その318）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（75）
長野県長野市元善町 善光寺 文化元年(1804)
正八角形の対角線で区切られた三角形に内接する円の直径を求めよ。
問では，間に一つの頂点を挟んだ2つの頂点間の対角線の長さ（距斜）と正八角形の一辺（濶）を与えたときの各円の直径を求めよと書かれているが，この 2 つの量はいずれも外円の半径により求められるので，外円の半径 R を適当に定めれば各円の直径を求めることができる。ここでは，R = 1 とする。
外円の半径と中心座標 R, (0, 0)
甲円の半径と中心座標 r1, (R/√2 + r1, 0)
乙円の半径と中心座標 r2, (R/√2 - r2, y2)
丙円の半径と中心座標 r3, (x3, y3)
丁円の半径と中心座標 r3, (-y3, -x3)
戊円の半径と中心座標 r2, (-y2, -R/√2 + r2)
以下の連立方程式を解く。3 組の連立方程式 [eq1], [eq2, eq3], [eq4, eq5, eq6] は互いに独立。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive;
@syms a, b
R = 1
a = R/sqrt(Sym(2))
b = sqrt((R - R/sqrt(Sym(2)))^2 + (R/sqrt(Sym(2)))^2)
eq1 = distance(R/sqrt(Sym(2)), R/sqrt(Sym(2)), R, 0, R/sqrt(Sym(2)) + r1, 0) - r1^2
res1 = solve(eq1, r1)[1]
res1 |> println
(-sqrt(2 - sqrt(2)) + sqrt(2)*(2 - sqrt(2)))/(2 - sqrt(2))^(3/2)
@syms r2::positive, y2::positive;
eq2 = distance(0, R, R/sqrt(Sym(2)), R/sqrt(Sym(2)), R/sqrt(Sym(2)) - r2, y2) - r2^2
eq3 = distance(0, R, R/sqrt(Sym(2)), -R/sqrt(Sym(2)), R/sqrt(Sym(2)) - r2, y2) - r2^2
res2 = solve([eq2, eq3], (r2, y2))
2-element Vector{Tuple{Sym, Sym}}:
(sqrt(2)*(-sqrt(2)/2 - (1 - sqrt(2)/2)^(3/2) + 1/2 + sqrt(1 - sqrt(2)/2)), sqrt(1 - sqrt(2)/2))
(sqrt(2)*(-(sqrt(2)/2 + 1)^(3/2) + 1/2 + sqrt(2)/2 + sqrt(sqrt(2)/2 + 1)), sqrt(sqrt(2)/2 + 1))
最初の解が適解である。
@syms r3::positive, x3::negative, y3::positive;
eq4 = distance(0, R, R/sqrt(Sym(2)), -R/sqrt(Sym(2)), x3, y3) - r3^2
eq5 = distance(0, R, -R/sqrt(Sym(2)), R/sqrt(Sym(2)), x3, y3) - r3^2
eq6 = distance(R/sqrt(Sym(2)), -R/sqrt(Sym(2)), -R/sqrt(Sym(2)), R/sqrt(Sym(2)), x3, y3) - r3^2
res3 = solve([eq4, eq5, eq6], (r3, x3, y3))
2-element Vector{Tuple{Sym, Sym, Sym}}:
(sqrt(-sqrt(4 - 2*sqrt(2)) - sqrt(2)/2 + 2), sqrt(2)*(-sqrt(2)*sqrt(2 - sqrt(2)) - sqrt(2 - sqrt(2)) + 1)/2, sqrt(1 - sqrt(2)/2) + sqrt(2)/2)
(sqrt(-sqrt(2*sqrt(2) + 4) + sqrt(2)/2 + 2), sqrt(2)*(-sqrt(sqrt(2) + 2) - 1 + sqrt(2)*sqrt(sqrt(2) + 2))/2, -sqrt(2)/2 + sqrt(sqrt(2)/2 + 1))
2番めの解が適解である。
甲円の直径 = 0.281305
乙円（戊円）の直径 = 0.496606
丙円（丁円）の直径 = 0.613126
# 正八角形の頂点座標
R = 1
x = zeros(9)
y = zeros(9)
for i in 1:8
θ = (i - 1)*pi/4
x[i] = cos(θ)R
y[i] = sin(θ)R
end
x[9] = x[1]
y[9] = y[1]
function symmetry(x, y, r, color=:red)
circle(x, y, r, color)
circle(-y, -x, r, color)
end;
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = (-sqrt(2 - sqrt(2)) + sqrt(2)*(2 - sqrt(2)))/(2 - sqrt(2))^(3/2)
(r2, y2) = (sqrt(2)*(-sqrt(2)/2 - (1 - sqrt(2)/2)^(3/2) + 1/2 + sqrt(1 - sqrt(2)/2)), sqrt(1 - sqrt(2)/2))
(r3, x3, y3) = (sqrt(-sqrt(2*sqrt(2) + 4) + sqrt(2)/2 + 2), sqrt(2)*(-sqrt(sqrt(2) + 2) - 1 + sqrt(2)*sqrt(sqrt(2) + 2))/2, -sqrt(2)/2 + sqrt(sqrt(2)/2 + 1))
@printf("甲円の直径 = %.6f\n乙円（戊円）の直径 = %.6f\n丙円（丁円）の直径 = %.6f\n", 2r1, 2r2, 2r3)
plot(x, y, color=:black, lw=0.5)
circle(0, 0, R)
symmetry(R/√2 + r1, 0, r1, :green)
symmetry(R/√2 - r2, y2, r2, :orange)
symmetry(x3, y3, r3, :magenta)
for i in 2:6
segment(x[i], y[i], x[8], y[8], :blue)
end
if more
point(R/√2 + r1, 0, "甲:r1\n(R/√2,0)", :green, :center, :bottom)
point(R/√2 - r2, y2, "乙:r2\n(R/√2-r2,y2)\n", :orange, :center, :bottom)
point(x3, y3, "丙:r3\n(x3,y3)\n", :magenta, :center, :bottom)
point(-y3, -x3, " 丁", :magenta)
point(-y2, -R/√2 + r2, " 戊", :orange)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
