算額（その409）

算額（その409）

兵庫県高砂市 生石神社 明治9年(1876)
兵庫県高砂市 高砂神社

算額(その8) : 高砂神社への算額奉納
https://toyo.repo.nii.ac.jp/?action=repository_action_common_download&item_id=2486&item_no=1&attribute_id=18&file_no=1

第3問　図のように，外円に甲円が内接している。その接点を通る任意の点を引き，外円，甲円，弦に接する乙円を描く。さらに，弦に関して乙円と反対側で，甲円と乙円に接する丙円を描く（すなわち，弦と乙円と丙円は，共通の接点を持つ）。
ここで，外円，甲円，乙円の直径が与えられたときに，丙円の径を求めよ。

「弦と乙円と丙円は，共通の接点を持つ」ということは，乙円と丙円の中心を結ぶ線分は弦と垂直に交わるということである。
また，甲円と丙円の中心を結ぶ線分と甲円と丙円の共通接線も直角に交わるのであるから，甲円と丙円の中心を結ぶ直線と，乙円と丙円の中心を結ぶ直線も直行するので，3つの円の中心を結んでできる三角形は直角三角形である（間違えてないですよね）。ということで，普通の4つの円の内外接問題 eq2~eq4 に，中心間の距離におけるピタゴラスの定理を適用する eq1 を加えて連立方程式を解けばよいはず。ただし，これによる解はちゃんと図に描けるが，数値が術と一致しない。

外円の中心を原点に置く。
外円の半径と中心座標を　r0, (0, 0)
甲円の半径と中心座標を　r1, (0, r0 - r1)
乙円の半径と中心座標を　r2, (x2, y2)
丙円の半径と中心座標を　r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive,
     r2::positive, x2::positive, y2::negative,
     r3::positive, x3::positive, y3::negative;

(r0, r1, r2) = (10, 6, 3)
eq1 = (r1 + r3)^2 + (r2 + r3)^2 - (r1 + r2)^2
eq2 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq3 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq4 = x3^2 + (r0 - r1 - y3)^2 - (r1 + r3)^2
eq5 = x2^2 + y2^2 - (r0 - r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r3, x3, y3, x2, y2))

   1-element Vector{NTuple{5, Sym}}:
    (-17*sqrt(8*sqrt(85) + 97)/42 - 9/2 + 2*sqrt(680*sqrt(85) + 8245)/21, -3*sqrt(40*sqrt(85) + 485)/14 + 3*sqrt(5)/2 + 2*sqrt(136*sqrt(85) + 1649)/21, 1 - sqrt(8*sqrt(85) + 97)/3, 3*sqrt(5), -2)

   r3 = 1.68466;  x3 = 2.22403;  y3 = -3.35579;  x2 = 6.7082;  y2 = -2

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r2) = (10, 6, 3)
   (r3, x3, y3, x2, y2) = res[1]
   @printf("r3 = %g;  x3 = %g;  y3 = %g;  x2 = %g;  y2 = %g\n", r3, x3, y3, x2, y2)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x3, y3, r3, :green)
   plot!([0, x2, x3, 0], [r0 - r1, y2, y3, r0 - r1], color=:black, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0, "r0 ", :red, :right, :bottom)
       point(0, r0 - r1, "甲円 \nr0-r1 ", :red, :right, :vcenter)
       point(x2, y2, "乙円:r2\n(x2,y2)", :blue)
       point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, delta=-delta)
       segment(0, r0, (r3*x2 + r2*x3)/(r2 + r3), (r3*y2 + r2*y3)/(r2 + r3), :orange)
       point((r3*x2 + r2*x3)/(r2 + r3), (r3*y2 + r2*y3)/(r2 + r3), "")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その408）

算額（その408）

千葉県成田市 成田山新勝寺 明治30年
https://toyo.repo.nii.ac.jp/?action=repository_action_common_download&item_id=2534&item_no=1&attribute_id=18&file_no=1

三角形内に界斜を引き，2 区分された領域に等円を入れる。中斜，小斜，界斜それぞれが 11 寸，6 寸，6 寸のとき，大斜はいかほどか。

数値が与えれて，求めるものは違うが，算額（その406），算額（その407）と同じである。

大斜，中斜，小斜，界斜を変数名として，
等円の半径と中心座標を r, (x1, r), (x2, r)
三角形の頂点座標を (x, y)
界斜と大斜の交点座標を (a, 0) として，以下の連立方程式を立て，nlsolve() により数値解を求める。

include("julia-source.txt");

using SymPy

@syms r::positive, a::positive, x::positive, y::positive, x1::positive,
     x2::positive, 大斜::positive, 中斜::positive, 小斜::positive, 界斜::positive;

