算額（その99）

算額（その99）

東京都渋谷区渋谷 金王八幡神社 安政6年(1859)4月
http://www.wasan.jp/tokyo/konnou2.html

参照先の画像では問題文が全く読めない

下図のようjに記号を定め，方程式を解く。

しかし，解は複素数になる。もっとも，虚部はほとんど 0 なのだが，厳密に言えば「@syms R::positive」などとした場合には，「解なし」ということになる。

解の式は恐ろしく複雑になるが，数値にすると以下のようになる。最初の組が適切な解である。

using SymPy
@syms R::positive, r::positive, x::positive, y::positive;
@syms R, r, x, y;
eq1 = r^2 + y^2 - (r + 593)^2;
eq2 = r^2 + (y - R)^2 - (R - r)^2;
eq3 = (593 - r)^2 + R^2 - (R + r)^2;
res = solve([eq1, eq2, eq3], (R, r, y))

   3-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-1186*sqrt(3) - 2965*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/6 - 593*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/6 + 1186*I/3 + 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/6 + 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/6)/(3*sqrt(3) - I), (593*sqrt(3) - 593*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/3 - 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/6 - 593*I/3 + 593*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/12 + 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/12)/(3*sqrt(3) - I), -593/3 + 4151/(9*(7/54 + 7*sqrt(3)*I/18)^(1/3)) + 593*(7/54 + 7*sqrt(3)*I/18)^(1/3))
    ((-1186*sqrt(3) - 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/18 - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/18 - 1186*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/9 - 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/9 + 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 + 1186*I/3 + 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/18 + 2965*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/18)/(3*sqrt(3) - I), (593*sqrt(3) - 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18 - 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/36 - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/36 - 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 - 593*I/3 + 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/9 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/36 + 1186*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/9 + 2965*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/36)/(3*sqrt(3) - I), -593/3 + 593*(-1/2 - sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3) + 4151/(9*(-1/2 - sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3)))
    ((-1186*sqrt(3) - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/9 - 1186*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/9 - 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18 - 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/18 + 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 + 1186*I/3 + 2965*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/18 + 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/18 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/18)/(3*sqrt(3) - I), (593*sqrt(3) - 2965*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/18 - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/18 - 593*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/9 - 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 - 593*I/3 + 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/18 + 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/36 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/36 + 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18)/(3*sqrt(3) - I), -593/3 + 4151/(9*(-1/2 + sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3)) + 593*(-1/2 + sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3))

for j = 1:length(res)
   println("解 $j")
   println("R = $(res[j][1].evalf())")
   println("r = $(res[j][2].evalf())")
   println("y = $(res[j][3].evalf())")
   println()
end

   解 1
   R = 475.549077332269 - 6.93889390390723e-18*I
   r = 164.545086163906
   y = 739.458905004458 + 0.e-20*I

   解 2
   R = -1332.45890500446 + 1.38777878078145e-17*I
   r = -237.774538666135 + 3.46944695195361e-18*I
   y = -263.909827672189 + 0.e-17*I

   解 3
   R = -329.090172327811
   r = 666.229452502229 - 6.93889390390723e-18*I
   y = -1068.54907733227 + 0.e-18*I

解説したページ https://sites.google.com/site/seijiimura/home/09-san-gaku-no-mondai-ni-chousen/7-jin-wang-ba-fan-shen-sheno-suan-e は
凡ミスがあるが，解としては正しい。ただし，右上の小円の中心座標はさらに求める必要がある。

ここで示した SymPy による方法では，必要な数値は（一度で）全て得られる。

算額では 463 寸ばかりとなっているが，中円の径は 475.549077332269 寸に間違いない。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
   R = 475.549077332269
   r = 164.545086163906
   y = 739.458905004458
    plot()
    circle(0, 0, 593, :green)
    circle(0, R, R, :blue)
    circle(0, -R, R, :blue)
   circle(593 - r, 0, r)
   circle(-593 + r, 0, r)
   circle(r, y, r)
   circle(-r, y, r)
   circle(r, -y, r)
   circle(-r, -y, r)
    if more
        point(0, R, "R ", :blue, :right)
        point(0, 593, "593 ", :green, :right)
        point(593-r, 0, "593-r ", :red, :top)
        point(r, y, "(r,y)", :red, :top)
        point(0, 0, "0 ", :brown, :top, :right)
        hline!([0], color=:black, lw=0.5)
        vline!([0], color=:black, lw=0.5)
    end
end;