算額（その1392）

算額（その1392）

十九 羽生市中手子林 八幡神社 文政元年(1818)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，三角形，正六角形
#Julia, #SymPy, #算額, #和算

（不等辺）三角形の中に，天円，地円，人円と正六角形を容れる。天円，地円の直径が 4 寸，6 寸のとき，人円の直径はいかほどか。

三角形の頂点座標を (3a, 0), (x0, y0), (b, 0)
正六角形の一辺の長さと中心座標を 2a, (0, √3a)
天円の半径と中心座標を r1, (x1, 2√3a + r1)
地円の半径と中心座標を r2, (2a, r2)
人円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms a::positive, b::negative,
     r1::positive, x1::positive,
     r2::positive,
     r3::positive, x3::negative
# a = √Sym(3)r2

eq1 = dist(-a, 0, -2a, √Sym(3)a, x3, r3) - r3^2
eq2 = dist(b, 0, -a, 2√Sym(3)a, x3, r3) - r3^2
eq3 = dist(b, 0, -a, 2√Sym(3)a, x1, 2√Sym(3)a + r1) - r1^2
eq4 = dist(3a, 0, 2a, √Sym(3)a, x1, 2√Sym(3)a + r1) - r1^2;

一度には解けないので，eq1, eq3, eq4 を連立方程式として，b, x1, x3 を求める。
なお，a = √3r2 であるが，a はそのままにしておく。

res = solve([eq1, eq3, eq4], (b, x1, x3))[3]  # 3 of 4

   (a*(24*sqrt(3)*a^3 - 48*a^2*r1 - 4*sqrt(3)*a*r1^2 + 3*r1^3)/(r1*(-12*a^2 + 4*sqrt(3)*a*r1 + 3*r1^2)), a - sqrt(3)*r1, -a - sqrt(3)*r3)

