算額（その425）

算額（その425）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

長方形内に，甲，乙，丙，丁，戊，己の 6 個の円がある。それぞれ長方形に内接したり，互いに外接しあっている。
己円の直径を丁円と戊円の直径で表わせ。

長方形の長辺の長さを a とする（短辺の長さは 甲円の直径に等しい）。
甲円の半径と中心座標 r1, (a - r1, r1)
乙円の半径と中心座標 r2, (r2, r2)
丙円の半径と中心座標 r3, (x3, 2r1 - r3)
丁円の半径と中心座標 r4, (a - r4, r4)
戊円の半径と中心座標 r5, (r5, r5)
己円の半径と中心座標 r6, (x6, y6)

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
     x3::positive, r4::positive, r5::positive, r6::positive,
     x6::positive, y6::positive;

eq1 = (a - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (2r1 - r3 - r1)^2 - (r1 + r3)^2
eq3 = 2(r1 - r4)^2 - (r1 + r4)^2
eq4 = (a - r1 - x6)^2 + (r1 - y6)^2 - (r1 + r6)^2
eq5 = (x3 - r2)^2 + (2r1 - r3 - r2)^2 - (r2 + r3)^2
eq6 = 2(r2 - r5)^2 - (r2 + r5)^2
eq7 = (r2 - x6)^2 + (r2 - y6)^2 - (r2 + r6)^2
eq8 = (x3 - x6)^2 + (2r1 - r3 - y6)^2 - (r3 + r6)^2
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (a, r1, r2, r3, x3, r6, x6, y6))

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, r1, r2, r3, x3, r6, x6, y6) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (a - r1 - r2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - r3)^2,  # eq2
       2*(r1 - r4)^2 - (r1 + r4)^2,  # eq3
       -(r1 + r6)^2 + (r1 - y6)^2 + (a - r1 - x6)^2,  # eq4
       (-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2,  # eq5
       2*(r2 - r5)^2 - (r2 + r5)^2,  # eq6
       -(r2 + r6)^2 + (r2 - x6)^2 + (r2 - y6)^2,  # eq7
       -(r3 + r6)^2 + (x3 - x6)^2 + (2*r1 - r3 - y6)^2,  # eq8
   ]
end;

(r4, r5) = (2, 1)
iniv = [big"180.0", 59, 35, 21, 55, 5, 65, 62] ./2
res = nls(H, ini=iniv);
println(res);

  (BigFloat[33.97056274847714058562026469051637694283606250452337687812015685588878974154529, 11.65685424949238019520675489683879231427868750150779229270671895196292991384848, 5.828427124746190097603377448419396157139343750753896146353359475981464956924239, 5.828427124746190097603377448419396157139343750753896146353359475981464956924308, 5.82842712474619009760948182491957344489198207867602126024818280529793350621324, 1.093836321356054313599895375610329707159527028223527974012807033844387485155129, 9.563017928136325881606859521228462607119160473284264318693911918118542430249862, 11.65685424949238019520479919353753625575412910164858126097863477985370829598065], true)

r4 = 2;  r5 = 1
a = 33.9706;  r1 = 11.6569;  r2 = 5.82843;  r3 = 5.82843;  x3 = 5.82843;  r6 = 1.09384;  x6 = 9.56302;  y6 = 11.6569

己円の直径は 2.18767 であるが，術では以下のように 7.31370849898476 であり，答えが違う。中村が術を読み違えているのか，式に起こすときに誤記があるのか。

(丁, 戊) = (2, 1) .* 2
天 = 丁 * 戊
地 = 丁 + 戊
法 = 2(sqrt((4地 * 戊 + 丁^2)天) + 2戊^2) + 地*丁
法*丁/(3 + √8)天

   7.31370849898476

中村は下記の x1 式を簡約化して x2 を提示し術と同じとしている。
しかし，x1 は正しいが x2 は式の記述ミスがある。正しくは x3（分母の式の一箇所の - が + でなければならない）。

@syms d, e
x1 = (3 + 2sqrt(Sym(2)))*d^2*e*((d^2 + d*e + 4e^2) - 2(d + 2e)*sqrt(d*e))/(d^4 - 2d^3*e - 7d^2*e^2 - 8d*e^3 + 16e^4)
x2 = (3 + 2sqrt(Sym(2)))*d^2*e / (d^2 + d*e + 4e^2 - 2(d + 2e)*sqrt(d*e))
x3 = (3 + 2sqrt(Sym(2)))*d^2*e / (d^2 + d*e + 4e^2 + 2(d + 2e)*sqrt(d*e))
x1(d => 4, e => 2).evalf() |> println
x2(d => 4, e => 2).evalf() |> println
x3(d => 4, e => 2).evalf() |> println

   2.18767264271211
   -35.4929704984046
   2.18767264271211

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, r5) = (2, 1)
   (a, r1, r2, r3, x3, r6, x6, y6) = res[1]
   @printf("r4 = %g;  r5 = %g\n", r4, r5)
   @printf("a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  r6 = %g;  x6 = %g;  y6 = %g\n", a, r1, r2, r3, x3, r6, x6, y6)
   @printf("己円の直径 = %g\n", 2r6)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(x3, 2r1 - r3, r3, :green)
   circle(a - r4, r4, r4, :magenta)
   circle(r5, r5, r5, :orange)
   circle(x6, y6, r6, :purple)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a - r1, r1, " 甲円:r1,(a-r1,r1)", :red, :left, :vcenter)
       point(r2, r2, " 乙円:r2,(r2,r2)", :blue, :center, :top, delta=-delta)
       point(x3, 2r1 - r3, " 丙円:r3,(x3,2r1-r3)", :green, :center, :top, delta=-delta)
       point(a - r4, r4, "丁円:r4,(a-r4,r4) ", :magenta, :right, :vcenter)
       point(r5, r5, "    戊円:r5,(r5,r5) ", :orange, :left, :vcenter)
       point(x6, y6, "  己円:r6,(x5,y5) ", :purple, :left, :vcenter)
       point(a, 2r1, "(a,2r1) ", :black, :right, :top, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;