算額（その667）

算額（その667）

長野市元善町 善光寺 寛政8年(1796)
中村信弥「改訂増補 長野県の算額」(p.56)
http://www.wasan.jp/zoho/zoho.html

台形内に甲円が内接する直角三角形を入れる。乙円，丙円も台形の辺と直角三角形の辺に内接する。甲円，乙円，丙円の直径がそれぞれ 20寸，24寸，10寸のとき，甲円が内接する直角三角形の直覚を挟む二辺のうち短い方の辺（鈎）の長さはいかほどか。

甲円が内接する直角三角形の頂点座標を (b, 0), (0, c), (a, d) とする。
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (a - r3, r3)
とおき，以下の連立方程式を解く。
求めるべき，「甲円が内接する直角三角形の直覚を挟む二辺のうち短い方の辺（鈎）の長さ」は (b, 0), (a, d) を端点とする直線の長さ sqrt((a - b)^2 + d^2) である。

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

using SymPy
@syms a::positive, b::positive, c::positive, d::positive,
     x1::positive, y1::positive,
     r1::positive, r2::positive, r3::positive,
     鈎::positive, 股::positive, 弦::positive
弦 = sqrt(a^2 + (c - d)^2)
鈎 = sqrt((a - b)^2 + d^2)
股 = sqrt(b^2 + c^2)
eq1 = 鈎^2 + 股^2 - 弦^2 |> simplify
eq2 = 鈎 + 股 - 弦 - 2r1
eq3 = b + c - 股 - 2r2
eq4 = d + (a - b) - 鈎 -2r3
eq5 = dist(0, c, a, d, x1, y1) - r1^2
eq6 = dist(b, 0, a, d, x1, y1) - r1^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (a, b, c, d, x1, y1))

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, b, c, d, x1, y1) = u
   return [
       -2*a*b + 2*b^2 + 2*c*d,  # eq1
       -2*r1 - sqrt(a^2 + (c - d)^2) + sqrt(b^2 + c^2) + sqrt(d^2 + (a - b)^2),  # eq2
       b + c - 2*r2 - sqrt(b^2 + c^2),  # eq3
       a - b + d - 2*r3 - sqrt(d^2 + (a - b)^2),  # eq4
       -r1^2 + (-a*(a*x1 + (-c + d)*(-c + y1))/(a^2 + (-c + d)^2) + x1)^2 + (-c + y1 - (-c + d)*(a*x1 + (-c + d)*(-c + y1))/(a^2 + (-c + d)^2))^2,  # eq5
       -r1^2 + (-d*(d*y1 + (a - b)*(-b + x1))/(d^2 + (a - b)^2) + y1)^2 + (-b + x1 - (a - b)*(d*y1 + (a - b)*(-b + x1))/(d^2 + (a - b)^2))^2,  # eq6
   ]
end;

(r1, r2, r3) = (20, 24, 10) .// 2
iniv = BigFloat[59, 46, 37, 22, 43, 15]
res = nls(H, ini=iniv)

   (BigFloat[63.00000000000000000000000000000000000000000000000000000000000000635922298913147, 48.00000000000000000000000000000000000000000000000000000000000001088653549950012, 35.99999999999999999999999999999999999999999999999999999999999998913806908318885, 20.00000000000000000000000000000000000000000000000000000000000000452357741848939, 46.00000000000000000000000000000000000000000000000000000000000000725087619460416, 13.99999999999999999999999999999999999999999999999999999999999999999791841656566], true)

「甲円が内接する直角三角形の直覚を挟む二辺のうち短い方の辺（鈎）の長さ」は (b, 0), (a, d) を端点とする線分の長さ sqrt((a - b)^2 + d^2) である。
a = 63, b = 48, d = 20 なので，線分の長さは 25 である。

その他のパラメータは以下の通りである。

   a = 63;  b = 48;  c = 36;  d = 20;  x1 = 46;  y1 = 14
   鈎 = 25;  股 = 60;  弦 = 65

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, d, x1, y1) = res[1]
   @printf("a = %g;  b = %g;  c = %g;  d = %g;  x1 = %g;  y1 = %g\n", a, b, c, d, x1, y1)
   鈎 = sqrt((a - b)^2 + d^2)
   股 = sqrt(b^2 + c^2)
   弦 = sqrt(a^2 + (c - d)^2)
   @printf("鈎 = %g;  股 = %g;  弦 = %g\n", 鈎, 股, 弦)
   plot([0, a, a, 0, 0], [0, 0, d, c, 0], color=:black, lw=0.5)
   plot!([0, b, a], [c, 0, d], color=:orange, lw=0.5)
   circle(x1, y1, r1)
   circle(r2, r2, r2, :blue)
   circle(a - r3, r3, r3, :green)
   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(x1, y1, "甲円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(r2, r2, "乙円:r2,(r2,r2)", :blue, :center, delta=-delta/2)
       point(a - r3, r3, "丙円:r3,(a-r3,r3)", :black, :right, delta=-2delta)
       point(a, 0, " a", :black, :center, delta=-1.5delta)
       point(b, 0, " b", :black, :center, delta=-1.5delta)
       point(0, c, "c ", :black, :right, :vcenter)
       point(a, d, "(a,d)", :black, :center, :bottom, delta=2delta)
       plot!(xlims=(-3, 66), ylims=(-3, 37))
   end
end;