b, x1, x3 を eq2 に代入し(eq12 とする），r3 を求める。

eq12 = eq2(b => res[1], x1 => res[2], x3 => res[3]) |> simplify;
res2 = solve(eq12, r3)[1];  # 1 of 2

ここで，a = √Sym(3)r2 を代入し，因数分解する。
結果は，長精度整数を含む長い式になる。

ans_r3 = res2(a => √Sym(3)r2) |> factor;
ans_r3 |> println

   -(r1^16 + 20*r1^15*r2 + 70*r1^14*r2^2 - 876*r1^13*r2^3 - 5148*r1^12*r2^4 + 16974*r1^11*r2^5 + 110601*r1^10*r2^6 - 293274*r1^9*r2^7 - 1113426*r1^8*r2^8 + 4030560*r1^7*r2^9 + 2711232*r1^6*r2^10 - 27387072*r1^5*r2^11 + 34572096*r1^4*r2^12 + 22394880*r1^3*r2^13 - 95738112*r1^2*r2^14 + 90699264*r1*r2^15 - 30233088*r2^16 + sqrt(r1^32 + 34*r1^31*r2 + 321*r1^30*r2^2 - 1406*r1^29*r2^3 - 37418*r1^28*r2^4 - 54378*r1^27*r2^5 + 1866549*r1^26*r2^6 + 5556726*r1^25*r2^7 - 62431686*r1^24*r2^8 - 193177494*r1^23*r2^9 + 1667661669*r1^22*r2^10 + 3566941110*r1^21*r2^11 - 36334402551*r1^20*r2^12 - 22380759756*r1^19*r2^13 + 596476955268*r1^18*r2^14 - 616148677152*r1^17*r2^15 - 6288496227900*r1^16*r2^16 + 18622534127616*r1^15*r2^17 + 24103256904576*r1^14*r2^18 - 214678876124160*r1^13*r2^19 + 294385908328704*r1^12*r2^20 + 800063891564544*r1^11*r2^21 - 3753687721439232*r1^10*r2^22 + 5063157982642176*r1^9*r2^23 + 4824832605788160*r1^8*r2^24 - 33676529335271424*r1^7*r2^25 + 73507309933658112*r1^6*r2^26 - 99895501823016960*r1^5*r2^27 + 93672133367169024*r1^4*r2^28 - 61376067146612736*r1^3*r2^29 + 27116508430467072*r1^2*r2^30 - 7312316880125952*r1*r2^31 + 914039610015744*r2^32))/(r1*(r1 - 2*r2)^3*(r1 + 6*r2)^6*(2*r1 - 3*r2)*(r1^2 - 3*r1*r2 + 3*r2^2)^2)

式を簡約化しなくても，r1, r2 に数値を代入すれば，数値解が得られる。

ans_r3(r1 => 4//2, r2 => 6//2) |> println

   11/2

天円，地円の直径が 6 寸，4 寸のとき，人円の直径は 11 寸になる。

簡約化

SymPy では結果の式を自動では簡約化できないので，手作業を交えて簡約化を試みる。
r3 = (項1 - sqrt(項2)) / 項3 の形をしている。

まず，分子の (項1 - sqrt(項2)) を取り出す

num = numerator(ans_r3);

   -r1^16 - 20*r1^15*r2 - 70*r1^14*r2^2 + 876*r1^13*r2^3 + 5148*r1^12*r2^4 - 16974*r1^11*r2^5 - 110601*r1^10*r2^6 + 293274*r1^9*r2^7 + 1113426*r1^8*r2^8 - 4030560*r1^7*r2^9 - 2711232*r1^6*r2^10 + 27387072*r1^5*r2^11 - 34572096*r1^4*r2^12 - 22394880*r1^3*r2^13 + 95738112*r1^2*r2^14 - 90699264*r1*r2^15 + 30233088*r2^16 - sqrt(r1^32 + 34*r1^31*r2 + 321*r1^30*r2^2 - 1406*r1^29*r2^3 - 37418*r1^28*r2^4 - 54378*r1^27*r2^5 + 1866549*r1^26*r2^6 + 5556726*r1^25*r2^7 - 62431686*r1^24*r2^8 - 193177494*r1^23*r2^9 + 1667661669*r1^22*r2^10 + 3566941110*r1^21*r2^11 - 36334402551*r1^20*r2^12 - 22380759756*r1^19*r2^13 + 596476955268*r1^18*r2^14 - 616148677152*r1^17*r2^15 - 6288496227900*r1^16*r2^16 + 18622534127616*r1^15*r2^17 + 24103256904576*r1^14*r2^18 - 214678876124160*r1^13*r2^19 + 294385908328704*r1^12*r2^20 + 800063891564544*r1^11*r2^21 - 3753687721439232*r1^10*r2^22 + 5063157982642176*r1^9*r2^23 + 4824832605788160*r1^8*r2^24 - 33676529335271424*r1^7*r2^25 + 73507309933658112*r1^6*r2^26 - 99895501823016960*r1^5*r2^27 + 93672133367169024*r1^4*r2^28 - 61376067146612736*r1^3*r2^29 + 27116508430467072*r1^2*r2^30 - 7312316880125952*r1*r2^31 + 914039610015744*r2^32)

項1，項2 を因数分解した結果は，かなり短い式になる。

term1 = -r1^16 - 20*r1^15*r2 - 70*r1^14*r2^2 + 876*r1^13*r2^3 + 5148*r1^12*r2^4 - 16974*r1^11*r2^5 - 110601*r1^10*r2^6 + 293274*r1^9*r2^7 + 1113426*r1^8*r2^8 - 4030560*r1^7*r2^9 - 2711232*r1^6*r2^10 + 27387072*r1^5*r2^11 - 34572096*r1^4*r2^12 - 22394880*r1^3*r2^13 + 95738112*r1^2*r2^14 - 90699264*r1*r2^15 + 30233088*r2^16;
term1 = term1 |> factor
term1 |> println

   -(r1 - 2*r2)^2*(r1 + 6*r2)^6*(r1^2 - 3*r2^2)*(r1^2 - 6*r1*r2 + 6*r2^2)*(r1^2 - 3*r1*r2 + 3*r2^2)^2

term2 = r1^32 + 34*r1^31*r2 + 321*r1^30*r2^2 - 1406*r1^29*r2^3 - 37418*r1^28*r2^4 - 54378*r1^27*r2^5 + 1866549*r1^26*r2^6 + 5556726*r1^25*r2^7 - 62431686*r1^24*r2^8 - 193177494*r1^23*r2^9 + 1667661669*r1^22*r2^10 + 3566941110*r1^21*r2^11 - 36334402551*r1^20*r2^12 - 22380759756*r1^19*r2^13 + 596476955268*r1^18*r2^14 - 616148677152*r1^17*r2^15 - 6288496227900*r1^16*r2^16 + 18622534127616*r1^15*r2^17 + 24103256904576*r1^14*r2^18 - 214678876124160*r1^13*r2^19 + 294385908328704*r1^12*r2^20 + 800063891564544*r1^11*r2^21 - 3753687721439232*r1^10*r2^22 + 5063157982642176*r1^9*r2^23 + 4824832605788160*r1^8*r2^24 - 33676529335271424*r1^7*r2^25 + 73507309933658112*r1^6*r2^26 - 99895501823016960*r1^5*r2^27 + 93672133367169024*r1^4*r2^28 - 61376067146612736*r1^3*r2^29 + 27116508430467072*r1^2*r2^30 - 7312316880125952*r1*r2^31 + 914039610015744*r2^32;
term2 = term2 |> factor
term2 |> println

   (r1 - 2*r2)^4*(r1 + 6*r2)^12*(r1^2 - 6*r1*r2 + 6*r2^2)^2*(r1^2 - 3*r1*r2 + 3*r2^2)^6

分子を書き直し，因数分解する。

num = term1 -(r1 - 2*r2)^2*(r1 + 6*r2)^6*(r1^2 - 6*r1*r2 + 6*r2^2)*(r1^2 - 3*r1*r2 + 3*r2^2)^3 |> factor;
num |> println

   -r1*(r1 - 2*r2)^2*(r1 + 6*r2)^6*(2*r1 - 3*r2)*(r1^2 - 6*r1*r2 + 6*r2^2)*(r1^2 - 3*r1*r2 + 3*r2^2)^2

分子/分母 を計算し，簡約化する。

term3 = denominator(ans_r3);

ans__r3 = num/term3 |> simplify
ans__r3 |> println

   (-r1^2 + 6*r1*r2 - 6*r2^2)/(r1 - 2*r2)

得られた式 (-r1^2 + 6*r1*r2 - 6*r2^2)/(r1 - 2*r2) は術と同じである。

ans__r3(r1 => 4//2, r2 => 6//2) |> println

   11/2

術は，
人円径 = 地円径*(6地円径 - 4天円径)/(2地円径 - 天円径) - 天円径
であるが，簡約化すると
人円径 = (6*地円径^2 - 6*地円径*天円径 + 天円径^2)/(2*地円径 - 天円径)
になる。

@syms 天円径, 地円径
人円径 = 地円径*(6地円径 - 4天円径)/(2地円径 - 天円径) - 天円径;
人円径 |> simplify |> println
人円径(天円径 => 4, 地円径 => 6) |> println

   (6*地円径^2 - 6*地円径*天円径 + 天円径^2)/(2*地円径 - 天円径)
   11

# (ox, oy) を中心，r を半径とする円に内接する正 n 多角形。fcol を指定すると塗りつぶし。
function polygon(ox, oy, r, n; φ=90, color=:black, lt=:solid, lw=0.5, fillcolor="")
   θ = range(0, 2π, length=n+1) .+ pi*φ/180
   if fillcolor == ""
       plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
             linecolor=color, linestyle=lt, linewidth=lw)
   else
       plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
             linecolor=color, linestyle=lt, linewidth=lw,
             seriestype=:shape, fillcolor=fillcolor)
   end
end;

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = √3r2
   r3 = (r1^2 - 6r1*r2 + 6r2^2)/(2r2 - r1)
   (b, x1, x3) = (a*(24*sqrt(3)*a^3 - 48*a^2*r1 - 4*sqrt(3)*a*r1^2 + 3*r1^3)/(r1*(-12*a^2 + 4*sqrt(3)*a*r1 + 3*r1^2)), a - sqrt(3)*r1, -a - sqrt(3)*r3)
   (x0, y0) = float.(intersection(b, 0, -a, 2√3a, 3a, 0, 2a, √3a))
   @printf("天円，地円の直径が %g, %g のとき，人円の直径は %g である。\n", 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  a = %g;  b = %g;  x1 = %g;  x3 = %g\n", r1, r2, r3, a, b, x1, x3)
   plot([b, x0, 3a, b], [0, y0, 0, 0], color=:green, lw=0.5)
   polygon(0, √3a, 2a, 6, φ=60, color=:magenta)
   circle(x1, 2√3a + r1, r1, :orange)
   circle(2a, r2, r2)
   circle(x3, r3, r3, :blue)
   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(a, 0, "a", :green, :center, delta=-2delta)
       point(2a, 0, "2a", :green, :center, delta=-2delta)
       point(3a, 0, "3a", :green, :center, delta=-2delta)
       point(0, 2√3a, "2√3a", :magenta, :center, delta=-2delta)
       point(x1, 2√3a + r1, " 天円:r2,(x1,2√3a+r1)", :black, :left, :vcenter)
       point(2a, r2, " 地円:r2,(2a,r2)", :red, :left, :vcenter)
       point(x3, r3, "人円:r3\n(x3,r3)", :red, :center, delta=-delta)
       point(x0, y0, " (x0,y0)", :green, :left, :vcenter)
       point(b, 0, "b", :green, :center, delta=-2delta)
       plot!(ylims=(-8delta, y0 + 3delta), xlims=(b - 3delta, 3a + 25delta))
   end
end;

draw(4/2, 6/2, true)