算額（その937）

算額（その937）

千葉県富津市西大和田（旧大佐和町絹） 吾妻神社 文久2年(1862)
「算額」第三集 全国調査，香川県算額研究会

正五角形の中に対角線 2 本を引き，区分された領域に 甲円 2 個，乙円 2 個，丙円 1 個を入れる。甲円，乙円の直径がそれぞれ 1.6180 寸，1.3819 寸のとき，丙円の直径はいかほどか。

座標の定義に使う三角関数の値
(x18, y18) = (cosd(18), sind(18))
(x54, y54) = (cosd(54), sind(54))

正五角形が内接する外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1,(0, R*y18 - r1)
乙円の半径と中心座標を r2, (x2, 2R*y18 - R)
丙円の半径と中心座標を r3, (0, r3 - 2R*y54)
とする。

算額の答えを得るには，実は乙円は関係ない（外円の大きさ R も関係ない）。

2 つの三角形 ABC と DEC は相似で，相似比は R*x18 : R*x54 = cosd(18) : cosd(54) である。

したがって，それに内接する甲円と丙円も相似で，その半径の比も r1:r3 = cosd(18) : cosd(54) なので，
r3 = r1 * cosd(54)/cosd(18)
  = r1 * sqrt(5/8 - sqrt(5)/8) / sqrt(sqrt(5)/8 + 5/8)
  = 1.6180/2 * (-1/2 + sqrt(5)/2)
  = 0.499989496898665

すなわち，甲円の直径が 1.6180 寸のとき，丙円の直径は 0.99997899379733 寸（ほぼほぼ 1 寸）である。

1.6180/2 * (-1/2 + sqrt(5)/2)

   0.499989496898665

以下では，図を描くためにパラメータを得る手順を示す。

include("julia-source.txt");

using SymPy
@syms x18::positive, y18::positive, x54::positive, y54::positive
(x18, y18) = (cosd(Sym(18)), sind(Sym(18)))
(x54, y54) = (cosd(Sym(54)), sind(Sym(54)));

まず，外円の大きさを求める。外円の半径 R　は，甲円の半径 r1 の (1 + √5) 倍である。 

@syms R::positive, r1::positive
eq1 = dist2(0, R, R*x18, R*y18, 0, R*y18 + r1, r1)
R = solve(eq1, R)[1] |> sympy.sqrtdenest |> simplify
R |> println

   r1*(1 + sqrt(5))

次に，丙円の大きさを求める。丙円の半径 r3 は（前にも示したが），甲円の半径 r1 の (√5 - 1)/2 倍である。

@syms r3::positive
eq4 = dist2(R*x18, R*y18, -R*x54, -R*y54, 0, r3 - R*y54, r3);
eq4 = eq4 |> simplify |> numerator
r3 = solve(eq4, r3)[2] |> simplify |> sympy.sqrtdenest |> simplify
r3 |> println

   r1*(-1 + sqrt(5))/2

最後に，乙円の中心の x 座標を求める。

@syms r2::positive, x2::positive, y2::negative
eq2 = dist2(R*x18, R*y18, R*x54, -R*y54, x2, 2R*y18 - R, r2) |> simplify |> numerator;
ans_x2 = solve(eq2, x2)[1] |> simplify;
ans_x2 |> println

   r1*(15*sqrt(2)*sqrt(sqrt(5) + 5) + 11*sqrt(10)*sqrt(5 - sqrt(5)) + 25*sqrt(2)*sqrt(5 - sqrt(5)) + 7*sqrt(10)*sqrt(sqrt(5) + 5))/(2*(11*sqrt(5) + 25)) - 2*sqrt(-850*r1^2 - 380*sqrt(5)*r1^2 + 190*r1^2*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 85*sqrt(5)*r1^2*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 152*sqrt(5)*r2^2 + 340*r2^2)/(11*sqrt(5) + 25)

function draw(more=false)
   pyplot(size=(600, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x18, y18) = (cosd(18), sind(18))
   (x54, y54) = (cosd(54), sind(54))
   r1 = 1.6180/2
   R = r1*(1 + √5)
   r3 = r1*(√5 - 1)/2
   r2 = 1.3819/2
   s = sqrt(5 - √5)
   t = sqrt(√5 + 5)
   x2 = r1*((15√2 + 7√10)t + (11√10 + 25√2)s)/(22√5 + 50) - 2sqrt(-(850 + 380√5)r1^2 + (190 + 85√5)r1^2*s*t + (152√5 + 340)r2^2)/(11√5 + 25)
   y2 = 2R*y18 - R
   @printf("2r1 = %g;  2R = %g;  2r3 = %g\n", 2r1, 2R, 2r3)
   x = R .* [x18, 0, -x18, -x54, x54, x18]
   y = R .* [y18, 1, y18, -y54, -y54, y18]
   plot(x, y, color=:blue, lw=0.5)
   segment(-R*x18, R*y18, R*x18, R*y18, :green)
   segment(-R*x54, -R*y54, R*x18, R*y18, :green)
   segment(R*x54, -R*y54, -R*x18, R*y18, :green)
   circle(0, R*y18 + r1, r1)
   circle(0, R*y18 - r1, r1)
   circle(0, r3 - R*y54, r3, :brown)
   circle2(x2, y2, r2, :magenta)
   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(R*x18, R*y18, " A", :black, :left, :vcenter)
       point(-R*x18, R*y18, "B ", :black, :right, :vcenter)
       point(0, 2R*y18 - R, "C  ", :black, :right, :vcenter)
       point(-R*x54, -R*y54, "D ", :black, :right, :vcenter)
       point(R*x54, -R*y54, " E", :black, :left, :vcenter)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R*x18, R*y18, "(R*x18,R*y18)  ", :blue, :right, :bottom, delta=delta/2)
       point(R*x54, -R*y54, " (R*x54,-R*y54)   ", :blue, :right, :bottom, delta=delta/2)
       point(0, R*y18 + r1, "甲円:r1\n(0,R*y18+r1)", :red, :center, delta=-delta/2)
       point(0, R*y18 - r1, "甲円:r1\n(0,R*y18-r1)", :red, :center, delta=-delta/2)
       point(x2, 2R*y18 - R, "乙円:r2\n(x2,2R*y18-R)", :magenta, :center, delta=-delta/2)
       point(0, r3 - R*y54, "丙円:r3\n(0,r3-2R*y54)", :brown, :center, :bottom, delta=delta/2)
       plot!(xlims=(-1.1R*x18, 1.1R*x18))
   end
end;