算額（その2030）

算額（その2030）

福島県小野町 東堂山満福寺観音堂 天保8年(1837)
https://www.nippon.com/ja/japan-topics/c12803/

キーワード：円4個，直角三角形

外円の中に，直角三角形と，その 3 辺の中点を中心とする 3 個の円を容れる。鈎，股が 3 寸，4 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (x, y)
鈎，股，弦を a, b, c
それぞれの中点を中心とする円の半径を r1, r2, r3
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, R::positive, r1::positive, r2::positive, r3::positive, x::positive, y::positive
@syms a, b, c, R, r1, r2, r3, x, y
#c = sqrt(a^2 + b^2)
r1 = a/2
r2 = b/2
r3 = c/2
eq1 = x^2 + (r1 - y)^2 - (R - r1)^2
eq2 = (r2 - x)^2 + y^2 - (R - r2)^2
eq3 = (r2 - x)^2 + (r1 - y)^2 - (R - r3)^2
res = solve([eq1, eq2, eq3], (R, x, y))[1]  # 1 of 2

   ((a^5 - a^4*c - a^3*b^2 - a^3*c^2 - a^2*b^3 + a^2*c^3 - sqrt(2)*a^2*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - sqrt(2)*a*b*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + sqrt(2)*a*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + b^5 - b^4*c - b^3*c^2 + b^2*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)), (a^4*b^2 + a^3*b^3 - 3*a^3*b^2*c + sqrt(2)*a^3*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - 3*a^2*b^3*c + 3*a^2*b^2*c^2 + sqrt(2)*a^2*b*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - 2*sqrt(2)*a^2*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - a*b^5 + a*b^4*c + a*b^3*c^2 - a*b^2*c^3 - sqrt(2)*a*b*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + sqrt(2)*a*c^2*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + b^6 - b^5*c - b^4*c^2 + b^3*c^3)/(4*b*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)), a*(a^4 - a^3*b - a^3*c + a^2*b*c - a^2*c^2 + a*b^3 - 3*a*b^2*c + a*b*c^2 + a*c^3 + b^4 - 3*b^3*c + 3*b^2*c^2 - b*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)) + sqrt(2)*sqrt(-a*b^3*(a - b - c)*(a - b + c))*(b - c)*(a + b - c)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)))

R, x, y を a, b, c を含む式で表すと，SymPy では簡約化できない複雑な長い式（省略）になる。
鈎，股として（弦は別途計算），具体的な数値を代入すればよい。

鈎，股が 3 寸，4 寸のとき，外円の直径は 2*72/23 = 6 + 6/23 である。

res[1](a => 3, b => 4, c => 5) |> println 
res[2](a => 3, b => 4, c => 5) |> println 
res[3](a => 3, b => 4, c => 5) |> println 

   72/23
   36/23
   24/23

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = sqrt(a^2 + b^2)
   (R, x, y) = ((a^5 - a^4*c - a^3*b^2 - a^3*c^2 - a^2*b^3 + a^2*c^3 - sqrt(2)*a^2*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - sqrt(2)*a*b*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + sqrt(2)*a*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + b^5 - b^4*c - b^3*c^2 + b^2*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)), (a^4*b^2 + a^3*b^3 - 3*a^3*b^2*c + sqrt(2)*a^3*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - 3*a^2*b^3*c + 3*a^2*b^2*c^2 + sqrt(2)*a^2*b*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - 2*sqrt(2)*a^2*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - a*b^5 + a*b^4*c + a*b^3*c^2 - a*b^2*c^3 - sqrt(2)*a*b*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + sqrt(2)*a*c^2*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + b^6 - b^5*c - b^4*c^2 + b^3*c^3)/(4*b*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)), a*(a^4 - a^3*b - a^3*c + a^2*b*c - a^2*c^2 + a*b^3 - 3*a*b^2*c + a*b*c^2 + a*c^3 + b^4 - 3*b^3*c + 3*b^2*c^2 - b*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)) + sqrt(2)*sqrt(-a*b^3*(a - b - c)*(a - b + c))*(b - c)*(a + b - c)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)))
   @printf("a = %g;  b = %g;  c = %g;  R = %g;  x = %g;  y = %g\n", a, b, c, R, x, y)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:green, lw=0.5)
   circle(x, y, R)
   circle(0, a/2, a/2, :blue)
   circle(b/2, 0, b/2, :magenta)
   circle(b/2, a/2, c/2, :orange)
   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(x, y, "外円:R,(x,y)", :red, :center, delta=-delta/2)
       point(0, a/2, "小円:r1,(0,a/2)", :blue, :center, delta=-delta/2)
       point(b/2, 0, "中円:r2,(b/2,0)", :magenta, :center, delta=-delta/2)
       point(b/2, a/2, "大円:r3,(b/2,a/2)", :orange, :center, delta=-delta/2)
       point(0, a, "a ", :green, :right, :bottom, delta=delta)
       point(b, 0, " b", :green, :left, delta=-delta)
   end
end;

draw(3, 4, true)