算額（その659）

算額（その659）

長野県上田市別所温泉 北向観音堂 慶応元年(1865)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

外円内に正三角形が内接し，天円，地円，人円が入っている。外円の直径の平方根が xxx.625，天円の直径の平方根が yyy.96875 のようになっている（小数部がわかっている）。条件を満たすような最小の外円の直径を求めよ。

基本的には r1 と r3 だけの関係式が分かれば問題は解けるのだが，それだけにこだわると逆に落とし穴にハマる。R がわかっているときに，残りすべてのパラメータを求めよう。
R がわかっているときに，r1, r2, r3, x2 の 4 個のパラメータを推定するには，条件が 4 個必要である。eq1 〜 eq3 は比較的簡単に導けるが，eq4 はちょっと盲点かもしれない。

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

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r2::positive, x2::positive, r3::positive,
     constant::positive

eq1 = x2^2 + (R - r1 - (r2 - R/2))^2 - (r1 + r2)^2
eq2 = x2^2 + (r2 - r3)^2 - (r2 + r3)^2
eq3 = 2r1 + 2r3 - 3R/2
eq4 = 2sqrt(r2*r3) + sqrt(Sym(3))*r2 - R*sqrt(Sym(3))/2

res = solve([eq1, eq2, eq3], (r2, x2, r3))
res = solve([eq1, eq2, eq3, eq4], (r1, r2, x2, r3))

   1-element Vector{NTuple{4, Sym}}:
    (9*R/16, R/4, sqrt(3)*R/4, 3*R/16)

つまり，天円の半径 r1 は外円の半径 R の 9/16 である。
解を求めるには以下の3条件を満たす整数 a, b を求めればよい。
いままで，r1, R は半径としてきたが，直径の場合は単に半径を2敗するだけであるから関係式は変わらない。
これ以降 r1, R は直径として読む。
r1/R = 9/16
sqrt(R) = a + 0.625
sqrt(r1) = b + 0.96875

ブルートフォースで該当する解を探索する。

function search()
   for a = 1:10
       R = (a + 0.625)^2
       for b = 1:a
           r1 = (b + 0.96875)^2
           r1 > R && break
           if isapprox(r1/R, 9/16)
               @printf("R = %g;  r1 = %g\n", R, r1)
               return
           end
       end
   end
end
search()    

   R = 6.89062;  r1 = 3.87598

条件を満たす外円の直径は 6.89062 である。
ちなみに，そのときの，各パラメータは以下の通り。
R = 6.89062;  r1 = 3.87597;  r2 = 1.72266;  r3 = 1.29199;  x2 = 2.98373

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (6.89062, 3.87598)
   (r1, r2, x2, r3) = R .* (9/16, 1/4, √3/4, 3/16)
   @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x2 = %g\n", R, r1, r2, r3, x2)
   plot([√3R/2, 0, -√3R/2, √3R/2], [-R/2, R, -R/2], color=:green, lw=0.5)
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x2, r2 - R/2, r2, :magenta)
   circle(-x2, r2 - R/2, r2, :magenta)
   circle(0, r3 - R/2, r3, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, R - r1, " 天円:r1\n (0,R-r1)", :black, :left, :vcenter)
       point(x2, r2 - R/2, " 地円:r2\n (√3R/2,-R/2)", :black, :left, :vcenter)
       point(0, r3 - R/2, " 人円:r3\n (0,r3-R/2)", :black, :left, :vcenter)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
   end
end;