算額（その503）

算額（その503）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8
https://dl.ndl.go.jp/pid/3508431/1/9

長方形内に2本の斜線を引き，分割された領域に甲円，乙円を入れる。
長方形の長辺と短辺の長さが与えられたとき，乙円の直径を求めよ。

長方形の長辺と短辺の長さを 2a, 2r1 とする。
甲円の半径と中心座標を r1, (a - r1, 0)
乙円の半径と中心座標を r2, (0, r1 - r2)
斜線と長方形の上辺の交点を (b, r1)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, b, r1, r2
eq1 = b/sqrt(b^2 + r1^2) - r2/(r1 - r2)
eq2 = r1/(a - r1) - r1/sqrt(b^2 + r1^2);
eq2 = (a - r1) - sqrt(b^2 + r1^2);
res = solve([eq1, eq2], (r2, b))

   2-element Vector{Tuple{Sym, Sym}}:
    (r1*sqrt(a*(a - 2*r1))/(sqrt(a*(a - 2*r1)) - sqrt(a^2 - 2*a*r1 + r1^2)), -sqrt(a*(a - 2*r1)))
    (r1*sqrt(a*(a - 2*r1))/(sqrt(a*(a - 2*r1)) + sqrt(a^2 - 2*a*r1 + r1^2)), sqrt(a*(a - 2*r1)))

2 組の解が得られるが，2 番目のものが適解である。
SymPy では式を簡約化できないが，容易に以下のようにできる。

b = sqrt(a^2 - 2a*r1)
r2 = (r1*b)/(b + a - r1)

長方形の長辺と短辺が 150, 72 のとき，乙円の半径は 10，直径は 20 である。

長 = 150;  平 = 72;  a = 75;  r1 = 36;  b = 15;  r2 = 10
乙円の直径 = 20

using Plots

