算額（その1375）

算額（その1375）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円1個，外円，正五角形
#Julia, #SymPy, #算額, #和算

外円の中に 3 個の合同な正五角形を入れる。外円の直径が 10 寸のとき，正五角形の一辺の長さはいかほどか。

問題文では，「円の中に 3 個の正五角形を容れる」と書いているが，3 個ではなく位置と方向を明示すれば「1 個を容れる」だけで十分である。もちろん図のように「10 個を容れる」としてもよいが，それではヒントを与えすぎると考えたのだろうか。

ここでは，青で描いた上部の正五角形について考える。赤で描いた外円の半径と中心座標を R, (0, 0) とする。正五角形が内接する灰色で描いた円の半径を AB = BD = r とする。
与えられるのは OA = R，求めるのは x = CD である。
∠CBD = 36°，∠COD = 18° である。
BC = r*cos(36°)
CD = r*sin(36°) より，r = CD/sin(36°) = x/sin(36°)
OC*tan(18°) = CD
OA = AB + BC + OC
R = r + r*cos(36°) + x/tan(18°)
R = x/sin(36°)*(1 + cos(36°)) + x/tan(18°)
これを x について解いて
x = 2R*sqrt(35 - 15√5)/(√10sqrt(5 - 2√5) + 5√2*sqrt(5 - 2√5) + 2√5*sqrt(5 - √5))
を得る。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r, x
eq = x/sind(Sym(36))*(1 + cosd(Sym(36))) + x/tand(Sym(18)) - R;
ans_x = solve(eq, x)[1]
ans_x |> println

   2*R*sqrt(35 - 15*sqrt(5))/(sqrt(10)*sqrt(5 - 2*sqrt(5)) + 5*sqrt(2)*sqrt(5 - 2*sqrt(5)) + 2*sqrt(5)*sqrt(5 - sqrt(5)))

SymPy では変数が残っている式は簡約化されにくいという特徴（欠点）を持っているので，R の係数部分（下記）だけを簡約化する。

ans_x/R |> println

   2*sqrt(35 - 15*sqrt(5))/(sqrt(10)*sqrt(5 - 2*sqrt(5)) + 5*sqrt(2)*sqrt(5 - 2*sqrt(5)) + 2*sqrt(5)*sqrt(5 - sqrt(5)))

@syms d
apart(ans_x/R, d) |> simplify |> println

   sqrt(25 - 10*sqrt(5))/10

長い式は sqrt(25 - 10√5)/10 に簡約化される。
術では sqrt(0.25 - sqrt(0.05)) としているが 1/10 を外側の平方根の中に入れると同じである。

正五角形の一辺の長さ 2x は 外円の半径が R のとき　2R*sqrt(25 - 10√5)/10 である。
外円の直径が 10 寸のとき，正五角形の一辺の長さは sqrt(25 - 10√5) = 1.6245984811645313 寸である。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # x = 2R*sqrt(35 - 15√5)/(√10sqrt(5 - 2√5) + 5√2*sqrt(5 - 2√5) + 2√5*sqrt(5 - √5))
   x = R*sqrt(25 - 10√5)/10
   r = x/sind(36)
   θ = 18:72:378
   xa = @. r*cosd(θ)
   ya = @. r*sind(θ) + (R - r)
   plot([-r*sind(36), 0, r*sind(36), 0], [R - r*(1 + cosd(36)), 0, R - r*(1 + cosd(36)), R - r], color=:gray70, lw=0.5)
   circle(0, 0, R, :red)
   circle(0, R - r, r, :gray70)
   for i = 1:10
       plot!(xa, ya, color=:green, lw=1)
       (xa, ya) = transform(xa, ya, deg=36)
   end
   plot!(xa, ya, color=:blue, lw=1)
   if more == true
       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, R, " A", :blue, :left, :bottom, delta=delta/2)
       point(0, R - r, " B", :blue, :left, :bottom, delta=delta/2)
       point(0, R - r*(1 + cosd(36)), " C", :blue, :left, :bottom, delta=delta/2)
       point(r*sind(36), R - r*(1 + cosd(36)), " D", :blue, :left, :bottom, delta=delta/2)
       point(0, 0, " O", :blue, :left, :bottom, delta=delta/2)
   end
end;

draw(10/2, true)