算額（その442）

算額（その442）

武州豐嶋郡王子稲荷社 寛政元年己酉3月(1624)
關流神谷定令門人 東都湯嶋天神御徒町 大川三之助久菶

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

長方形内に甲円 2 個，乙円 1 個，丙円 2 個があり，長方形に内接し，互いに外接している。
長方形の短辺が 451 寸のとき，長辺の長さはいかほどか。

長方形の下辺の中点を原点とし，長方形の短辺，長辺をそれぞれ a, 2b とおく。
図形の右半分を対象とする。
甲円の半径，中心座標を r1, (b - r1, a - r1) ただし b = 2r1
乙円の半径，中心座標を r2, (0, r2)
丙円の半径，中心座標を r3, (b - r3, r3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive;
@syms a, b, r1, r2, r3
b = 2r1
eq1 = (b - r1)^2 + (a - r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (a - r1 - r3)^2 - (r1 + r3)^2
eq3 = (b - r3)^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))

   6-element Vector{Tuple{Sym, Sym, Sym}}:
    (a*(3 - 2*sqrt(2)), 2*a*(3 - 2*sqrt(2)), 2*a)
    (a*(2*sqrt(2) + 3), 2*a*(2*sqrt(2) + 3), 2*a)
    (a*(14*sqrt(2) + 2*sqrt(82 + 80*sqrt(2)) + 105 + 13*sqrt(41 + 40*sqrt(2)))/(14*(-3 + 4*sqrt(2) + sqrt(41 + 40*sqrt(2)))), a*(-301 - 17*sqrt(41 + 40*sqrt(2)) + 20*sqrt(82 + 80*sqrt(2)) + 224*sqrt(2))/(56*(-3 + 4*sqrt(2) + sqrt(41 + 40*sqrt(2)))), a*(5/4 + sqrt(2) + sqrt(41/4 + 10*sqrt(2))/2))
    (a*(-105 + 14*sqrt(2) - 13*sqrt(41 - 40*sqrt(2)) + 2*sqrt(82 - 80*sqrt(2)))/(14*(3 + 4*sqrt(2) - sqrt(41 - 40*sqrt(2)))), a*(301 + 224*sqrt(2) + 17*sqrt(41 - 40*sqrt(2)) + 20*sqrt(82 - 80*sqrt(2)))/(56*(3 + 4*sqrt(2) - sqrt(41 - 40*sqrt(2)))), a*(-sqrt(2) + 5/4 + sqrt(41/4 - 10*sqrt(2))/2))
    (a*(-105 + 14*sqrt(2) - 2*sqrt(82 - 80*sqrt(2)) + 13*sqrt(41 - 40*sqrt(2)))/(14*(3 + 4*sqrt(2) + sqrt(41 - 40*sqrt(2)))), a*(301 + 224*sqrt(2) - 20*sqrt(82 - 80*sqrt(2)) - 17*sqrt(41 - 40*sqrt(2)))/(56*(3 + 4*sqrt(2) + sqrt(41 - 40*sqrt(2)))), a*(-sqrt(2) + 5/4 - sqrt(41/4 - 10*sqrt(2))/2))
    (a*(-105 - 14*sqrt(2) + 2*sqrt(82 + 80*sqrt(2)) + 13*sqrt(41 + 40*sqrt(2)))/(14*(-4*sqrt(2) + 3 + sqrt(41 + 40*sqrt(2)))), a*(-224*sqrt(2) - 17*sqrt(41 + 40*sqrt(2)) + 20*sqrt(82 + 80*sqrt(2)) + 301)/(56*(-4*sqrt(2) + 3 + sqrt(41 + 40*sqrt(2)))), a*(-sqrt(41/4 + 10*sqrt(2))/2 + 5/4 + sqrt(2)))

r1 > r2 > r3 なので，6 番目の組が適解である。
r1, r2, r3 は以下のように簡約化される。

res[6][1] |> simplify |> apart |> simplify |> println
res[6][2] |> simplify |> apart |> simplify |> println
res[6][3] |> simplify |> apart |> simplify |> println

   a*(-2*sqrt(82 + 80*sqrt(2)) + 7 + sqrt(41 + 40*sqrt(2)) + 14*sqrt(2))/28
   a*(-11*sqrt(41 + 40*sqrt(2)) - 28*sqrt(2) + 63 + 8*sqrt(82 + 80*sqrt(2)))/112
   a*(-sqrt(41 + 40*sqrt(2)) + 5 + 4*sqrt(2))/4

長方形の長辺の長さは 4r1 = 4a*(-2*sqrt(82 + 80*sqrt(2)) + 7 + sqrt(41 + 40*sqrt(2)) + 14*sqrt(2))/28

短辺の長さ = a = 451 寸のとき,長辺の長さは約 563 寸である。

a = 451
4a*(-2*sqrt(82 + 80*sqrt(2)) + 7 + sqrt(41 + 40*sqrt(2)) + 14*sqrt(2))/28

   563.0009311237039

ついでに，r1, r2, r3 はそれぞれ 140.75023278092598, 106.71276946728344, 87.85200899283763 である。

a*(-2*sqrt(82 + 80*sqrt(2)) + 7 + sqrt(41 + 40*sqrt(2)) + 14*sqrt(2))/28 |> println
a*(-11*sqrt(41 + 40*sqrt(2)) - 28*sqrt(2) + 63 + 8*sqrt(82 + 80*sqrt(2)))/112 |> println
a*(-sqrt(41 + 40*sqrt(2)) + 5 + 4*sqrt(2))/4 |> println

   140.75023278092598
   106.71276946728344
   87.85200899283763

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 451
   temp = sqrt(41 + 40*sqrt(2))
   r1 = a*((1 - 2*√2)temp + 7 + 14*sqrt(2))/28
   r2 = a*((8*√2 - 11)temp - 28*sqrt(2) + 63)/112
   r3 = a*(-temp + 5 + 4*sqrt(2))/4
   b = 2r1
   @printf("a = %g;  r1 = %g;  r2 = %g;  r3 = %g\n", a, r1, r2, r3)
   @printf("長方形の長辺の長さ = %g\n", 4r1)
   plot([b, b, -b, -b, 2r1], [0, a, a, 0, 0], color=:black, lw=0.5)
   circle(r1, a - r1, r1, :green)
   circle(-r1, a - r1, r1, :green)
   circle(0, r2, r2, :magenta)
   circle(b - r3, r3, r3, :blue)
   circle(-b + r3, r3, r3, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(b - r1, a - r1, "甲円:r1,(b-r1,a-r1)", :green, :center, :top, delta=-delta)
       point(0, r2, "乙円:r2,(0,r2) ", :magenta, :center, :top, delta=-delta)
       point(b - r3, r3, "丙円:r3,(b-r3,r3) ", :blue, :center, :top, delta=-delta)
       point(b, a, "(b,a)", :black, :right, :bottom, delta=delta)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;