eq1 = (界斜 + 中斜 + a)r + (界斜 + 小斜 + 大斜 - a)r - 大斜*y
eq2 = x^2 + y^2 - 中斜^2
eq3 = (大斜 - x)^2 + y^2 - 小斜^2
eq4 = distance(a, 0, x, y, x1, r) - r^2
eq5 = distance(a, 0, x, y, x2, r) - r^2
eq6 = distance(0, 0, x, y, x1, r) - r^2
eq7 = distance(大斜, 0, x, y, x2, r) - r^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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)
   (大斜, r, a, x, y, x1, x2) = u
   return [
       r*(a + 中斜 + 界斜) + r*(-a + 大斜 + 小斜 + 界斜) - y*大斜,  # eq1
       x^2 + y^2 - 中斜^2,  # eq2
       y^2 - 小斜^2 + (-x + 大斜)^2,  # eq3
       -r^2 + (r - y*(a^2 - a*x - a*x1 + r*y + x*x1)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x1 - (x1*(a^2 - 2*a*x + x^2 + y^2) - y*(a*r - a*y - r*x + x1*y))/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq4
       -r^2 + (r - y*(a^2 - a*x - a*x2 + r*y + x*x2)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x2 - (x2*(a^2 - 2*a*x + x^2 + y^2) - y*(a*r - a*y - r*x + x2*y))/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq5
       -r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq6
       -r^2 + (r - y*(r*y + x*x2 - x*大斜 - x2*大斜 + 大斜^2)/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2 + (x2 - (x2*(x^2 - 2*x*大斜 + y^2 + 大斜^2) + y*(r*x - r*大斜 - x2*y + y*大斜))/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2,  # eq7
   ]
end;

(中斜, 小斜, 界斜) = (11, 6, 4)
iniv = [big"14.0", 1.5, 10, 10.2, 4, 7.5, 11]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[14.99999999999999999999999999999999999999999999999999999999989878381102949411463, 1.414213562373095048801688724209698078569671875376948073176727771489368647748441, 8.999999999999999999999999999999999999999999999999999999999912693616410060723697, 10.33333333333333333333333333333333333333333333333333333333329944227882266388349, 3.771236166328253463471169931225861542852458334338528195137946162035145601912081, 7.999999999999999999999999999999999999999999999999999999999882812339426680153511, 10.99999999999999999999999999999999999999999999999999999999998515232648292937302], true)

中斜, 小斜, 界斜 がそれぞれ 11, 6, 4 寸のとき，
大斜 = 15;  r = 1.41421;  a = 9;  x = 10.3333;  y = 3.77124;  x1 = 8;  x2 = 11 である。
算額の術では 大斜 = sqrt((中斜 + 小斜)^2 - 4界斜^2) となっており，正しい答えが求まっていることがわかる。

(中斜, 小斜, 界斜) = (11, 6, 4)
sqrt((中斜 + 小斜)^2 - 4界斜^2)

   15.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (中斜, 小斜, 界斜) = (11, 6, 4)
   (大斜, r, a, x, y, x1, x2) = res[1]
   @printf("大斜 = %g;  r = %g;  a = %g;  x = %g;  y = %g;  x1 = %g;  x2 = %g\n", 大斜, r, a, x, y, x1, x2)
   plot([0, 大斜, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x1, r, r)
   circle(x2, r, r)
   segment(a, 0, x, y, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(大斜, 0, " 大斜", :black, :left, :bottom, delta=delta)
       point(x, y, " (x,y)", :black, :left, :vcenter, delta=-delta)
       point(x1, r, "(x1,r)", :red, :center, delta=-delta)
       point(x2, r, "(x2,r)", :red, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(x/2, 2y/3, "中斜  ", :black, :left, delta=3delta, mark=false)
       point((大斜 + x)/2, 2y/3, "小斜", :black, :left, mark=false)
       point((a + x)/2, 4y/5, "  界斜", :blue, :right, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       #plot!(xlims=(0, 350))
   else
       plot!(showaxis=false)
   end
end;

算額（その407）

算額（その407）

埼玉県神川町 光明寺 大正3年(1914)
https://toyo.repo.nii.ac.jp/?action=repository_action_common_download&item_id=2534&item_no=1&attribute_id=18&file_no=1

三角形内に界斜を引き，2 区分された領域に等円を入れる。大斜，中斜，小斜それぞれが 315 寸，257 寸，68 寸のとき，界斜はいかほどか。

数値が与えれているが，算額（その406）と同じである。

大斜，中斜，小斜，界斜を変数名として，
等円の半径と中心座標を r, (x1, r), (x2, r)
三角形の頂点座標を (x, y)
界斜と大斜の交点座標を (a, 0) として，以下の連立方程式を立て，nlsolve() により数値解を求める。

include("julia-source.txt");

using SymPy

@syms r::positive, a::positive, x::positive, y::positive, x1::positive,
     x2::positive, 大斜::positive, 中斜::positive, 小斜::positive, 界斜::positive;

eq1 = (界斜 + 中斜 + a)r + (界斜 + 小斜 + 大斜 - a)r - 大斜*y
eq2 = x^2 + y^2 - 中斜^2
eq3 = (大斜 - x)^2 + y^2 - 小斜^2
eq4 = distance(a, 0, x, y, x1, r) - r^2
eq5 = distance(a, 0, x, y, x2, r) - r^2
eq6 = distance(0, 0, x, y, x1, r) - r^2
eq7 = distance(大斜, 0, x, y, x2, r) - r^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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)
   (界斜, r, a, x, y, x1, x2) = u
   return [
       r*(a + 中斜 + 界斜) + r*(-a + 大斜 + 小斜 + 界斜) - y*大斜,  # eq1
       x^2 + y^2 - 中斜^2,  # eq2
       y^2 - 小斜^2 + (-x + 大斜)^2,  # eq3
       -r^2 + (r - y*(a^2 - a*x - a*x1 + r*y + x*x1)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x1 - (x1*(a^2 - 2*a*x + x^2 + y^2) - y*(a*r - a*y - r*x + x1*y))/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq4
       -r^2 + (r - y*(a^2 - a*x - a*x2 + r*y + x*x2)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x2 - (x2*(a^2 - 2*a*x + x^2 + y^2) - y*(a*r - a*y - r*x + x2*y))/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq5
       -r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq6
       -r^2 + (r - y*(r*y + x*x2 - x*大斜 - x2*大斜 + 大斜^2)/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2 + (x2 - (x2*(x^2 - 2*x*大斜 + y^2 + 大斜^2) + y*(r*x - r*大斜 - x2*y + y*大斜))/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2,  # eq7
   ]
end;

(大斜, 中斜, 小斜) = (315, 257, 68)
iniv = [big"35.0", 15, 220, 250, 30, 220, 250]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[40.00000000000000000000000000000000000000000000000000000000000000000000012868471, 13.99999999999999999999999999999999999999999999999999999999999999999999999479218, 231.0000000000000000000000000000000000000000000000000000000000000000000001855419, 255.0, 31.99999999999999999999999999999999999999999999999999999999999999999999999999945, 223.9999999999999999999999999999999999999999999999999999999999999999999999166771, 259.0000000000000000000000000000000000000000000000000000000000000000000000208086], true)

大斜, 中斜, 小斜 がそれぞれ 4315, 257, 68 寸のとき，
界斜 = 40;  r = 14;  a = 231;  x = 255;  y = 32;  x1 = 224;  x2 = 259である。
算額の術では 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2 となっており，正しい答えが求まっていることがわかる。

(大斜, 中斜, 小斜) = (315, 257, 68)
sqrt((中斜 + 小斜)^2 - 大斜^2)/2

   40.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (大斜, 中斜, 小斜) = (315, 257, 68)
   (界斜, r, a, x, y, x1, x2) = res[1]
   #(界斜, r, a, x, y, x1, x2) = [40, 50, 150, 150, 200, 103, 210]
   @printf("界斜 = %g;  r = %g;  a = %g;  x = %g;  y = %g;  x1 = %g;  x2 = %g\n", 界斜, r, a, x, y, x1, x2)
   plot([0, 大斜, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x1, r, r)
   circle(x2, r, r)
   segment(a, 0, x, y, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(大斜, 0, " 大斜", :black, :left, :bottom, delta=delta)
       point(x, y, " (x,y)", :black, :left, :vcenter, delta=-delta)
       point(x1, r, "(x1,r)", :red, :center, delta=-delta)
       point(x2, r, "(x2,r)", :red, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(x/2, 2y/3, "中斜", :black, :left, mark=false)
       point((大斜 + x)/2, 2y/3, "小斜", :black, :left, mark=false)
       point((a + x)/2, 2y/3, "  界斜", :blue, :left, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       #plot!(xlims=(0, 5))
   else
       plot!(showaxis=false)
   end
end;

算額（その406）

算額（その406）

兵庫県伊丹市 昆陽寺 嘉永7年(1854)，昭和43年(1968)復元奉納
https://toyo.repo.nii.ac.jp/?action=repository_action_common_download&item_id=2534&item_no=1&attribute_id=18&file_no=1

第7問，三角形内に界斜を引き，2 区分された領域に等円を入れる。大斜，中斜，小斜が与えられたとき，界斜はいかほどか。

大斜，中斜，小斜，界斜を変数名として，
等円の半径と中心座標を r, (x1, r), (x2, r)
三角形の頂点座標を (x, y)
界斜と大斜の交点座標を (a, 0) として，以下の連立方程式を立て，nlsolve() により数値解を求める。

include("julia-source.txt");

using SymPy

@syms r::positive, a::positive, x::positive, y::positive, x1::positive,
     x2::positive, 大斜::positive, 中斜::positive, 小斜::positive, 界斜::positive;

eq1 = (界斜 + 中斜 + a)r + (界斜 + 小斜 + 大斜 - a)r - 大斜*y
eq2 = x^2 + y^2 - 中斜^2
eq3 = (大斜 - x)^2 + y^2 - 小斜^2
eq4 = distance(a, 0, x, y, x1, r) - r^2
eq5 = distance(a, 0, x, y, x2, r) - r^2
eq6 = distance(0, 0, x, y, x1, r) - r^2
eq7 = distance(大斜, 0, x, y, x2, r) - r^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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)
   (界斜, r, a, x, y, x1, x2) = u
   return [
       r*(a + 中斜 + 界斜) + r*(-a + 大斜 + 小斜 + 界斜) - y*大斜,  # eq1
       x^2 + y^2 - 中斜^2,  # eq2
       y^2 - 小斜^2 + (-x + 大斜)^2,  # eq3
       -r^2 + (r - y*(a^2 - a*x - a*x1 + r*y + x*x1)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x1 - (a^2*x1 - a*r*y - 2*a*x*x1 + a*y^2 + r*x*y + x^2*x1)/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq4
       -r^2 + (r - y*(a^2 - a*x - a*x2 + r*y + x*x2)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x2 - (a^2*x2 - a*r*y - 2*a*x*x2 + a*y^2 + r*x*y + x^2*x2)/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq5
       -r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq6
       -r^2 + (r - y*(r*y + x*x2 - x*大斜 - x2*大斜 + 大斜^2)/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2 + (x2 - (x2*(x^2 - 2*x*大斜 + y^2 + 大斜^2) + y*(r*x - r*大斜 - x2*y + y*大斜))/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2,  # eq7   
   ]
end;

(大斜, 中斜, 小斜) = (9//2, 4, 3)
iniv = [big"3.0", 1.1, 2.5, 3, 3, 2, 4]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[2.680951323690902076203539918038571176272552963887644954668323575268635315149453, 0.6975859391932934831098526753294690788534179473670310788732592770458445561358907, 2.43201081695757731639921335154698418305054378580278996601587232707696243601583, 3.027777777777777777777777777777777777777777777777777777777777777777777777777774, 2.613916932191048363498370393355647818939596565227754137625588266018177169955276, 1.875529746633337620097836716754206503388995410957580351169949221547302218391068, 3.306481070324239696301376634792777679661548374845176108710264222377909400696189], true)

大斜, 中斜, 小斜 がそれぞれ 4.5, 4, 3 のとき，
界斜 = 2.68095;  r = 0.697586;  a = 2.43201;  x = 3.02778;  y = 2.61392;  x1 = 1.87553;  x2 = 3.30648
である。
算額の術では 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2 となっており，正しい答えが求まっていることがわかる。

(大斜, 中斜, 小斜) = (9//2, 4, 3)
sqrt((中斜 + 小斜)^2 - 大斜^2)/2

   2.680951323690902

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (大斜, 中斜, 小斜) = (9//2, 4, 3)
   (界斜, r, a, x, y, x1, x2) = res[1]
   @printf("界斜 = %g;  r = %g;  a = %g;  x = %g;  y = %g;  x1 = %g;  x2 = %g\n", 界斜, r, a, x, y, x1, x2)
   plot([0, 大斜, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x1, r, r)
   circle(x2, r, r)
   segment(a, 0, x, y, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(大斜, 0, " 大斜", :black, :left, :bottom, delta=delta)
       point(x, y, " (x,y)", :black, :left, :vcenter, delta=-delta)
       point(x1, r, "(x1,r)", :red, :center, delta=-delta)
       point(x2, r, "(x2,r)", :red, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(x/2, 2y/3, "中斜", :black, :left, mark=false)
       point((大斜 + x)/2, 2y/3, "小斜", :black, :left, mark=false)
       point((a + x)/2, 2y/3, "  界斜", :blue, :left, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その405）

算額（その405）

長野県北佐久郡軽井沢町峠 熊野神社 明治5年(1872)
山口正義：「やまぶき3」
https://yamabukiwasan.sakura.ne.jp/ymbk42.pdf

外円の中に，元，利，貞，亨，無名の円が入っている。外円の直径が 16 寸のとき，亨円の直径を求めよ。

外円の中心を原点に置く。
外円の半径と中心座標を r0, (0, 0)
元円の半径と中心座標を r1, (0, r0 - r1) と (0, r0 - 3r1)
利円の半径と中心座標を r2, (x2, y2)
貞円の半径と中心座標を r3, (x3, y3)
亨円の半径と中心座標を r4, (0, r4 - r0)
無名円の半径と中心座標を r5, (x5, y5)
として連立方程式を立て nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive, r4::positive,
     r5::positive, x5::positive, y5::positive;

r0 = 16//2
eq1 = x2^2 + y2^2 - (r0 - r2)^2
eq2 = x3^2 + y3^2 - (r0 - r3)^2
eq3 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq4 = x5^2 + (r0 - 3r1 - y5)^2 - (r1 + r5)^2
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq7 = x3^2 + (r4 - r0 - y3)^2 - (r3 + r4)^2
eq8 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq9 = x5^2 + (r4 - r0 - y5)^2 - (r4 + r5)^2
eq10 = 4r1 + 2r4 - 2r0
eq11 = r0 - 2r1 - y2
# res = solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11))

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)
   (r1, r2, x2, y2, r3, x3, y3, r4, r5, x5, y5) = u
   return [
x2^2 + y2^2 - (8 - r2)^2,  # eq1
x3^2 + y3^2 - (8 - r3)^2,  # eq2
x2^2 - (r1 + r2)^2 + (-r1 - y2 + 8)^2,  # eq3
x5^2 - (r1 + r5)^2 + (-3*r1 - y5 + 8)^2,  # eq4
-(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq5
-(r2 + r5)^2 + (x2 - x5)^2 + (y2 - y5)^2,  # eq6
x3^2 - (r3 + r4)^2 + (r4 - y3 - 8)^2,  # eq7
-(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq8
x5^2 - (r4 + r5)^2 + (r4 - y5 - 8)^2,  # eq9
4*r1 + 2*r4 - 16,  # eq10
-2*r1 - y2 + 8,  # eq11
      ]
end;
r0 = 16//2
iniv = [big"1.9", 2.49, 4.3, 3.37, 2.29, 5.61, -1.28, 3.9, 1.1, 2.35, 0.34]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.750000000000000000000000000000000000000000001683866081422877969875418056451415, 2.243589743589743589743589743589743589743589745370092093352019975896094213312947, 3.589743589743589743589743589743589743589743594636276696631029402293443881749929, 4.499999999999999999999999999999999999999999996632267837154244060249163887074371, 2.348993288590604026845637583892617449664429532022065401438839115069164780141786, 5.637583892617449664429530201342281879194630878160421974850153126915301979093254, 0.3892617449664429530201342281879194630872483133777124317852969056382365325617382, 4.499999999999999999999999999999999999999999996632267837154244060249163887074371, 1.017441860465116279069767441860465116279069768415921597062855408025815041846539, 2.44186046511627906976744186046511627906976744452530995904551667578153181660868, 1.447674418604651162790697674418604651162790691173766512477690728560896338607614], true)

r1 = 1.75;  r2 = 2.24359;  x2 = 3.58974;  y2 = 4.5;  r3 = 2.34899;  x3 = 5.63758;  y3 = 0.389262;  r4 = 4.5;  r5 = 1.01744;  x5 = 2.44186;  y5 =1.44767

亨円の直径は 9 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 16//2
   (r1, r2, x2, y2, r3, x3, y3, r4, r5, x5, y5) = res[1]
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  r5 = %g;  x5 = %g;  y5 =%g\n", r1, r2, x2, y2, r3, x3, y3, r4, r5, x5, y5)
   @printf("亨円の直径 = %g\n", 2r4)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r0 - 3r1, r1, :blue)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
   circle(0, r4 - r0, r4, :orange)
   circle(x5, y5, r5, :brown)
   circle(-x5, y5, r5, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0 - r1, " r0-r1", :blue)
       point(0.1r1, r0 - r1, " 元", :blue, :left, :bottom, delta=delta, mark=false)
       point(0, r0 - 3r1, " r0-3r1", :blue)
       point(x2, y2, " 利:r2,(x2,y2)", :green, :center, delta=-delta)
       point(x3, y3, " 貞:r3,(x3,y3)", :magenta, :center, :bottom, delta=delta/2)
       point(0, r4 - r0, " r4-r0  亨", :orange)
       point(x5, y5, " 無:r5,(x5,y5)", :brown, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その404）

算額（その404）

法道寺善の算変法「観術」
永井信一：早数教分科会レポート「和算と反転法について」
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/hanten.pdf

外円の中に 5 個の円が入っている。東円，西円，南円の直径がそれぞれ 3 寸，1 寸，2 寸のとき，外円，北円，心円の直径を求めよ。

外円の中心を原点に置き，西円の中心が x 軸上に来るように円を配置する。
東円の半径と中心座標を r1, (x1, y1)
西円の半径と中心座標を r2, (r0 - r2, 0)
南円の半径と中心座標を r3, (x3, y3)
北円の半径と中心座標を r4, (x4, y4)
心円の半径と中心座標を r5, (x5, y5)
外円の半径と中心座標を r0, (0, 0)
として連立方程式を立て nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, x1::negative, y1::negative,
     r2::positive, r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::negative,
     r5::positive, x5::positive, y5::negative;

(r1, r2, r3) = (3, 1, 2) .// 2
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x3^2 + y3^2 - (r0 - r3)^2
eq3 = x4^2 + y4^2 - (r0 - r4)^2
eq4 = (x3 - x1)^2 + (y3 - y1)^2 - (r3 + r1)^2
eq5 = (x4 - x1)^2 + (y4 - y1)^2 - (r4 + r1)^2
eq6 = (x5 - x1)^2 + (y5 - y1)^2 - (r5 + r1)^2
eq7 = (r0 - r2 - x3)^2 + y3^2 - (r2 + r3)^2
eq8 = (r0 - r2 - x4)^2 + y4^2 - (r2 + r4)^2
eq9 = (r0 - r2 - x5)^2 + y5^2 - (r2 + r5)^2
eq10 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq11 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2;
# res = solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11))

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)
   (r0, x1, y1, x3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
       x1^2 + y1^2 - (r0 - 3/2)^2,  # eq1
       x3^2 + y3^2 - (r0 - 1)^2,  # eq2
       x4^2 + y4^2 - (r0 - r4)^2,  # eq3
       (-x1 + x3)^2 + (-y1 + y3)^2 - 25/4,  # eq4
       -(r4 + 3/2)^2 + (-x1 + x4)^2 + (-y1 + y4)^2,  # eq5
       -(r5 + 3/2)^2 + (-x1 + x5)^2 + (-y1 + y5)^2,  # eq6
       y3^2 + (r0 - x3 - 1/2)^2 - 9/4,  # eq7
       y4^2 - (r4 + 1/2)^2 + (r0 - x4 - 1/2)^2,  # eq8
       y5^2 - (r5 + 1/2)^2 + (r0 - x5 - 1/2)^2,  # eq9
       -(r5 + 1)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq10
       -(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2,  # eq11
   ]
end;
iniv = [big"2.5", -0.9, -0.4, 0.9, 1.2, 0.6, 1.2, -1.3, 0.32, 1.0, -0.3]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[2.508041569829081421977001552141168300225929489745528932614456499113765038860902, -0.485951405752271286255384816114094927579869826307774264025836820627077860874849, -0.8831755418663212924359755945019721217093601316808355690943080396556492904135786, 1.010043911301963852566206096056080557623999906394429728983586117023084099186487, 1.119821715084321284656363526042374037999885992006369003023642241853344206655946, 0.6000000000000000000000000000000000000000012240791396860081807414352057058378827, 1.609242974712810880330524278490115654664770306212569723396079629646188139868433, -1.025163245797121287768208353426213271483677234175174482387089368453727261649156, 0.3262233880108996037951833703229864820441137963024028723844040481034973892347998, 1.20445443902395988951867252073114035193644842415907539462550103774599845159441, -0.1920750116506622924019613444332251808267586158530185578363832744817135812429724], true)

   r0 = 2.50804;  x1 = -0.485951;  y1 = -0.883176;  x3 = 1.01004;  y3 = 1.11982;  r4 = 0.6;  x4 = 1.60924;  y4 = -1.02516;  r5 = 0.326223;  x5 = 1.20445;  y5 = -0.192075
   北径 = 1.2;  外径 = 5.01608;  心径 = 0.652447

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (3, 1, 2) .// 2
   (r0, x1, y1, x3, y3, r4, x4, y4, r5, x5, y5) = res[1]
   @printf("r0 = %g;  x1 = %g;  y1 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", r0, x1, y1, x3, y3, r4, x4, y4, r5, x5, y5)
   @printf("北径 = %g;  外径 = %g;  心径 = %g\n", 2r4, 2r0, 2r5)  
   plot()
   circle(0, 0, r0, :black)
   circle(x1, y1, r1)
   circle(r0 - r2, 0, r2, :blue)
   circle(x3, y3, r3, :orange)
   circle(x4, y4, r4, :magenta)
   circle(x5, y5, r5, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, " 東:r1,(x1,y1)", :red, :center, delta=-delta)
       point(r0 - r2, 0, " 西:r2\n(r0-r2,0)", :blue, :center, :bottom, delta=delta)
       point(x3, y3, " 南:r3,(x3,y3)", :orange, :center, delta=-delta)
       point(x4, y4, " 北:r4,(x4,y4)", :magenta, :center, delta=-delta)
       point(x5, y5, " 心:r5\n(x5,y5)", :brown, :center, :vcenter, delta=delta/2)
       point(0, r0, " r0", :black, :left, :top, delta=-delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その403）

算額（その403）

東京都あきる野市二宮　二宮神社 寛政6年(1794)
山口正義：やまぶき，第15号
https://yamabukiwasan.sakura.ne.jp/ymbk15.html

鈎股弦の中に内接する全円と 5 個の等円がある。全円の直径が 44寸，弦が 110 寸のとき，等円の直径はいかほどか。

全円の半径と中心座標を r1, (r1, r1)
等円の半径を r2 として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鈎::positive,  股::positive, 弦::positive, r1::positive,
   r2::positive, x::positive, y::positive;

eq1 = 鈎 + 股 - 弦 ~ 2r1
eq2 = 鈎^2 + 股^2 ~ 弦^2
eq3 = 弦 ~ 股 - x + 鈎 - y + 8r2
eq4 = r1/(鈎 - r1) ~ r2/(鈎 - y)
eq5 = r1/(股 - r1) ~ r2/(股 - x)
res = solve((eq1, eq2, eq3, eq4, eq5), (鈎, 股, r2, x, y))

   2-element Vector{NTuple{5, Sym}}:
    (r1 + 弦/2 - sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1 + 弦/2 + sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1*弦/(8*r1 + 弦), r1*(8*r1 + 5*弦 + 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦), r1*(8*r1 + 5*弦 - 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦))
    (r1 + 弦/2 + sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1 + 弦/2 - sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1*弦/(8*r1 + 弦), r1*(8*r1 + 5*弦 - 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦), r1*(8*r1 + 5*弦 + 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦))

最初の解が適解である。

等円の直径は 2r1*弦/(8*r1 + 弦) である。

r1 = 44//2, 弦 = 110 のとき，等円の直径は 16.923077 である。

   鈎 = 66;  股 = 88;  r2 = 8.46154;  x = 62.6154;  y = 49.0769
   等円の直径 = 16.923077

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, 弦) = (44//2, 110)
   (鈎, 股, r2, x, y) = (r1 + 弦/2 - sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1 + 弦/2 + sqrt(-4*r1^2 - 4*r1*弦 + 弦^2)/2, r1*弦/(8*r1 + 弦), r1*(8*r1 + 5*弦 + 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦), r1*(8*r1 + 5*弦 - 4*sqrt(-4*r1^2 - 4*r1*弦 + 弦^2))/(8*r1 + 弦))
   @printf("鈎 = %g;  股 = %g;  r2 = %g;  x = %g;  y = %g\n",鈎, 股, r2, x, y)
   @printf("等円の直径 = %.8g\n", 2r2)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   #circle(x, r2, r2, :blue)
   #circle(r2, y, r2, :blue)
   for i = 0:4
       circle(x - i*2r2*股/弦, r2 + i*2r2*鈎/弦, r2, :blue)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, r1, "全円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x, r2, "(x,r2)", :blue)
       point(r2, y, "(r2,y)", :blue)
       point(0, 鈎, "鈎 ", :black, :right, :vcenter)
       point(股, 0, "股", :black, :left, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       plot!(xlims=(-5, 90))
   else
       plot!(showaxis=false)
   end
end;

 

算額（その402）

算額（その402）

愛知県岡崎市 六所神社(1779)
http://www.st.nanzan-u.ac.jp/info/ma-thesis/2015/ss/m14ss006.pdf

甲，乙，丙の正方形がある。その面積の和は 61 歩，甲乙の一辺の差は 2 寸，乙丙の一辺の差は 1 寸であるとき，丙の一辺の長さを求めよ。

甲，乙，丙の一辺の長さをそれぞれ 「甲」，「乙」，「丙」とし，以下の方程式を解く。

using SymPy

@syms 甲::positive, 乙::positive, 丙::positive;

eq1 = 甲^2 + 乙^2 + 丙^2 ~ 61
eq2 = 甲 - 乙 ~ 2
eq3 = 乙 - 丙 ~ 1
res = solve((eq1, eq2, eq3), (甲, 乙, 丙))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (6, 4, 3)

正方形 乙 の一辺の長さは 3 寸である。

甲，乙の一辺の長さはそれぞれ 6 寸, 4 寸である。

算額（その401）

算額（その401）

埼玉県所沢市下新井 熊野神社 明治6年(1873)
山口正義：毛呂周辺の算額
https://yamabukiwasan.sakura.ne.jp/22moroshuuhen.pdf

大円内に正7角形と小円が入っている。大円と小円の直径を与えたときに正7角形の一辺の長さを求めよ。

大円と小円の半径を r0, r1，正7角径が内接する円の半径を r2 として，以下の方程式を解く。

当時 sin(π/7) とか sin(π/7) どうやって計算していたのだろうか。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive,
     x2::positive, y2::positive,
     x3::positive, y3::positive,
     x4::positive, y4::positive;

eq = 2r1 + r2*(1 + cos(PI/7))- 2r0
res = solve(eq, r2)[1]
res |> println

   (2*r0 - 2*r1)/(cos(π/7) + 1)

半径は r2 = 2(r0 - r1)/(cos(π/7) + 1) である。
一辺の長さは 2r2*sin(π/7) = 5.477843385363599

r0 = 10; r1 = 4
r2 = 2(r0 - r1)/(cos(π/7) + 1)
2r2*sin(π/7)

   5.477843385363599

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (10, 4)
   r2 = 2(r0 - r1)/(cos(π/7) + 1)
   x = []
   y = []
   for (i, θ) = enumerate(π/2:2π/7:5π/2)
       append!(x, r2*cos(θ))
       append!(y, r2*sin(θ) + r0 - r2)
   end
   @printf("r2 = %g;  r0 = %g;  r1 = %g\n", r2, r0, r1)
   @printf("正7角径の一辺の長さは %g\n", 2r2*sin(π/7))
   plot(x, y, color=:blue, lw=0.5)
   circle(0, 0, r0, :black)
   circle(0, r1 - r0, r1)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(0, r0 - r2, " r0 - r2", :blue)
       point(0, r0-r2*(1+cos(π/7)), " r0-r2*(1+cos(π/7))", :blue, :left, :bottom, delta=delta/2)
       point(0, r1 - r0, " r1-r0", :red)
       point(0, -2, "")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その400）

算額（その400）

埼玉県川島町下小見野 光西寺観音堂 明治25年(1892)
https://yamabukiwasan.sakura.ne.jp/kousaiji.pdf

鈎股弦の中に菱形2個，甲円，乙円，丙円，全円を入れる。甲円，乙円，丙円の直径がそれぞれ 12.5, 7.5, 12 寸のとき，全円の直径はいかほどか。

鈎の上に ya，股の上に xa, xb をとる。その他 yb，yc を定める。
甲円の半径と中心座標を r1, (xb + r1, r1)
乙円の半径と中心座標を r2, (xa + r2, r2)
丙円の半径と中心座標を r3, (r3, r3)

以下の連立方程式を立て解を求める。

include("julia-source.txt");

using SymPy

@syms xa::positive, xb::positive, ya::positive, yb::positive, yc::positive,
     鈎::positive, 股::positive, 弦::positive,
     r0::positive, r1::positive, r2::positive, r3::positive;

(r1, r2, r3) = (125, 75, 120) .// 20
eq1 = ya/xa - 鈎/股
eq2 = yb/(xb - xa) - 鈎/股
eq3 = yc/(股 -xb) - 鈎/股
eq4 = ya + xa - sqrt(ya^2 + xa^2) - 2r3
eq5 = yb + (xb - xa) - sqrt(yb^2 + (xb - xa)^2) - 2r2
eq6 = yc + (股 - xb) - sqrt(yc^2 + (股 - xb)^2) - 2r1
eq7 = ya + yb + yc - 鈎
eq8 = 鈎 + 股 - 弦 - 2r0
eq9 = 鈎^2 + 股^2 - 弦^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

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)
   (xa, xb, ya, yb, yc, 鈎, 股, 弦, r0) = u
   return [
       -鈎/股 + ya/xa,  # eq1
       yb/(-xa + xb) - 鈎/股,  # eq2
       yc/(-xb + 股) - 鈎/股,  # eq3
       xa + ya - sqrt(xa^2 + ya^2) - 12,  # eq4
       -xa + xb + yb - sqrt(yb^2 + (-xa + xb)^2) - 15/2,  # eq5
       -xb + yc + 股 - sqrt(yc^2 + (-xb + 股)^2) - 25/2,  # eq6
       ya + yb + yc - 鈎,  # eq7
       -2*r0 - 弦 + 股 + 鈎,  # eq8
       -弦^2 + 股^2 + 鈎^2,  # eq9
   ]
end;

iniv = [big"48.7", 79.1, 14.0, 8.7, 14.5, 37.2, 129.9, 135.1, 16.0]
res = nls(H, ini=iniv);
println(res);

  (BigFloat[50.14137786721619150424442148054351701570094903819146137553978665178239848983374, 81.47973903422631119439718490588321515051404218706112472703481843278690018275391, 13.88771365970725570233435600656382262425686180304594397465088775657750893885396, 8.679821037317034813958972504102389140160538626903714984058309678348674823548669, 14.46636839552839135659828750683731523360089771150619163995304430858705258296429, 37.03390309255268187289161601750352699801829814145585059866224174351323634536706, 133.7103409792431773446517906147827120418691974351772302874933122640360306862704, 138.7442440717958592175434066322862390398874955766330808838267099554706218187839, 16.00000000000000000000000000000000000000000000000000000116442202603932260642599], true)

   xa = 50.1414;  xb = 81.4797;  ya = 13.8877;  yb = 8.67982;  yc = 14.4664
   鈎 = 37.0339;  股 = 133.71;  弦 = 138.744;  r0 = 16
   全円の直径 = 32;  甲円の直径 = 12.5;  乙円の直径 = 7.5;  丙円の直径 = 12

全円の半径は 16 寸（直径 32 寸）である。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (125, 75, 120) .// 20
   (xa, xb, ya, yb, yc, 鈎, 股, 弦, r0) = res[1]
   @printf("xa = %g;  xb = %g;  ya = %g;  yb = %g;  yc = %g\n鈎 = %g;  股 = %g;  弦 = %g;  r0 = %g\n", xa, xb, ya, yb, yc, 鈎, 股, 弦, r0)
   @printf("全円の直径 = %g;  甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g\n", 2r0, 2r1, 2r2, 2r3)
   x1 = xb + r1
   x2 = xa + r2
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(r3, r3, r3, :green)
   circle(r0, r0, r0, :orange)
   segment(0, ya, xa, 0)
   segment(xa, yb, xb, 0)
   segment(xa, 0, xa, yb + yc)
   segment(xb, 0, xb, yc)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, 鈎, " 鈎=ya+yb+yc", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom, delta=2delta)
       point(xa, 0, "xa", :black, :center, :top, delta=-2delta)
       point(xb, 0, "xb", :black, :center, :top, delta=-2delta)
       point(0, ya, "ya ", :black, :right, :vcenter)
       point(xa, yb, "(xa,yb)", :black, :right, :vcenter)
       point(xa, yb + yc, "(xa,yb+yc)", :black, :left, :bottom, delta=2delta)
       point(xb, yc, "(xb,yc)", :black, :left, :bottom, delta=2delta)
       point(xb + r1, r1, " 甲円:r1,(xb+r1,r1)", :red, :left, :bottom)
       point(xa + r2, r2, " 乙円:r2,(xa+r2,r2)", :blue, :left, :top, delta=2delta)
       point(r3, r3, " 丙円:r3,(r3,r3)", :green, :left, :top, delta=2delta)
       point(r0, r0, " 全円:r0,(r0,r0)", :orange, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-8, 140), ylims=(-6, 42))
   else
       plot!(showaxis=false)
   end
end;

算額（その399）

算額（その399）

埼玉県川島町下小見野 光西寺観音堂 明治25年(1892)
https://yamabukiwasan.sakura.ne.jp/kousaiji.pdf

外円の中に正方形および大円と小円が入っている。外円の直径が16寸，大円の直径が 8 寸のとき小円の直径はいかほどか。

外円の半径と中心座標を r0, (0, r0)
大円の半径と中心座標を r1, (0, 0)
右上の小円の半径と中心座標を r2, (x2, x2)

以下の連立方程式を立て解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive;

eq1 = (r0 - sqrt(Sym(2))r1)/2 - r2
eq2 = (sqrt(Sym(2))r1 + r2/2)/sqrt(Sym(2)) - x2;
res = solve([eq1, eq2], (r2, x2))

   Dict{Any, Any} with 2 entries:
     r2 => r0/2 - sqrt(2)*r1/2
     x2 => sqrt(2)*r0/8 + 3*r1/4

   r2 = 1.17157;  x2 = 4.82843
   小円の直径 = 2.34315

小円の半径は，「外円の半径/2 - 大円の半径√2/2」である。
外円の直径が 8 寸，大円の直径が 4 寸のとき，小円の直径は，8 - 4√2　= 2.3431457505076194...である。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (16, 8) .// 2
   (r2, x2) = (r0/2 - sqrt(2)*r1/2, sqrt(2)*r0/4 + r1/2)
   @printf("r2 = %g;  x2 = %g\n", r2, x2)
   @printf("小円の直径 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :black)
   rect(-r1, -r1, r1, r1, :blue)
   circle(0, 0, r1)
   circle4(x2, x2, r2, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r0, 0, "r0 ", :black, :right, :bottom, delta=delta)
       point(r1, 0, " r1", :blue, :left, :bottom, delta=delta)
       point(x2, x2, "小円:r2,(x2,y2) ", :green, :right, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その398）

算額（その398）

東京都府中市 大國魂神社 明治18年(1885)
http://www.f.waseda.jp/takezawa/sousuukyo/paper/2011ronbun.pdf

半円内に甲円，乙円，丙円，丁円，戊円の 9 円をいれる。外円の直径が 10 寸のとき，丙円の直径を求めよ。

甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5, y5)
以下の連立方程式を立て，nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive, r5::positive, x5::positive, y5::positive;

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = x5^2 + (y5 - r1)^2 - (r1 + r5)^2
eq4 = (x4 - x2)^2 + (y4 - r2)^2 - (r2 + r4)^2
eq5 = (x2 - x5)^2 + (y5 - r2)^2 - (r2 + r5)^2
eq6 = (x4 - x3)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq7 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = x2^2 + r2^2 - (2r1 - r2)^2
eq10 = x3^2 + y3^2 - (2r1 - r3)^2
eq11 = x4^2 + y4^2 - (2r1 - r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r1, r3, x3, r4, x41, r42, x42, y42))

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-14")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-14")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r2, x2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
       x2^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       x3^2 + (-r1 + y3)^2 - (r1 + r3)^2,  # eq2
       x5^2 + (-r1 + y5)^2 - (r1 + r5)^2,  # eq3
       (-r2 + y4)^2 - (r2 + r4)^2 + (-x2 + x4)^2,  # eq4
       (-r2 + y5)^2 - (r2 + r5)^2 + (x2 - x5)^2,  # eq5
       -(r3 + r4)^2 + (-x3 + x4)^2 + (y3 - y4)^2,  # eq6
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq7
       -(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2,  # eq8
       r2^2 + x2^2 - (2*r1 - r2)^2,  # eq9
       x3^2 + y3^2 - (2*r1 - r3)^2,  # eq10
       x4^2 + y4^2 - (2*r1 - r4)^2,  # eq11
   ]
end;

r1 = 5//2
iniv = [big"1.3", 3.45, 0.52, 2.76, 3.36, 0.34, 3.62, 2.93, 0.26, 2.84, 2.5]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.25, 3.535533905932737622004221810526050308609409471587999437373189411634045676321378, 0.5000000000000000000914962217596418053047474197754641655637114736626653004349011, 2.828427124746190097802451876811358797663566934419370521520740609153604311115202, 3.49999999999999999972551133472107458408575774067360750330886560977611543609143, 0.4166666666666666664757008581209034940070734535154511572023238905891655189934959, 3.535533905932737622365050944667911778662710884062463557643917323522673436115187, 2.916666666666666666600917687737661690052746729593517424887121308994452716368222, 0.3333333333333333328134030564312976635069374013893610682603843011810001642030419, 2.828427124746190097243674263049636031780028552049905775368715116058220501682123, 2.666666666666666666169204262275820633241336297105830855213052541023599670593682], true)

r2 = 1.25;  x2 = 3.53553;  r3 = 0.5;  x3 = 2.82843;  y3 = 3.5;  r4 = 0.416667;  x4 = 3.53553;  y4 = 2.91667;  r5 = 0.333333;  x5 = 2.82843;  y5 = 2.66667
丙円の直径 = 1

外円の直径が 10 寸のとき，丙円の直径は 1 寸である。

using Plots

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

算額（その397）

算額（その397）

長野県北佐久郡軽井沢町峠 熊野神社 明治5年(1872)
山口正義：「やまぶき3」
https://yamabukiwasan.sakura.ne.jp/ymbk43.pdf

外円内に，甲円，乙円，丙円，丁円が入っている。乙円の直径がわかっているとき，甲円の直径を求めよ。

算額の図では 丁円が 8 個あるとされているが，そのような仮定のもとでは未知数が 7 個であるにもかかわらず条件数が 8 個となる。実際に解を求める場合に矛盾が生じる（図が描けない）。甲円と接する円は x 軸に接する丁円とは半径が異なるとすれば，未知数と条件数が等しくなる。そこで，前者の円を戊円とする。
甲円の半径と中心座標を r1, (0, r2 + r1)
乙円の半径と中心座標を r2, (0, 0)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r41, (x41, y41)
戊円の半径と中心座標を r42, (x42, y42)
以下の連立方程式を立て，nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, x3::positive,
     r41::positive, x41::positive, r42::positive, x42::positive, y42::positive;
@syms r1, r2, r3, x3, r4, x41, r42, x42, y42

r0 = 2r1 + r2
eq1 = x3^2 + (r2 + r1 - r3)^2 - (r1 + r3)^2
eq2 = x42^2 + (r2 + r1 - y42)^2 - (r1 + r42)^2
eq3 = x41^2 + r41^2 - (r2 + r41)^2
eq4 = x42^2 + y42^2 - (r2 + r42)^2
eq5 = (x3 - x41)^2 + (r3 - r41)^2 - (r3 + r41)^2
eq6 = (x3 - x42)^2 + (r3 - y42)^2 - (r3 + r42)^2
eq7 = (x41 - x42)^2 + (y42 - r41)^2 - (r41 + r42)^2
eq8 = x3^2 + r3^2 - (r0 - r3)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r1, r3, x3, r4, x41, r42, x42, y42))

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r1, r3, x3, r41, x41, r42, x42, y42) = u
   return [
       x3^2 - (r1 + r3)^2 + (r1 + r2 - r3)^2,  # eq1
       x42^2 - (r1 + r42)^2 + (r1 + r2 - y42)^2,  # eq2
       r41^2 + x41^2 - (r2 + r41)^2,  # eq3
       x42^2 + y42^2 - (r2 + r42)^2,  # eq4
       (r3 - r41)^2 - (r3 + r41)^2 + (x3 - x41)^2,  # eq5
       -(r3 + r42)^2 + (r3 - y42)^2 + (x3 - x42)^2,  # eq6
       (-r41 + y42)^2 - (r41 + r42)^2 + (x41 - x42)^2,  # eq7
       r3^2 + x3^2 - (2*r1 + r2 - r3)^2,  # eq8
   ]
end;

r2 = 1
iniv = [big"19.0", 11, 30, 4, 17, 4, 13, 11] ./ 14
res = nls(H, ini=iniv);
println(res);
   (BigFloat[1.590635214223556758596745603528478178766990154309840919228137335063082400883928, 1.295317607111778379298372825977315420000085315906557774300910306055315016133428, 2.578929232003002706842677554771274123131823604105670170499680855552520301272628, 0.323829401777944594824592995962549474245268966811601990866971647246166910750795, 1.283611624891224327544304513220693406771330824235714117692914677600658970155483, 0.3457573000500555096803689076368287502783069832673107490451059925743197817098452, 0.9810739873474257565739818948867657209770995516904925318028651309533017094462853, 0.9211712880828615044789766470750869333771167444223464888129631384631261089220609], true)

   r1 = 1.59064;  r3 = 1.29532;  x3 = 2.57893;  r4 = 0.323829;  x41 = 1.28361;  r42 = 0.345757;  x42 = 0.981074;  y42 = 0.921171

乙円の半径が 1 のとき 甲円の半径は 1.5906352142235

術では，甲円の直径は乙円の直径の (sqrt(1152)+44)/49 = 1.5906352142235 倍である。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #r2 = 1
   (r1, r3, x3, r41, x41, r42, x42, y42) = res[1]
   @printf("r1 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x41 = %g;  r42 = %g;  x42 = %g;  y42 = %g\n", r1, r3, x3, r41, x41, r42, x42, y42)
   plot()
   circle(0, 0, r2 + 2r1, :black)
   circle(0, r2 + r1, r1)
   circle(0, -r2 - r1, r1)
   circle(0, 0, r2, :blue)
   circle4(x3, r3, r3, :orange)
   circle4(x41, r41, r41, :green)
   circle4(x42, y42, r42, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r2+r1, " 甲円\n r2＋r1", :red, :left, :vcenter)
       point(0, 0, " 乙円:r2\n (0,0)", :blue, :left, :vcenter)
       point(x3, r3, " 丙円:r3\n (x3,y3)", :orange, :left, :vcenter)
       point(x41, r41, " 丁円:r41,(x41,r41)", :green, :left, :vcenter)
       point(x42, y42, " 戊円:r42,(x42,r42)", :magenta, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その396）

算額（その396）

東京都渋谷区 金王神社（金王八幡宮） 元治元年(1864)
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf

金王八幡宮③「宝物館」（渋谷散歩③）
https://wheatbaku.exblog.jp/30049403/

直線上に大円，中円，小円が載っており，互いに接している。中円，小円の直径がそれぞれ 9 寸，4 寸のとき，大円の直径はいくつか。

大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (x3, r3)
として，以下の連立方程式を立て解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, x1::positive, r2::positive,
     r3::positive, x3::positive;

eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, x1, x3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-2*r2^(11/2)*r3^(3/2) + 8*r2^(9/2)*r3^(5/2) - 12*r2^(7/2)*r3^(7/2) + 8*r2^(5/2)*r3^(9/2) - 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2, 2*r2^(3/2)*sqrt(r3)/(r2 - r3) - 2*r2*r3/(r2 - r3), 2*sqrt(r2)*sqrt(r3))
    ((2*r2^(11/2)*r3^(3/2) - 8*r2^(9/2)*r3^(5/2) + 12*r2^(7/2)*r3^(7/2) - 8*r2^(5/2)*r3^(9/2) + 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2, 2*r2^(3/2)*sqrt(r3)/(r2 - r3) + 2*r2*r3/(r2 - r3), 2*sqrt(r2)*sqrt(r3))

2 番目の解が適解である。

   r1 = 18;  x1 = 18;  x3 = 6
   大円の直径 = 36

r1 を求めるために得られた式は長く複雑である。

(r2, r3) = (9, 4) .// 2
(2*r2^(11/2)*r3^(3/2) - 8*r2^(9/2)*r3^(5/2) + 12*r2^(7/2)*r3^(7/2) - 8*r2^(5/2)*r3^(9/2) + 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2

   18.000000000000014

術では，中円径 ÷ (sqrt(中円径/小円径) - 1)^2 と簡単である。

(中円径, 小円径) = (9, 4)
中円径 / (sqrt(中円径/小円径) - 1)^2

   36.0

この式を得るには，直線上に載っていて外接する 2 円の中心の水平距離に関する公式 d12 = 2sqrt(r1*r2) を用いて以下のように展開する。

r1, r3 から r2 を求める場合も，r1, r2 から r3 を求める場合も同じである。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (9, 4) .// 2
   (r1, x1, x3) = ((2*r2^(11/2)*r3^(3/2) - 8*r2^(9/2)*r3^(5/2) + 12*r2^(7/2)*r3^(7/2) - 8*r2^(5/2)*r3^(9/2) + 2*r2^(3/2)*r3^(11/2) + r2^6*r3 - 3*r2^5*r3^2 + 2*r2^4*r3^3 + 2*r2^3*r3^4 - 3*r2^2*r3^5 + r2*r3^6)/(r2^3 - 3*r2^2*r3 + 3*r2*r3^2 - r3^3)^2, 2*r2^(3/2)*sqrt(r3)/(r2 - r3) + 2*r2*r3/(r2 - r3), 2*sqrt(r2)*sqrt(r3))
   @printf("r1 = %g;  x1 = %g;  x3 = %g\n", r1, x1, x3)
   @printf("大円の直径 = %g\n", 2r1)
   plot()
   circle(x1, r1, r1)
   circle(0, r2, r2, :blue)
   circle(x3, r3, r3, :orange)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r1, " 大:r1,(x1,r1)", :red, :center, :bottom, delta=delta)
       point(0, r2, " 中:r2,(0,r2)", :blue, :center, :bottom, delta=delta)
       point(x3, r3, "小:r3,(x3,r3) ", :orange, :right, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その395）

算額（その395）

東京都府中市 大國魂神社 明治18年(1885)
佐藤健一『多摩の算額』
山口正義『やまぶき』第8号
https://yamabukiwasan.sakura.ne.jp/ymbk8.pdf

半円の中に大円，小円，甲円，乙円，丙円が入っている。甲円，乙円の直径がそれぞれ 36 寸，22 寸のとき，丙円の直径を求めよ。

半円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (x2, r2)
甲円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (x4, r4)
丙円の半径と中心座標を r5, (x5, y5)
として，以下の連立方程式を立て，nlsove() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, x1::positive, r2::positive,
     x2::negative, r3::positive, x3::positive, r4::positive,
     x4::negative, r5::positive, x5::negative;

(r3, r4) = (Sym(36), Sym(22)) .// 2
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2
eq4 = (x4 - x2)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq5 = (x2 - x5)^2 + (r2 - r5)^2 - (r2 + r5)^2
eq6 = x1^2 + r1^2 - (r0 - r1)^2
eq7 = x2^2 + r2^2 - (r0 - r2)^2
eq8 = x3^2 + r3^2 - (r0 - r3)^2
eq9 = x5^2 + r5^2 - (r0 - r5)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

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-63")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-63")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r0, r1, x1, r2, x2, x3, x4, r5, x5) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq1
       (r1 - 18)^2 - (r1 + 18)^2 + (-x1 + x3)^2,  # eq2
       (r1 - 11)^2 - (r1 + 11)^2 + (x1 - x4)^2,  # eq3
       (r2 - 11)^2 - (r2 + 11)^2 + (-x2 + x4)^2,  # eq4
       (r2 - r5)^2 - (r2 + r5)^2 + (x2 - x5)^2,  # eq5
       r1^2 + x1^2 - (r0 - r1)^2,  # eq6
       r2^2 + x2^2 - (r0 - r2)^2,  # eq7
       x3^2 - (r0 - 18)^2 + 324,  # eq8
       r5^2 + x5^2 - (r0 - r5)^2,  # eq9
   ]
end;

iniv = [big"76.0", 34.0, 24.0, 32.0, -36.0, 60.0, -4.0, 12.0, -70.0]
res = nls(H, ini=iniv);
println(res);

  (BigFloat[109.2759607781864242781305327440617590523959228887413843855648291428112275807799, 50.78632222480445137700303498597358678562931717751496146783878544338869690961004, 29.0135708416527090497508430387449613589229832428047177392913207362097294847525, 38.48830848061307241285693076285644633867531725256882355802544542966271346448629, -59.40994721511970583145682953216998977779327580589532321869496230814697858668088, 89.48352371236297305829391802881040814933023875202264524826791200242400134239185, -18.25796581510720064785824833255286994844925646259226400277467085047377787714211, 10.0405052747635702369533443385048324890519598382837919173533224253202548503772, -98.72620666671672609506542347825094825128246061329118501178770677156669483447183], true)

r0 = 109.276;  r1 = 50.7863;  x1 = 29.0136;  r2 = 38.4883;  x2 = -59.4099;  x3 = 89.4835;  x4 = -18.258;  r5 = 10.0405;  x5 = -98.7262
丙円の直径 = 20.081

術を読み解くと以下のようである。
『8の平方根から1を引き天とする。乙円径を甲円径で割ったものの平方根をとり天を引き二乗したもので，乙円径を割る』

(乙円径, 甲円径) = (22, 36)
天 = √8 - 1
乙円径/(√(乙円径/甲円径) - 天)^2

   20.08101054952714

確かに数値解と同じになる。

2res[1][8]

   20.08101054952714047390668867700966497810391967656758383470664485064050970075439

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, x1, r2, x2, x3, x4, r5, x5) = res[1]
   (r3, r4) = (36, 22) .// 2
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  x3 = %g;  x4 = %g;  r5 = %g;  x5 = %g\n", r0, r1, x1, r2, x2, x3, x4, r5, x5)
   @printf("丙円の直径 = %g\n", 2r5)
   plot()
   circle(0, 0, r0, :black, beginangle=0, endangle=180)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :orange)
   circle(x4, r4, r4, :green)
   circle(x5, r5, r5, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r1, " 大:r1,(x1,r1)", :red, :center, :bottom, delta=delta)
       point(x2, r2, " 小:r2,(x2,r2)", :blue, :center, :bottom, delta=delta)
       point(x3, r3, "甲:r3,(x3,r3) ", :orange, :right, :vcenter)
       point(x4, r4, " 乙:r4,(x4,r4)", :green, :left, :vcenter)
       point(x5, r5, " 丙:r5,(x5,r5)", :magenta, :left, :vcenter)
       point(0, r0, " r0", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;