算額（その770）
神壁算法 天明8年戊申 二月 關流藤田貞資門人 細川能登守家士 渡邉金八源耀
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
正五角形に5本の斜線を入れ，内部に小さな正五角形を作り，できた 6 つの区画にそれぞれ同じ直径を持つ円を内接させる。
正五角形の一辺の長さが 533 寸 6 分のとき，等円の直径はいかほどか。
正五角形の一辺の長さを L，正五角形が内接する円の半径を R とする。
内側の正五角形の一辺の長さを l とする。
斜線と正五角形の一辺で作られる三角形 BAR において第二余弦定理を使う。
AR, AB, BR の長さを AR, AB, BR，また，A, B, R の座標を (x, y), (x0, y0), (0, R) とする。
ポイントは，斜線は外側の五角形の一辺と平行ではないということである。
include("julia-source.txt");
using SymPy
@syms L::positive, R::positive,
x0::positive, y0::positive,
x1::positive, y1::positive,
x::negative, y::positive,
AR::positive, l::positive, r2::positive
d18 = Sym(18)
R = L/2sind(2d18)
x0 = R*cosd(d18)
y0 = R*sind(d18)
eq1 = (x0 - x)^2 + (y0 - y)^2 - (AR + l)^2
eq2 = x^2 + (R - y)^2 - AR^2
eq3 = AR^2 + (AR + l)^2 - 2AR*(AR + l)*cosd(4d18) - L^2
eq4 = dist(x0, y0, R*cosd(-3d18), R*sind(-3d18), x1, y1) - r2^2
eq5 = dist(x0, y0, x, y, x1, y1) - r2^2
eq6 = dist(x0, y0, x, y, 0, 0) - r2^2
eq7 = dist(c*cosd(atand(y, x) - 4d18), c*sind(atand(y, x) - 4d18), R*cosd(-3d18), R*sind(-3d18), x1, y1) - 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 Float64.(v), r.f_converged
end;
function H(par)
(AR, l, x, y, x1, y1, r2) = par
a = 2sqrt(5/8 - √5/8)
b = √5/4 - 1/4
g = L*sqrt(5/8 + √5/8)
c = x - g/a
d = L/2 - g/a
e = x1 - g/a
f = (atan(y/x)/pi + 1) - 2/5
h = -L/a*((√5 + 1)/4 + b)
i = sqrt(x^2 + y^2)
j = i*cos(pi*f)
k = i*sin(pi*f)
m = -L*(√5 + 1)/4a - k
n = L/2 - j
o = (x1 - j)*n + (y1 - k)*m
p = y - L*b/a
q = y1 - L*b/a
r = h^2 + d^2
s = p^2 + c^2
t = n^2 + m^2
u = (L*b*p + g*c)/(s*a)
return [
(L*b/a - y)^2 + (g/a - x)^2 - (AR + l)^2, # eq1
x^2 + (L/a - y)^2 - AR^2, # eq2
AR^2 - 2AR*b*(AR + l) - L^2 + (AR + l)^2, # eq3
(y1 - L*b/a - (q*h + d*e)*h/r)^2 + (x1 - g/a - (q*h + d*e)*d/r)^2 - r2^2, # eq4
(y1 - L*b/a - (p*q + c*e)*p/s)^2 + (x1 - g/a - (p*q + c*e)*c/s)^2 - r2^2, # eq5
(p*u - L*b/a)^2 + (c*u - g/a)^2 - r2^2, # eq6
(x1 - n*o/t - j)^2 + (y1 - o*m/t - k)^2 - r2^2, # eq7
]
end;
L = 5336//10
iniv = BigFloat[347.1, 179.5, -94.4, 119.9, 259.5, 10.0, 123.5]
res = nls(H, ini=iniv)
([347.0551765561487, 179.4560303257959, -94.42748155686841, 119.94507841468496, 259.51771893633804, 10.023199995813941, 123.50001782997607], true)
正五角形の一辺の長さが 533 寸 6 分のとき，等円の直径は 247 寸である。
その他のパラメータは以下のとおりである。
AR = 347.055; x = -94.4275; y = 119.945; x1 = 259.518; y1 = 10.0232; r2 = 123.5
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
L = 5336//10
#l = 100
R = L/2sind(2*18)
x0 = R*cosd(18)
y0 = R*sind(18)
(AR, l, x, y, x1, y1, r1) = res[1]
@printf("等円の直径 = %g\n", 2r1)
@printf("AR = %g; x = %g; y = %g; x1 = %g; y1 = %g; r2 = %g\n", AR, x, y, x1, y1, r1)
r = sqrt(x^2 + y^2)
θ = atand(y, x)
ix = [r*cosd(72i + θ) for i in 0:5]
iy = [r*sind(72i + θ) for i in 0:5]
ox = [R*cosd(72i + 18) for i in 0:5]
oy = [R*sind(72i + 18) for i in 0:5]
plot()
circle(0, 0, R, :gray80)
circle(0, 0, r, :gray90)
circle(0, 0, r1)
rotate(x1, y1, r1, angle=72)
for i in 1:5
segment(ox[i], oy[i], ix[i], iy[i], :blue)
segment(ox[i], oy[i], ox[i + 1], oy[i + 1], :blue)
end
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(x, y, "A:(x,y) ", :blue, :right, :vcenter)
point(x0, y0, " B:(x0,y0)", :blue, :left, :vcenter)
point(x1, y1, "等円:r1\n(x1,y1)", :red, :center, :bottom, delta=delta)
point(0, R, " R:(0,R)", :black, :left, :bottom, delta=delta)
point(r, 0, " r", :gray50, :left, :bottom, delta=delta)
end
end;
