算額（その297）
中村信弥「算額への招待 追補2」
http://www.wasan.jp/syotai/syotai.html
長野県千曲市八幡 武水別神社（八幡神社） 天保8年(1837)
5個の正方形と1個の円が直角三角形に内接している。円の直径が 7 寸，一番小さい正方形の一辺の長さが 6 寸のとき，鉤，股，弦の長さを求めよ。
(復元）算額の図が小さく（あるいは誤解？）分かりづらいが，円と一番小さい正方形の関係は，以下の拡大図のようになっている（算額では円が正方形に内接しているように描かれている）。
また，「算額への招待 追補2」は算額の術の解説をしており，上に示す 1 つの解のみを提示しているが，解は二通りある。もう一つの解の拡大図は以下のとおりである。
記号の説明のために全体図も二通り示しておく。
この 2 つの解は，図を裏返してみれば相似であることがわかる。
何れにせよ，SymPy では二通りの解を与える。
また，図に現れる直角三角形は全て相似なので，x3, y3 以降は SymPy で多元連立方程式を解くことでも得られるが，逐次的に繰り返して求めることができる。直角三角形の頂点 (a, b) と弦の長さ c も計算によって求める。
図に示すように記号を定め，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms x1::positive, x2::positive, y2::positive, r1::positive;
r1 = 7 // 2
x2 = x1 + 6
eq1 = (x2 - x1) / x1 - y2 / x2
eq2 = x2 + y2 - sqrt(x2^2 + y2^2) - 2r1
res = solve([eq1, eq2], (x1, y2))
2-element Vector{Tuple{Sym, Sym}}:
(9/2, 14)
(8, 21/2)
● x1, y2 が (9/2, 14) の場合の解
x1 = 9/2; x2 = 21/2; y2 = 14
鉤 = 1023.63621399177
股 = 767.727160493827
弦 = 1279.54526748971
● x1, y2 が (8, 21/2) の場合の解
x1 = 8; x2 = 14; y2 = 21/2
鉤 = 182.185253906250
股 = 242.913671875000
弦 = 303.642089843750
using Plots
function draw(zoomin=false, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(s, r) = (6, 7/2)
(x1, y2) = res[2]
println((x1, y2))
x2 = x1 + s
tanθ = y2/x2
sinθ = y2/sqrt(x2^2 + y2^2)
cosθ = x2/sqrt(x2^2 + y2^2)
xs = zeros(7)
ys = zeros(7)
xs[1:2] = [x1, x2]
ys[1:2] = [s, y2]
for i = 3:7
xs[i] = xs[i - 1] + ys[i - 1]
ys[i] = xs[i]*tanθ
end
println("ys[6]= $(ys[6])")
c = xs[7] + ys[6] * tanθ
a = xs[6] + ys[6]*cosθ^2
b = ys[6] + ys[6]*cosθ*sinθ
println("x1 = $x1; x2 = $(x1 + 6); y2 = $y2")
println("鉤 = $(sqrt(b^2 + (c - a)^2))")
println("股 = $(sqrt(b^2 + a^2))")
println("弦 = $c")
if zoomin
plot()
else
plot([0, c, a, 0], [0, 0, b, 0], color=:black, lw=0.5)
end
circle(xs[2] - r, r, r)
for i = 1:6
println((xs[i], ys[i], xs[i+1], 0))
rect(xs[i], 0, xs[i+1], ys[i], :blue)
end
abline(0, 0, tanθ, 0, 60)
if more
if !zoomin
point(c, 0, " c", :black, :left, :bottom)
point(a, b, "(a,b) ", :black, :right, :bottom)
end
point(x1, 0, "x1 ", :blue, :right, :bottom)
point(x1, s, "(x1,s) ", :blue, :right, :bottom)
point(x2-r, r, "(x2-r,r)")
point(x2, 0, " x2", :blue, :left, :bottom)
point(x2, y2, "(x2,y2) ", :blue, :right, :bottom)
point(xs[3], 0, " x3", :blue, :left, :bottom)
point(xs[3], ys[3], "(x3,y3) ", :blue, :right, :bottom)
if zoomin
txt = @sprintf(" x1 = %.1f\n x2 = %.1f\n y2 = %.1f", x1, x1 + 6, y2)
annotate!(xs[3], ys[3]/2, text(txt, 10, :left))
else
point(xs[4], 0, " x4", :blue, :left, :bottom)
point(xs[4], ys[4], "(x4,y4) ", :blue, :right, :bottom)
point(xs[5], 0, " x5", :blue, :left, :bottom)
point(xs[5], ys[5], "(x5,y5) ", :blue, :right, :bottom)
point(xs[6], 0, " x6", :blue, :left, :bottom)
point(xs[6], ys[6], "(x6,y6) ", :blue, :right, :bottom)
point(xs[7], 0, " x7", :blue, :left, :bottom)
end
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
zoomin && plot!(showaxis=true, xlims=(0, 47), ylims=(0, 35))
else
plot!(showaxis=false)
end
end;