function draw(長, 平, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1) = (長, 平) ./ 2
   # (r2, b) = (r1*sqrt(a*(a - 2*r1))/(sqrt(a*(a - 2*r1)) + sqrt(a^2 - 2*a*r1 + r1^2)), sqrt(a*(a - 2*r1)))
   b = sqrt(a^2 - 2a*r1)
   r2 = (r1*b)/(b + a - r1)
   @printf("長 = %g;  平 = %g;  a = %g;  r1 = %g;  b = %g;  r2 = %g\n", 2a, 2r1, a, r1, b, r2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot([a, a, -a, -a, a], [-r1, r1, r1, -r1, -r1], color=:black, lw=0.5)
   circle(a - r1, 0, r1)
   circle(r1 - a, 0, r1)
   circle(0, r1 - r2, r2, :blue)
   circle(0, r2 - r1, r2, :blue)
   segment(b, r1, -b, -r1, :green)
   segment(-b, r1, b, -r1, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, "a ", :red, :right, :bottom, delta=delta)
       point(0, r1 - r2, " 乙円:r2,(0,r1-r2)", :blue, :left, :vcenter)
       point(a - r1, 0, "甲円:r1,(a-r1,0)", :red, :center, :bottom, delta=2delta)
       point(b, r1, "(b,r1)", :green, :center, :bottom, delta=2delta)
   else
       plot!(showaxis=false)
   end
end;

算額（その502）

算額（その502）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8

甲円，乙円，丙円が図のように外接しあっている。甲円，乙円の直径が与えられたとき，丙円の直径を求めよ。

乙円の半径と中心座標を r2, (0, 0 )
甲円の半径と中心座標を r1, (r1, y1)
丙円の半径と中心座標を r3, (0, r2 + r3)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r1, y1, r2, r3

eq1 = r1^2 + (y1 - r2 - r3)^2 - (r1 + r3)^2
eq2 = r1^2 + y1^2 - (r1 + r2)^2;
res = solve([eq1, eq2], (y1, r3));

2 組の解が得られるが，2 番目のものが適解である。
簡約化すると以下のようになる。

# y1
res[1][1] |> simplify |> println

   sqrt(r2^3*(2*r1 + r2))/r2

# r3: 乙円の半径
res[1][2] |> simplify |> println

   (r2*(r1 + 2*r2) - 2*sqrt(r2^3*(2*r1 + r2)))/(r1 - 4*r2)

   r1 = 5;  r2 = 7;  y1 = 10.9087;  r3 = 0.857477
   甲円の直径が 10，乙円の直径が 14 のとき，丙円の直径は 1.71495 である。

using Plots

function draw(r1, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y1, r3) = (sqrt(r2^3*(2*r1 + r2))/r2, (r2*(r1 + 2*r2) - 2*sqrt(r2^3*(2*r1 + r2)))/(r1 - 4*r2))
   @printf("r1 = %g;  r2 = %g;  y1 = %g;  r3 = %g\n", r1, r2, y1, r3)
   @printf("甲円の直径が %g，乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r2, 2r3)
   plot()
   circle(r1, y1, r1)
   circle(-r1, y1, r1)
   circle(0, 0, r2, :blue)
   circle(0, r2 + r3, r3, :green) 
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r1, y1, " 甲円:r1\n (r1,y1)", :red, :left, :vcenter)
       point(0, 0.4r2, " 乙円:r2,(0,0)", :blue, :left, mark=false)
       point(0, r2 + r3, " 丙円:r3,(0,r2+r3)", :green, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

算額（その501）

算額（その501）

山形県鶴岡市茨新田 鶴岡神明宮 文久2年(1862)
http://www.wasan.jp/yamagata/sinmei.html

外円内に 4 本の斜線を引き，区画された領域に大円，中円，小円を入れる。
外円，大円，小円の直径がそれぞれ 125寸，99寸，21寸のとき，中円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0 )
大円の半径と中心座標を r1, (0, r1 - r0)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (0, r0 - r3)
斜線と外円の交点座標を (a, sqrt(r0^2 - a^2), (b, sqrt(r0^2 - b^2)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms a, b, r0, r1, r2, x2, y2, r3

eq1 = distance(a, sqrt(r0^2 - a^2), -b, sqrt(r0^2 - b^2), 0, r1 - r0) - r1^2
eq2 = distance(a, sqrt(r0^2 - a^2), -b, sqrt(r0^2 - b^2), x2, y2) - r2^2
eq3 = distance(a, sqrt(r0^2 - a^2), -b, sqrt(r0^2 - b^2), 0, r0 -r3) - r3^2
eq4 = distance(-a, sqrt(r0^2 - a^2), b, sqrt(r0^2 - b^2), x2, y2) - r2^2
eq5 = distance(a, sqrt(r0^2 - a^2), b, sqrt(r0^2 - b^2), x2, y2) - r2^2;

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, r2, x2, y2) = u
   t1 = sqrt(-a^2 + r0^2)
   t2 = sqrt(-b^2 + r0^2)
   return [
       -r1^2 + (t1 - t2)^2*(a^2*b - a*b^2 + a*r0^2 - a*r0*t1 + a*r0*t2 + a*r1*t1 - a*r1*t2 - a*t1*t2 - b*r0^2 - b*r0*t1 + b*r0*t2 + b*r1*t1 - b*r1*t2 + b*t1*t2)^2/(4*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)^2) + (-r0 + r1 - (-a^3*b - a^2*r0^2 + a^2*r0*t1 - 3*a^2*r0*t2 - a^2*r1*t1 + 3*a^2*r1*t2 + a^2*t1*t2 + a*b^3 + b^2*r0^2 + 3*b^2*r0*t1 - b^2*r0*t2 - 3*b^2*r1*t1 + b^2*r1*t2 - b^2*t1*t2 - 4*r0^3*t1 + 4*r0^3*t2 + 4*r0^2*r1*t1 - 4*r0^2*r1*t2)/(2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)))^2,  # eq1
       -r2^2 + (x2 - (x2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) + (t1 - t2)*(a^2*b + a^2*x2 - a*b^2 + a*r0^2 + a*y2*t1 - a*y2*t2 - a*t1*t2 + b^2*x2 - b*r0^2 + b*y2*t1 - b*y2*t2 + b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2 + (y2 - (y2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) - (a + b)*(a^2*b + a^2*x2 - a*b^2 + a*r0^2 + a*y2*t1 - a*y2*t2 - a*t1*t2 + b^2*x2 - b*r0^2 + b*y2*t1 - b*y2*t2 + b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2,  # eq2
       -r3^2 + (t1 - t2)^2*(a^2*b - a*b^2 + a*r0^2 + a*r0*t1 - a*r0*t2 - a*r3*t1 + a*r3*t2 - a*t1*t2 - b*r0^2 + b*r0*t1 - b*r0*t2 - b*r3*t1 + b*r3*t2 + b*t1*t2)^2/(4*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)^2) + (r0 - r3 - (-a^3*b - a^2*r0^2 - a^2*r0*t1 + 3*a^2*r0*t2 + a^2*r3*t1 - 3*a^2*r3*t2 + a^2*t1*t2 + a*b^3 + b^2*r0^2 - 3*b^2*r0*t1 + b^2*r0*t2 + 3*b^2*r3*t1 - b^2*r3*t2 - b^2*t1*t2 + 4*r0^3*t1 - 4*r0^3*t2 - 4*r0^2*r3*t1 + 4*r0^2*r3*t2)/(2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)))^2,  # eq3
       -r2^2 + (x2 - (x2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) + (t1 - t2)*(-a^2*b + a^2*x2 + a*b^2 - a*r0^2 - a*y2*t1 + a*y2*t2 + a*t1*t2 + b^2*x2 + b*r0^2 - b*y2*t1 + b*y2*t2 - b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2 + (y2 - (y2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) + (a + b)*(-a^2*b + a^2*x2 + a*b^2 - a*r0^2 - a*y2*t1 + a*y2*t2 + a*t1*t2 + b^2*x2 + b*r0^2 - b*y2*t1 + b*y2*t2 - b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2,  # eq4
       -r2^2 + (x2 - (a*(a*b - r0^2 + t1*t2) - (a - b)*(a*b + a*x2 - b*x2 - r0^2 + y2*t1 - y2*t2 + t1*t2)/2)/(a*b - r0^2 + t1*t2))^2 + (y2 - (2*t1*(a*b - r0^2 + t1*t2) - (t1 - t2)*(a*b + a*x2 - b*x2 - r0^2 + y2*t1 - y2*t2 + t1*t2))/(2*(a*b - r0^2 + t1*t2)))^2,  # eq5
   ]
end;
(r0, r1, r3) = (125, 99, 21) .// 2
iniv = BigFloat[60.0, 33, 10, 29, 40]
res = nls(H, ini=iniv)

  (BigFloat[60.5769230769230769230769230769230769230769230769230769230769230769230769357718, 31.73076923076923076923076923076923076923076923076923076923076923076923419539377, 10.81730769230769230769230769230769230769230769230769230769230769230769343490884, 28.12500000000000000000000000000000000000000000000000000000000000000000245224407, 40.62499999999999999999999999999999999999999999999999999999999999999999867465516], true)

   r0 = 62.5;  r1 = 49.5;  r3 = 10.5
   a = 60.5769;  b = 31.7308;  r2 = 10.8173;  x2 = 28.125;  y2 = 40.625
   中円の直径 = 2r2 = 21.6346

中円の直径は 21.6346 ほどである。術では「二十一寸五十二分...」とあるので，一致しない。

実際の図は，算額に描かれたものとかなり異なる。
このようなことは算額においてはよくあることであるが，原因の一つは算額の図における各円の直径の割合が問に書かれたものと違うこと。
さらに，この算額の場合には外円の直径が「一百二十五寸」と書かれているのだが，「一」の部分が不鮮明であること。「二」でも「三」でもなさそうではあるのだが。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r3) = (125, 99, 21) .// 2
   (a, b, r2, x2, y2) = res[1]
   @printf("r0 = %g;  r1 = %g;  r3 = %g\n", r0, r1, r3)
   @printf("a = %g;  b = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", a, b, r2, x2, y2)
   @printf("中円の直径 = 2r2 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :black) 
   circle(0, r1 - r0, r1, :red) 
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(0, r0 - r3, r3, :magenta)
   plot!([a, b, -a], [sqrt(r0^2 - a^2), sqrt(r0^2 - b^2), sqrt(r0^2 - a^2)], color=:brown, lw=0.5)
   plot!([-a, -b, a], [sqrt(r0^2 - a^2), sqrt(r0^2 - b^2), sqrt(r0^2 - a^2)], color=:brown, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1 - r0, " 大円:r1,(0,r1-r0)", :red, :left, :vcenter)
       point(x2, y2, " 中円:r2,(x2,y2)", :blue, :center, :top, delta=-delta)
       point(0, r0 - r3, " 小円:r3\n (0,r0-r3)", :magenta, :left, :vcenter)
       point(r0, 0, "r0 ", :black, :right, :bottom)
       point(a, sqrt(r0^2 - a^2), "(a,√(r0^2-a^2)", :brown, :right)
       point(b, sqrt(r0^2 - b^2), " (b,√(r0^2-b^2)", :brown, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

算額（その500）

算額（その500）

群馬県高崎市 清水寺 文政11年(1828)
http://takasakiwasan.web.fc2.com/gunnsann/pdf/046_2monndai.pdf

長径 7 寸，短径 5 寸の菱形の中に楕円と小円 4 個が入っている。小円の直径はいかほどか。

菱形の中心を原点とし，菱形の長径，短径を 2x，2y とする。
楕円の長径，短径を 2a，2b とする。
等円の半径と中心座標を r1, (a + r, 0), (0, b + r)
として以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x, y, a, b, r, x1, y1

z = sqrt(x^2 + y^2)
eq1 = r/(x - a - r) - y/z  # 半径と短径の関係
eq2 = r/(y - b - r) - x/z  # 半径と長径の関係
eq3 = (b^2*x1)/(a^2*y1) - y/x  # 接線の傾きに関して
eq4 = b^2/y1 - y  # 接線の切片に関して
eq5 = x1^2/a^2 + y1^2/b^2 - 1;  # 接点が楕円状にあることに関して
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, r, 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, r, x1, y1) = u
   return [
       r/(-a - r + x) - y/sqrt(x^2 + y^2),  # eq1
       r/(-b - r + y) - x/sqrt(x^2 + y^2),  # eq2
       -y/x + b^2*x1/(a^2*y1),  # eq3
       b^2/y1 - y,  # eq4
       -1 + y1^2/b^2 + x1^2/a^2,  # eq5
   ]
end;
(x, y) = (7, 5) .// 2
iniv = BigFloat[2.55, 1.7, 0.35, 2, 1] .* (x/3.5)
res = nls(H, ini=iniv)

   (BigFloat[2.54414717679825272253216901118052540925755526167623360938083320066487243514035, 1.716860438920218883208445479192299273530396028742001559200850725830125354563721, 0.3513564057748656326259665107244486935073587532847506507495566147571622830042101, 1.849338530631577084618933179122660657165347568695981336767010069217909175976173, 1.179043906691730653843619157769528102024751736645727616594992765005302482977008], true)

x = 3.5;  y = 2.5;  a = 2.54415;  b = 1.71686;  r = 0.351356
円の直径 = 2r = 0.702713

x = 3.5, y = 2.5 のとき，円の直径は 0.702713 である。

算額の図を再現している（かどうかは怪しいが），元の記事に掲載されている図は，実際のものとはかなり異なっている。

一番の問題は，図では分かりづらいが実測してみると左右の円と上下の円の直径が異なっていることである。比は 1.125 程度も違う。
また，菱形の長径と短径の比は，問では 7/5=1.4 であるが，図では 1.64 程度である。
もし菱形の長径と短径が条件通り 7, 5 で，円の直径が 1 とすれば，楕円は菱形に内接しない。

いずれにしろ，得られた数値で図形を描くと条件を満たす図になっているので，得られた結果として「円の直径は 0.702713 である」は正しいと思われる。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, y) = (7, 5) ./ 2
   (a, b, r, x1, y1) = res[1]
   @printf("x = %g;  y = %g;  a = %g;  b = %g;  r = %g\n", x, y, a, b, r)
   @printf("円の直径 = 2r = %g\n", 2r)
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(0, b + r, r) 
   circle(0, -b - r, r) 
   circle(a + r, 0, r)
   circle(-a - r, 0, r)
   ellipse(0, 0, a, b, color=:blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, y1, " (x1,y1)", :blue, :left, :vcenter)
       point(x, 0, " x", :black, :left, :bottom, delta=delta/2)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta/2)
       point(a + r, 0, "a+r", :red, :center, :top, delta=-delta/2)
       point(0, y, " y", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :top, delta=-delta/2)
       point(0, b + r, " b+r", :red, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;