算額（その1440）
四 岩手県花巻市南笹間 東光寺 慶応2年(1866)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
寒河江市平塩 熊野神社 明治2年(1869)
徳竹亜紀子，谷垣美保，萬 伸介：天童市山元の若松寺観音堂算額の調査報告 ― 服部武右衛門が関わる算額 ー，仙台高等専門学校名取キャンパス 研究紀要，第 59 号，p.1-8，2023.
https://www.jstage.jst.go.jp/article/nitsendainatori/59/0/59_1/_pdf/-char/ja
キーワード：円1個，外円，正方形4個
#Julia, #SymPy, #算額, #和算
円の中に大正方形 1 個，小正方形 3 個を容れる。大正方形の一辺の長さが 1 寸のとき，小正方形の一辺の長さはいかほどか。
円の半径と中心座標を R, (0, 0)
大正方形の一辺の長さを a
小正方形の一辺の長さを b
正方形の頂点座標を (ax, ay), ..., (fx, fy)
とおき，方程式を解く。
座標点のうちいくつかは特定できる。結局，求めるべきパラメータは b, R, cx である。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, a::positive, b::positive,
ax::positive, ay::positive, bx::positive, by::positive,
cx::positive, cy::positive, dx::positive, dy::positive,
ex::positive, ey::positive, fx::positive, fy::positive
ax = a/2
ay = sqrt(R^2 - ax^2)
bx = ax
by = ay - a
cy = -sqrt(R^2 - cx^2)
dx = a//2 ### 数値解から得られた関係式
dy = -sqrt(R^2 - dx^2)
ex = b/2
ey = -a//2 ### 数値解から得られた関係式
fx = ex
fy = ey - b
eq1 = (cx - dx)^2 + (cy - dy)^2 - b^2 |> expand
eq2 = (ex - bx)^2 + (ey - by)^2 - b^2 |> expand
eq3 = (bx - dx)^2 + (by - dy)^2 - 2b^2 |> expand;
res = solve([eq1, eq2, eq3], (b, R, cx))[4] # 4 of 4
(a*sqrt(-3*sqrt(2) - sqrt(11 - 6*sqrt(2)) + 6), a*sqrt(12 - 6*sqrt(2))/2, a*(7*sqrt(34 - 24*sqrt(2)) + sqrt(22 - 12*sqrt(2)) + 14*sqrt(17 - 12*sqrt(2)) + 2*sqrt(11 - 6*sqrt(2)))/12)
# b 二重根号を外す
res[1] |> sympy.sqrtdenest |> println
a*(-1 + sqrt(2))
# R
res[2] |> println
a*sqrt(12 - 6*sqrt(2))/2
# cx 二重根号を外す
res[3] |> sympy.sqrtdenest |> simplify |> println
a*(3 - sqrt(2))/2
小正方形の一辺の長さ b は，大正方形の一辺の長さ a の (√2 - 1) 倍である。
大正方形の一辺の長さが 1 寸のとき，小正方形の一辺の長さは (√2 - 1) = 0.41421356237309515 寸である。
function draw(a, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
b = a*(√2 - 1)
R = a*sqrt(12 - 6√2)/2
cx = a*(3 - √2)/2
ax = a/2
ay = sqrt(R^2 - ax^2)
bx = ax
by = ay - a
cy = -sqrt(R^2 - cx^2)
dx = a//2 ### 数値解から得られた関係式
dy = -sqrt(R^2 - dx^2)
ex = b/2
ey = -a//2 ### 数値解から得られた関係式
fx = ex
fy = ey - b
@printf("a = %g; b = %g; R = %g; cx = %g\n", a, b, R, cx)
plot([bx, ax, -ax, -bx, bx], [by, ay, ay, by, by], color=:blue, lw=0.5)
circle(0, 0, R)
plot!([cx, bx, ex, dx, cx], [cy, by, ey, dy, cy], color=:green, lw=0.5)
plot!(-[cx, bx, ex, dx, cx], [cy, by, ey, dy, cy], color=:green, lw=0.5)
plot!([fx, ex, -ex, -fx, fx], [fy, ey, ey, fy, fy], color=:green, lw=0.5)
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(ax, ay, " (ax,ay)", :green, :left, :vcenter)
point(bx, by, " (bx,by)", :green, :left, :vcenter)
point(cx, cy, " (cx,cy)", :green, :left, :vcenter)
point(dx, dy, " (dx,dy)", :green, :left, :vcenter)
point(ex, ey, " (ex,ey)", :green, :left, :vcenter)
point(fx, fy, " (fx,fy)", :green, :left)
point(0, R, "R", :red, :center, :bottom, delta=delta/2)
end
end;
draw(1, true)
