算額（その441）

算額（その441）

東都下谷稲荷社 寛政元年己酉3月(1624)
黒川喜兵衛盛榮

藤田貞資(1789)：神壁算法
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

円弧の中に甲円 1 個，乙円 2 個，丙円 2 個があり，それぞれの円は隣同士は概説し，円弧を構成する外円には内接し，弦に接している。
甲円，丙円の直径がそれぞれ 12 寸，3 寸のとき，乙円の直径はいかほどか。

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

外円の中心を原点とする。
外円の半径，中心座標を r0, (0, 0)
甲円の半径，中心座標を r1, (0, a + r1)
乙円の半径，中心座標を r2, (x2, a + r2)
丙円の半径，中心座標を r3, (x3, a + r3)
弦が y 軸と交わる座標を (0, a)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, a::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive;

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq3 = x2^2 + (a + r2)^2 - (r0 -r2)^2
eq4 = x3^2 + (a + r3)^2 - (r0 -r3)^2
eq5 = a + 2r1 - r0

res = solve([eq1, eq2, eq3, eq4, eq5], (a, r0, r2, x2, x3))

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

2 番めのものが適解である。

names = ("a", "r0", "r2", "x2", "x3")
for (i, name) in enumerate(names)
   @printf("%s = %g;  ", name, res[2][i](r1 => 6, r3 => 1.5).evalf())
end

   a = 6;  r0 = 18;  r2 = 4;  x2 = 9.79796;  x3 = 14.6969;  

甲円，丙円の直径が 12 寸，3 寸のとき，乙円の直径は 8 寸である。

using Plots

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

算額（その440）

算額（その440）

東都下谷摩利支天堂 天明8年戊申3月(1788)
東都浅草新堀端住 三浦伴次郎成喜

藤田貞資(1789)：神壁算法
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

求めるものが小頭か大頭かの違いだけで，算額（その185）と同じ問題である。プログラミングスタイルを最近のものに変えた。

等脚台形の隔斜（対角線）で区切られた領域に甲，乙，丙の円を入れる。甲円の直径が 100 寸，乙円の直径が 28 寸，丙円の直径が 45 寸のとき，大頭（台形の下底）はいかほどか。

以下のように記号を定め方程式を解く。
台形の右下隅 A，右上隅 b， の座標をそれぞれ (a, 0)，(b, y) と置く。
甲円，乙円，丙円の半径をそれぞれ，r1，r2，r3 とする。
甲円，乙円，丙円の中心から隔斜（対角線）までの距離が円の半径に等しいという連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, x1::positive, y1::positive, y::positive,
     r1::positive, r2::positive, r3::positive;

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

(r1, r2, r3) =(100, 28, 45) ./ 2
iniv = [big"151.0", 44, 39, 112, 144]
res = nls(H, ini=iniv)
names = ("a", "b", "x1", "y1", "y")
if res[2]
   for (name, x) in zip(names, res[1])
       @printf("%s = %g;  ", name, round(Float64(x), digits=6))
   end
   println()
else
   println("収束していない")
end

   a = 150;  b = 42;  x1 = 37.5;  y1 = 112.5;  y = 144;  

大頭の長さは 2a = 300（寸）である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) =(100, 28, 45) ./ 2
   (a, b, x1, y1, y) = res[1]
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   @printf("a = %g; b = %g;  x1 = %g; y1 = %g; y = %g\n", a, b, x1, y1, y)
   plot([a, b, -b, -a, a], [0, y, y, 0, 0], color=:orange, lw=0.5)
   circle(0, r1, r1, :green)
   circle(x1, y1, r3, :blue)
   circle(-x1, y1, r3, :blue)
   circle(0, y - r2, r2, :red)
   segment(-b, y, a, 0, :orange)
   segment(-a, 0, b, y, :orange)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1, " 甲円:r1,(0,r1)", :black, :left, :vcenter)
       point(0, y - r2, " 乙円:r2,(0,y-r2)", :black, :left, :vcenter)
       point(x1, y1, " 丙円:r3,(x1,y1,r3)", :black, :left, :vcenter)
       point(b, y, " B(b,y)", :black, :left, :bottom, delta=delta/2)
       point(a, 0, " A(a,0)", :black, :left, :bottom, delta=delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その439）

算額（その439）

東都本郷真光寺天満宮 天明7年丁未11月(1787)
東都湯嶋御徒町 池田重次郎教政

藤田貞資(1789)：神壁算法
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円の一部分である円弧内に，大中小の 5 個の円が入っている。大円と小円は弦に接している。
弦の長さが 8 寸，小円の直径が 1 寸のとき，矢（弦と円の距離，図では r0 - a）の長さはいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, a + r1)
中円の半径と中心座標を r2, (0, r0 - r2)
小円の半径と中心座標を r3, (r3, a + r3)
弦と y 軸の交点座標を (0, a)
として以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive,
     r1::positive, x1::positive,
     r2::positive, r3::positive,
     弦::positive, a::positive;

eq1 = x1^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = (x1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = x1^2 + (r0 - r2 - a - r1)^2 - (r1 + r2)^2
eq4 = r3^2 + (r0 - r2 - a - r3)^2 - (r2 + r3)^2 
eq5 = (r0^2 - (弦/2)^2) - a^2;

res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, x1, r2, a))

   2-element Vector{NTuple{5, Sym}}:
    ((1536*r3^6 + 288*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 + 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), r3*(128*r3^4 - 8*r3^2*弦^2 - 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 + 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)), r3*弦*(2*弦 + sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2), r3*(96*r3^4 + 44*r3^2*弦^2 + 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 + 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), (-1536*r3^6 + 1440*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 + 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))
    ((1536*r3^6 + 288*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 - 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 + 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), r3*(128*r3^4 - 8*r3^2*弦^2 + 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 - 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)), r3*弦*(2*弦 - sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2), r3*(96*r3^4 + 44*r3^2*弦^2 - 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 - 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), (-1536*r3^6 + 1440*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 - 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 + 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))

最初の組が適解である。

矢の長さは 2*r3*(3*r3^2 - 4*sqrt(r3^2 + 20) - 24)/(9*r3^2 - 16)

res[1][1] - res[1][5] |> simplify |> println

   r3*(24*r3^2 - 3*弦^2 - 弦*sqrt(16*r3^2 + 5*弦^2))/(36*r3^2 - 弦^2)

弦の長さが 8 寸，小円の直径が 1 寸 のとき，矢の長さは 3 寸である。

(r3, 弦) = (1/2, 8)
r3*(24*r3^2 - 3*弦^2 - 弦*sqrt(16*r3^2 + 5*弦^2))/(36*r3^2 - 弦^2)

   3.0

   r0 = 4.16667;  r1 = 1.125;  x1 = 2;  r2 = 1.04167;  a = 1.16667
   矢 = 3

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, 弦) = (1/2, 8)
   (r0, r1, x1, r2, a) = ((1536*r3^6 + 288*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 + 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)),
       r3*(128*r3^4 - 8*r3^2*弦^2 - 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 + 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)),
       r3*弦*(2*弦 + sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2),
       r3*(96*r3^4 + 44*r3^2*弦^2 + 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 + 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)),
       (-1536*r3^6 + 1440*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 + 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  a = %g\n",
           r0, r1, x1, r2, a)
   @printf("矢 = %g\n", r0 - a)
   plot()
   circle(0, 0, r0, :gray)
   circle(x1, a + r1, r1)
   circle(-x1, a + r1, r1)
   circle(0, r0 - r2, r2, :blue)
   circle(r3, a + r3, r3, :orange)
   circle(-r3, a + r3, r3, :orange)
   b = sqrt(r0^2 - a^2)
   segment(-b, a, b, a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, a + r1, "大円:r1,(x1,a+r1)", :red, :center, :bottom, delta=delta/2)
       point(0, r0 - r2, "中円:r2,(0,r0-r2)", :blue, :center, :bottom, delta=delta/2)
       point(r3, a + r3, " 小円:r3,(r3,a+r3)", :black, :left, :vcenter)
       point(0, a, " a", :black, :left, :top, delta=-delta/2)
       point(弦/2, a, "(弦/2,a)", :black, :right, :top, delta=-delta/2)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(r0, 0, "r0 ", :black, :right, :bottom, delta=delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5, ylims=(-0.3, 4.3))
   else
       plot!(showaxis=false)
   end
end;

算額（その438）

算額（その438）

社丈右衛門正常 天明2年壬寅4月(1782)

藤田貞資(1789)：神壁算法 上巻
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

交差する 2 個の大円の中に，等円 3 個，甲円，乙円，丙円がそれぞれ 4 こずつ入っている。大円の直径が与えられたとき，大円の直径が与えられたとき，各円の直径を求めよ。

大円の半径と中心座標を 2r, (0, r), (0, -r)
等円の半径と中心座標を r, (0, 2r), (0, 0), (0, -2r)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
として以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive,
     r5::positive, x5::positive, y5::positive,
     r6::positive, x6::positive, y6::positive;

r = r0//2
eq1 = x1^2 + (y1 + r)^2 - (2r + r1)^2
eq2 = x1^2 + (y1 - r)^2 - (2r - r1)^2
eq3 = x1^2 + (2r - y1)^2 - (r1 + r)^2

eq4 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = x2^2 + (r - y2)^2 - (2r - r2)^2
eq6 = x2^2 + (y2 + r)^2 - (2r + r2)^2

eq7 = x3^2 + (r - y3)^2 - (2r - r3)^2
eq8 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq9 = x3^2 + (y3 + r)^2 - (2r + r3)^2;

solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (r1, x1, y1, r2, x2, y2, r3, x3, y3))

   2-element Vector{NTuple{9, Sym}}:
    (3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 27*r0/730, 182*sqrt(3)*r0/365, 27*r0/365)
    (3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5)

最初の組が適解である。
(3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 27*r0/730, 182*sqrt(3)*r0/365, 27*r0/365)

r0 = 36; r = 18;  r1 = 10.8; x1 = 24.9415; y1 = 21.6; r2 = 3.95122; x2 = 30.4165; y2 = 7.90244; r3 = 1.33151; x3 = 31.0915; y3 = 2.66301

円の直径だけをいえば，甲円，乙円，丙円の直径は大円の直径のそれぞれ 3/10, 9/82, 27/730 倍である。

大円の直径 = 72 のとき，等円の直径 = 36;  甲円の直径 = 21.6;  乙円の直径 = 7.90244;  丙円の直径 = 2.66301

72 .* [3/10, 9/82, 27/730] |> println

   [21.599999999999998, 7.902439024390244, 2.663013698630137]

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 72/2
   r = r0/2
   (r1, x1, y1, r2, x2, y2, r3, x3, y3) = (3*r0/10, 2*sqrt(3)*r0/5, 3*r0/5, 9*r0/82, 20*sqrt(3)*r0/41, 9*r0/41, 27*r0/730, 182*sqrt(3)*r0/365, 27*r0/365)
   @printf("r0 = %g; r = %g;  r1 = %g; x1 = %g; y1 = %g; r2 = %g; x2 = %g; y2 = %g; r3 = %g; x3 = %g; y3 = %g\n",
           r0, r, r1, x1, y1, r2, x2, y2, r3, x3, y3)
   @printf("大円の直径 = %g;  等円の直径 = %g;  甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g\n", 2r0, 2r, 2r1, 2r2, 2r3)
   plot()
   circle(0, r, r0, :gray)
   circle(0, -r, r0, :gray)
   circle(0, 2r, r, :blue)
   circle(0, -2r, r, :blue)
   circle(0, 0, r, :blue)
   circle4(x1, y1, r1)
   circle4(x2, y2, r2, :blue)
   circle4(x3, y3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, "  甲円:r1,(x1,y1)", :black, :left, :vcenter)
       point(x2, y2, "   乙円:r2,(x2,y2)", :black, :left, :vcenter)
       point(x3, y3, "  丙円:r3,(x3,y3)", :black, :left, :vcenter)
       point(0, r, " r")
       point(0, -r, " -r")
       point(0, 2r, " 2r")
       point(0, 3r, " 3r")
       point(0.08r, 2.4r, "等円", :blue, mark=false)
       point(0.8x1, 2.1r, "大円", :gray, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その437）

算額（その437）

橘田彌曾八元克 天明八年戊申2月(1788)

藤田貞資(1789)：神壁算法 上巻
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

大円の中に，甲，乙，丙，丁，戊，己の 6 個の円が入っている。丁円，戊円，己円の直径がそれぞれ 872.3 寸，671 寸，572 寸のとき，大円の直径はいかほどか。

（必要ならば）己円の中心がy 軸上に来るように回転する。中心の x 座標が 0 になる。
大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)　r4 = 872.3/2
戊円の半径と中心座標を r5, (x5, y5)  r5 = 671/2
己円の半径と中心座標を r6, (x6, y6)  r6 = 572/2, x6 = 0
として以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r0,
      r1, x1, y1, r2, x2, y2, r3, x3, y3,
      r4, x4, y4, r5, x5, y5, r6, x6, y6

eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = x3^2 + y3^2 - (r0 - r3)^2
eq4 = x4^2 + y4^2 - (r0 - r4)^2
eq5 = x5^2 + y5^2 - (r0 - r5)^2
eq6 = x6^2 + y6^2 - (r0 - r6)^2
eq7 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq8 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq9 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq10 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq11 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq12 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq13 = (x2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq14 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq15 = (x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^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)
   (r0, r1, x1, y1, r2, x2, y2, r3, x3, y3, x4, y4, x5, y5, y6) = u
   return [
       x1^2 + y1^2 - (r0 - r1)^2,  # eq1
       x2^2 + y2^2 - (r0 - r2)^2,  # eq2
       x3^2 + y3^2 - (r0 - r3)^2,  # eq3
       x4^2 + y4^2 - (r0 - r4)^2,  # eq4
       x5^2 + y5^2 - (r0 - r5)^2,  # eq5
       x6^2 + y6^2 - (r0 - r6)^2,  # eq6
       -(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2,  # eq7
       -(r1 + r3)^2 + (x1 - x3)^2 + (y1 - y3)^2,  # eq8
       -(r1 + r4)^2 + (x1 - x4)^2 + (y1 - y4)^2,  # eq9
       -(r1 + r5)^2 + (x1 - x5)^2 + (y1 - y5)^2,  # eq10
       -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq11
       -(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2,  # eq12
       -(r2 + r6)^2 + (x2 - x6)^2 + (y2 - y6)^2,  # eq13
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq14
       -(r3 + r6)^2 + (x3 - x6)^2 + (y3 - y6)^2,  # eq15
   ]
end;

(r4, r5, r6) = (872.3, 671, 572) ./2
x6 = 0
iniv = [big"95.0", 57, 12, 40, 41, -40, -43, 34, 36, -47, -72, 18, 75, -8, -79] .* (620.3/28)
res = nls(H, ini=iniv);
println([round(Float64(x), digits=6) for x in res[1]], " 収束:", res[2]);

   [1586.0, 830.761905, 265.84381, 706.902857, 726.916667, -668.763333, -539.24, 623.071429, 672.917143, -688.777143, -994.422, 577.304, 1248.06, 78.08, -1300.0] 収束:true

r0 = 1586
r1 = 830.762; x1 = 265.844;  y1 = 706.903
r2 = 726.917; x2 = -668.763; y2 = -539.24
r3 = 623.071; x3 = 672.917;  y3 = -688.777
r4 = 436.15;  x4 = -994.422; y4 = 577.304
r5 = 335.5;   x5 = 1248.06;  y5 = 78.08
r6 = 286;     x6 = 0;        y6 = -1300

大円の半径は 1586 寸（直径は 3172 寸）である。

using Plots

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

算額（その436）

算額（その436）

橘田彌曾八元克 天明八年戊申2月(1788)

藤田貞資(1789)：神壁算法 上巻
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円の中に，甲円 3 個，乙円 2 個，丙円 1 個が入っている。甲円は弦に接している。
外円の直径が 3 寸 6 分 のとき，甲円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (2r1, a + r1)
乙円の半径と中心座標を r2, (r1, y2)
丙円の半径と中心座標を r3, (0, r0 - r3)
として以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, a::positive, r1::negative,
     r3::positive, r2::positive, y2::positive;
r0 = 36//20
eq1 = (2r1)^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = r1^2 + (y2 - a - r1)^2 - (r1 + r2)^2
eq3 = r1^2 + y2^2 - (r0 - r2)^2
eq4 = r1^2 + (r0 - r3 - y2)^2  - (r3 + r2)^2
eq5 = 2r3 + 2r1 + a - r0;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, r1, r3, r2, y2))

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, r3, r2, y2) = u
   return [
       4*r1^2 + (a + r1)^2 - (r0 - r1)^2,  # eq1
       r1^2 - (r1 + r2)^2 + (-a - r1 + y2)^2,  # eq2
       r1^2 + y2^2 - (r0 - r2)^2,  # eq3
       r1^2 - (r2 + r3)^2 + (r0 - r3 - y2)^2,  # eq4
       a - r0 + 2*r1 + 2*r3,  # eq5
   ]
end;
r0 = 3.6/2
iniv = r0 .* [big"0.22", 0.24, 0.14, 0.12, 0.69]
res = nls(H, ini=iniv);
println([round(Float64(x), digits=6) for x in res[1]], " 収束:", res[2]);

   [0.330519, 0.500029, 0.234712, 0.282615, 1.43263] 収束:true

a = 0.330519;  r1 = 0.500029;  r3 = 0.234712;  r2 = 0.282615;  y2 = 1.43263

外円の直径が 3.6 のとき，甲円の直径 は 1.00006 である。

using Plots

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

算額（その435）

算額（その435）

東京都府中市 大國魂神社 明治18年(1885)

山口正義(2015)：やまぶき2 第33号
https://yamabukiwasan.sakura.ne.jp/ymbk33.pdf

外円とその平行な 2 本の弦に内接・外接する，甲円 2 個，乙円 1 個，丙円 3 個が入っている。甲円，乙円，丙円，外円も内接・外接している。
外円と甲円の直径の差が 4 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r1, a + r1)
乙円の半径と中心座標を r2, (0, r0 - r2)
丙円の半径と中心座標を r3, (x3, a + 2r1 + r3), (0, a + 2r1 - r3)
甲円と接する下側の弦と y 軸の交点の y 座標を a
外円と甲円の直径の差を A
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms A::positive, r0::positive, a::negative,
     r1::positive, r2::positive, r3::positive, x3::positive;
eq1 = r1^2 + (a + r1)^2 - (r0 -r1)^2
eq2 = r1^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = 2r0 - 2r1 - A
eq5 = x3^2 + (a + 2r1 + r3)^2 - (r0 - r3)^2
eq6 = r0 + a + 2r1 + 2r2 - 2r0
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r0, a, r1, r2, r3, x3))

   1-element Vector{NTuple{6, Sym}}:
    (A*(sqrt(3) + 2)/4, A*(1 - sqrt(3))/4, sqrt(3)*A/4, A/8, sqrt(3)*A/16, sqrt(2)*3^(1/4)*A/8)

   r0 = 3.73205;  a = -0.732051;  r1 = 1.73205;  r2 = 0.5;  r3 = 0.433013;  x3 = 0.930605

乙円の半径は A/8。すなわち，外円と甲円の直径の差の 1/8 である。乙円の直径は外円と甲円の直径の差の 1/4 である。
つまり，外円と甲円の直径の差が 4 寸のとき，乙円の直径は 1 寸である。

山口は「大國魂神社奉献 関流和算額」という小冊子に書かれている本問への答えに失望したと書いてるが，ここで示したのが正解だ（答と術は正しい）と伝えたいものだ。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   A = 4
   (r0, a, r1, r2, r3, x3) = (A*(sqrt(3) + 2)/4, A*(1 - sqrt(3))/4, sqrt(3)*A/4, A/8, sqrt(3)*A/16, sqrt(2)*3^(1/4)*A/8)
   @printf("r0 = %g;  a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g\n", r0, a, r1, r2, r3, x3)
   plot()
   circle(0, 0, r0, :gray)
   circle(r1, a + r1, r1)
   circle(-r1, a + r1, r1)
   circle(0, r0 - r2, r2, :blue)
   circle(x3, a + 2r1 + r3, r3, :green)
   circle(-x3, a + 2r1 + r3, r3, :green)
   circle(0, a + 2r1 - r3, r3, :green)
   b = sqrt(r0^2 - a^2)
   segment(-b, a, b, a, :magenta)
   c = sqrt(r0^2 - (a + 2r1)^2)
   segment(-c, a + 2r1, c, a + 2r1, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, a + r1, "甲円:r1,(r1,a+r1)", :red, :center, :top, delta=-delta)
       point(0, r0 - r2, "乙円:r2,(0,r0-r2) ", :black, :right, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その434）

算額（その434）

埼玉県大宮市 氷川神社 明治31年(1898)

山口正義(2015)：やまぶき2 第33号
https://yamabukiwasan.sakura.ne.jp/ymbk33.pdf

正方形内に，四分円 1 個，半円 2 個，円 1 子が入っている。円は四分円，半円と 3箇所で内接・外接している。円の直径が 1 寸のとき，正方形の一辺の長さを求めよ。

正方形の一辺の長さは 2r2 である。円の半径を r1，中心座標を (x1, y1) として，以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive,
     x1::positive, y1::positive, a::positive;
eq1 = (r2 - x1)^2 + y1^2 - (r1 + r2)^2
eq2 = (2r2 - x1)^2 + y1^2 - (2r2 - r1)^2
eq3 = x1^2 + (y1 - r2)^2 - (r2 - r1)^2
res = solve([eq1, eq2, eq3], (r2, x1, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (33*r1/8, 3*r1, 5*r1)

r2 = r1*33/8 ゆえ，正方形の一変の長さは 2r2 = 2r1*33/8 すなわち，円の直径が 1 寸のとき 33/8 = 4.125 寸である。

算額の答で（「埼玉の算額」より）は，四寸一分三厘五毛となっているようだが，誤植である。

   r1 = 0.5;  r2 = 2.0625;  x1 = 1.5;  y1 = 2.5;  正方形の一辺の長さ = 4.125

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1//2
   (r2, x1, y1) = (33*r1/8, 3*r1, 5*r1)
   a = 2r2
   @printf("r1 = %g;  r2 = %g;  x1 = %g;  y1 = %g;  正方形の一辺の長さ = %g\n", r1, r2, x1, y1, 2r2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(a, 0, 2r2, :blue, beginangle=90, endangle=180)
   circle(r2, 0, r2, :green, beginangle=0, endangle=180)
   circle(0, r2, r2, :magenta, beginangle=-90, endangle=90)
   circle(x1, y1, r1, :red)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, "(x1,y1)")
       point(0, r2, " r2")
       point(r2, 0, " r2", :green, :left, :bottom, delta=delta)
       point(2r2, 0, "2r2 ", :green, :right, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その433）

算額（その433）

埼玉県大宮市 氷川神社 明治31年(1898)

山口正義(2015)：やまぶき2 第33号
https://yamabukiwasan.sakura.ne.jp/ymbk33.pdf

外円の中に正方形 1 個，大円 2 個，小円 4 個が入っている。
小円の直径が 1 寸のとき，大円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
等円の半径と中心座標を r2, (r2, y2)
乙円の半径と中心座標を r3, (0, r3 - r0)
とし，以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive,
     x2::positive;
r0 = 2r1
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = r2/(r0 - x2) - (sqrt(Sym(2)) - 1)  # tand(45/2)
res = solve([eq1, eq2], (r1, x2))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(2)*r2 + 3*r2/2, sqrt(2)*r2 + 2*r2)

大円の半径は，小円の半径の (2*sqrt(2) + 3)/2 倍である。

sqrt(Sym(2))*r2 + 3*r2/2 |> simplify |> println

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

小円の直径が 1 寸のとき，大円の直径は 2.914213562373095 寸である。

r2 = 1
r2*(2*sqrt(2) + 3)/2

   2.914213562373095

   r1 = 1.45711;  x2 = 1.70711;  大円の直径 = 2.91421

using Plots

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

 

算額（その432）

算額（その432）

神奈川県横須賀市 浦賀叶大明神 文化11年(1814)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円に等円 2 個が内接し，その上下に甲円と乙円が大円に内接し，等円に概説している。大円の直径が 1 尺，甲円の直径が 5 寸のとき乙円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
等円の半径と中心座標を r2, (r2, y2)
乙円の半径と中心座標を r3, (0, r3 - r0)
とし，以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive,
     y2::negative, r3::positive;
eq1 = r2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = r2^2 + (y2 - r3 + r0)^2 - (r2 + r3)^2
eq3 = r2^2 + y2^2 - (r0 - r2)^2
res = solve([eq1, eq2, eq3], (r2, r3, y2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*r0*r1*(r0 - r1)/(r0 + r1)^2, r0*(r0 - r1)/(r0 + 3*r1), (r0^2 - 3*r0*r1)/(r0 + r1))

術および中村の現代解法は乙円の半径についてかなり複雑な式を提示している。
SymPy での結果は，乙円の半径は，r0*(r0 - r1)/(r0 + 3*r1) すなわち，「大円の半径と大円の半径と甲円の半径の差をかけたものを大円の半径と甲円の半径の3倍の和で割る」となり，藤井貞夫が「長野県和算研究会報No.5」で指摘したものになった。

大円の直径が 1 尺，甲円の直径が 5 寸のとき乙円の直径は 2 寸である。

2res[1][2](r0 => 10/2, r1 => 5/2).evalf() |> float |> println

   2.0

   r2 = 2.22222;  r3 = 1;  y2 = -1.66667;  乙円の直径 = 2

using Plots

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

 

算額（その431）

算額（その431）

東京都浅草 池田邸内 伽山

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円 2 個，中円 5 個，小円 4 個が図のように互いに内接・外接している。大円の直径が 21 寸のとき，小円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     x3::positive;
r2 = r1//2
eq1 = (2r2 - x3)^2 + x3^2 - (r2 + r3)^2
eq2 = x3^2 + (x3 + r2)^2 - (r1 + r3)^2
res = solve([eq1, eq2], (r3, x3))

   1-element Vector{Tuple{Sym, Sym}}:
    (3*r1/14, 4*r1/7)

大円の直径が 21 寸のとき，小円の直径は 21 * 3 / 14 = 4.5 寸である。

21 * 3 / 14

   4.5

r3 = 2.25;  x3 = 6;  小円の直径 = 4.5

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 21//2
   r2 = r1//2
   (r3, x3) = (3*r1/14, 4*r1/7)
   @printf("r3 = %g;  x3 = %g;  小円の直径 = %g\n", r3, x3, 2r3)
   plot()
   circle(0, r2, r1, :red)
   circle(0, -r2, r1, :red)
   circle(0, 0, r2, :green)
   circle42(0, 2r2, r2, :green)
   circle4(x3, x3, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r2, " r2\n 大円:r0", :red, :left, :top, delta=-delta)
       point(0, 2r2, " 2r2\n 中円:r2", :green, :left, :vcenter)
       point(x3, x3, " 小円:r3,(x3,x3)", :black, :left, :vcenter)
       point(2r2, 0, "2r2", :green, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

#算額

算額（その430）

算額（その430）

東京都浅草 池田邸内 伽山

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円の中に甲円 2 個と乙円 1 個，等円 3 個が入っている。等円は 1 本の元と長さの等しい 3 本の面（隣り合う 3 本の弦）に接している。面の長さが 1 寸のとき，乙円の直径を求めよ。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r1, y1) ただし y1 = a + r1
乙円の半径と中心座標を r2, (0, r0 - r2)
等円の半径と中心座標を r3, (0, a - r3), (2r3, a - r3)
弦と y 軸の交点座標を a
真ん中の面と y 軸の交点座標を b
面の長さの半分を c とし，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive,
     r3::positive, a::negative;
c = 1//2
b = a - 2r3
d = sqrt(r0^2 - a^2)
y1 = a + r1
eq1 = r1^2 + (r0 - r2 - y1)^2 - (r1 + r2)^2
eq2 = r1^2 + y1^2 - (r0 - r1)^2
eq3 = r0^2 - b^2 - c^2
eq4 = distance(c, b, d, a, 2r3, b+ r3) - r3^2
eq5 = (d - c)^2 + (a - b)^2 - (2c)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, r2, r3, a))

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, r1, r2, r3, a) = u
   return [
       r1^2 - (r1 + r2)^2 + (-a + r0 - r1 - r2)^2,  # eq1
       r1^2 + (a + r1)^2 - (r0 - r1)^2,  # eq2
       r0^2 - (a - 2*r3)^2 - 1/4,  # eq3
       -r3^2 + (2*r3 - 2*r3*(-4*a^2 + 4*r0^2 + 4*r3*sqrt(-a^2 + r0^2) + 2*r3 - 4*sqrt(-a^2 + r0^2) + 1)/(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))^2 + (a - r3 - (2*r3*(8*r3^2 + 8*r3*sqrt(-a^2 + r0^2) - 4*r3 - 2*sqrt(-a^2 + r0^2) + 1) + (a - 2*r3)*(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))/(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))^2,  # eq4
       4*r3^2 + (sqrt(-a^2 + r0^2) - 1/2)^2 - 1,  # eq5
   ]
end;

iniv = [big"65.0", 31, 27, 10, -32] ./ 40
res = nls(H, ini=iniv);
println(res);
   (BigFloat[2.176250899482821511100052865997767880198073191589329947230101745924818974139391, 0.9777711616568012201157129796477534133807314665018133584786185122552240330104753, 1.070019583175477769385089300350675253824657209000148549887247280185448490007896, 0.2236067977499789696409173668731276235440618359611525724270897245410529337440141, -1.670820393249936908922752100619382870632185507883457717281269173623140154153112], true)

r0 = 2.17625;  r1 = 0.977771; , r2 = 1.07002; , r3 = 0.223607; , a = -1.67082
乙円の直径 = 2.14004;  甲円の直径 = 1.95554;  等円の直径 = 0.447214
面 = 1;  弦 = 2d = 2.78885;  矢 = r0 + a = 0.505431

面の長さが 1 寸のとき，乙円の直径は 2.14004 寸である。

using Plots

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

術では，乙円径は 1.4221...寸としているが，この値は間違っている。
術の式通りに計算すると下記のように 2.14004...寸となる。

@syms 角, 元, 氏, 房, 面

面 = 1
角 = sqrt(80)
元 = sqrt(角 + 10)
氏 = 元 - sqrt(角 + 5.8 - sqrt(12.8))
房 = sqrt(2氏 * 元) - 氏
乙円径 = (元 * 房 * 面) / (2元 -(2sqrt((元 -2房)元) + 房))
@printf("面 = %g;  角 = %g;  元 = %g;  氏 = %g;  房 = %g\n", 面, 角, 元, 氏, 房)  
@printf("乙円径 = %g\n", 乙円径)  

   面 = 1;  角 = 8.94427;  元 = 4.3525;  氏 = 1.01086;  房 = 1.95554
   乙円径 = 2.14004

中村は，術と中村の導いた乙円径を表す式が同じであることを示しているが，(7)の分母の式において（またしても）符号を誤っているのでそれに基づいたのでは正しい答えにならない。"+ 房" は "- 房" である。

@printf("大円径 = %g\n", 元*面)

@printf("h(矢) = %g\n", 氏*面/2)

@printf("(7) の分子 = 甲円径 = %g\n", 房*面)
@printf("(7) の分母 = %g\n", (2元 - 2sqrt(元*(元 - 2房)) - 房)面)

x = (R*d/((2元 - 2sqrt(元*(元 - 2房)) - 房)面))
@printf("乙円径 = %g\n", x.evalf())

   大円径 = 4.3525
   h(矢) = 0.505431
   (7) の分子 = 甲円径 = 1.95554
   (7) の分母 = 3.97726
   乙円径 = 2.14004

以下は，中村による現代解法である。

@syms R,  # 大円径
     r,  # 等円径
     d,  # 甲円径
     x,  # 乙円径
     h,  # 矢
     a   # 面

a = 1

r = a/sqrt(Sym(5))
@printf("等円径 = %g\n", r.evalf())

R = sqrt(10 + 4sqrt(Sym(5)))a
@printf("大円径 = %g\n", R.evalf())

h = (sqrt(10 + 4sqrt(Sym(5))) - sqrt(4sqrt(Sym(5)) + 29//5 - 8sqrt(Sym(5))/5))a/2
@printf("h(矢) = %g\n", h.evalf())

eq3 = (d/2)^2 - h*(R - d) + h^2;
d = solve(eq3, d)[1]
@printf("甲円径 = %g\n", d.evalf())

x = (R*(2R - d) + 2R*sqrt(R*(R - 2d)))/(d + 4R)
@printf("乙円径 = %g\n", x.evalf())

   等円径 = 0.447214
   大円径 = 4.3525
   h(矢) = 0.505431
   甲円径 = 1.95554
   乙円径 = 2.14004

 

 

#算額

算額（その429）

算額（その429）

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

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

2 個の大円が交わっており，正方形 3 個，小円 2 個が入っている。大円の直径が 169.3 寸のとき，小円の直径を求めよ。

正方形の一辺の長さを 2a
大円の半径と中心座標を r1, (x1, 0), (-x1, 0)
小円の半径と中心座標を r2, (0, a + r2)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     a::positive;

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

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (2*r1*(4 - sqrt(6))/15, r1*(2*sqrt(6) + 7)/25, r1*(-4/25 + 6*sqrt(6)/25))

小円の半径は大円の半径の 2(4 - √6)/15 倍である。
小円の直径が 169.3 のとき，大円の直径は 35.00018487290774 である。

res[1][1](r1 => 169.3/2).evalf() |> println
res[1][2](r1 => 169.3/2).evalf() |> println
res[1][3](r1 => 169.3/2).evalf() |> println

   17.5000924364539
   40.2899445381277
   36.2198336143830

using Plots

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

#算額

算額（その428）

算額（その428）

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

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

外円の中に長方形が内接している。一辺の長さが 41.86 寸の正方形が 6個，大円が 2 個，小円が 4 個入っている。小円の直径を求めよ。

正方形の一辺の長さを a
大円の半径と中心座標を r1, (r0 - r1)
小円の半径と中心座標を r2, (r0 - 2r1 + r2, y)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r0::positive, r1::positive,
     r2::positive, y::positive;
@syms a, r0, r1, r2, y
eq1 = (sqrt(Sym(2))a/2 + a)^2 + (r0 - sqrt(Sym(2))a + a)^2 - r0^2
eq2 = (r0 - 2r1)^2 + (r0 - sqrt(Sym(2))a)^2 - r0^2
eq3 = (r0 - 2r1 + r2)^2 + y^2 - (r0 - r2)^2
eq4 = (r1 - r2)^2 + y^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (r0, r1, r2, y));

2 組の解が得られるが，2 番目のものが適解である。

names = ["r0", "r1", "r2", "y"]
for (i, name) in enumerate(names)
   println("\n", name)
   println(simplify(res[2][i]))
   println(res[1][i](a => 4186).evalf())
end

   r0
   a*(5 + 7*sqrt(2))/4
   15592.3214511641

   r1
   (-23922439916777731675282139353831061069516954103552970349841905*a^2 + 16915719487681281681493317896073504500369711065425523483418778*sqrt(2)*a^2 - 6986054008647705185868918429789392*sqrt(5)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4) + 4939886163250256057660601952272630*sqrt(10)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4))/(8*a*(4882633868649679478319412976906635182914530485393454527631129 - 3452543518573294933348514588454588356006400466465731641999825*sqrt(2)))
   1681.32809749726

   r2
   a*(205 + 151*sqrt(2))/1168
   1500.02961796760

   y
   -sqrt(-6848627247131292427665721780*sqrt(2)*a^2 + 9685421536530988262785792422*a^2 - 4*sqrt(10)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4))/(8*sqrt(8442762866313367764895228377 - 5969934874720155329318365400*sqrt(2)))
   3176.18761647798

小円の半径は a*(205 + 151*√2)/1168 なので，直径は a*(205 + 151*√2)/584 である。

a = 41.86 を代入すると，小円の直径は 30 寸あまりあり。

41.86*(205 + 151*√2)/584

   30.00059235935206

res2 = Float64.([res[2][i](a => 4186).evalf()/100 for i in 1:4])

   4-element Vector{Float64}:
    155.92321451164108
     16.81328097497264
     15.00029617967603
     31.761876164779796

r0 = 155.923;  r1 = 16.8133;  r2 = 15.0003;  y = 31.7619
小円の直径 = 30.0006

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r2, y) = res2
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  y = %g\n", r0, r1, r2, y)
   @printf("小円の直径 = %g\n", 2r2)
   a = 41.86
   (x0 , y0) = (r0 - 2r1, r0 - √2a)
   plot([x0, x0, -x0, -x0, x0], [-y0, y0, y0, -y0, -y0], color=:orange, lw=0.5)
   plot!([0, √2a/2, 0, -√2a/2, 0], [r0 - √2a, r0 - √2a + √2a/2, r0, r0 - √2a + √2a/2, r0 - √2a], color=:red, lw=0.5)
   rect(√2a/2, r0 - √2a, √2a/2 + a, r0 - √2a + a, :red)
   rect(-√2a/2, r0 - √2a, -√2a/2 - a, r0 - √2a + a, :red)
   plot!([0, √2a/2, 0, -√2a/2, 0], -[r0 - √2a, r0 - √2a + √2a/2, r0, r0 - √2a + √2a/2, r0 - √2a], color=:red, lw=0.5)
   rect(√2a/2, -r0 + √2a, √2a/2 + a, -r0 + √2a - a, :red)
   rect(-√2a/2, -r0 + √2a, -√2a/2 - a, -r0 + √2a - a, :red)
   circle(0, 0, r0, :black)
   circle(r0 - r1, 0, r1, :blue)
   circle(-r0 + r1, 0, r1, :blue)
   circle4(r0 - 2r1 + r2, y, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta)
       point(r0 - r1, 0, "r0-r1", :blue, :center, delta=-delta/2)
       point(r0 - 2r1 + r2, y, "(r0-2r1+r2,y) ", :green, :right, :vcenter)
       point(0, r0 - √2a, "r0-√2a ", :orange, :right, :top, delta=-delta/2)
       point(a/√2, r0-√2a, "(a/√2,r0-√2a)", :red, :center, delta=-3delta)
       point(a/√2 + a, r0-√2a, "(a/√2+a,r0-√2a)", :red, :center, delta=-delta/2)
       point(a/√2 + a, r0 - √2a + a, "(a/√2+a,r0-√2a+a)", :red, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その427）

算額（その427）

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

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

大円は鈎股弦（直角三角形）の，鈎と弦に接しており，中円，小円は弦と大円に接しかつそれぞれも互いに接している。鈎，股がそれぞれ　3000 寸，4000 寸，小円の直径が 456 寸のとき，大円の直径を求めよ。

鈎，股，弦の長さをそれぞれ鈎，股，弦
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, x2)
小円の半径と中心座標を r3, (x31, y31), (x32, y32)
として，連立方程式を解く。

solve() の能力的に一度では解けないので，まず最初に大円と中円の半径を求める。このためには，図を反時計回りに回転させ，弦（下図の紫の線）が水平になるように回転させて方程式を解くとよい。eq13 はr1 と鈎股弦の面積の関係式である（上図を参照）。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, r2::positive,
     r3::positive, x3::positive;

eq11 = x3^2 + (r1 - 2r2 + r3)^2 - (r1 - r3)^2
eq12 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq13 = r1*(鈎 + 股) + (r1 - 2r2)*弦 - 鈎*股
res1 = solve([eq11, eq12, eq13], (r1, r2, x3))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-2*r3*弦^2 - 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 - sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), -sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) + 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 - 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 - sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) + 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 + 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 + sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), -sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) - 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 + 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 + sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) - 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))

4 組の解が得られるが，そのうちの最初の 2 組が適解である。1 番目と 2 番目は，小円の位置（x3）が対称なもので実質的に同じである。

println("r1 = ", res1[2][1](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))
println("r2 = ", res1[2][2](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))
println("x3 = ", res1[2][3](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))

   r1 = 1250
   r2 = 300
   x3 = 120*sqrt(19)

算額の解としては大円の半径が 1250，直径は 2500 でここまででよい。

次いで，図を描くために，r1, r2 が既知として，小円の中心座標を求める。

using SymPy
@syms r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, x31::positive, y31::positive,
     x32::positive, y32::positive;

(r1, r2, r3) = (1250, 300, 456//2)
x2 = r1 + (r1 - r2)*3//5
y2 = r1 + (r1 - r2)*4//5
eq1 = (x2 - x31)^2 + (y2 - y31)^2 - (r2 + r3)^2
eq2 = (x2 - x32)^2 + (y2 - y32)^2 - (r2 + r3)^2
eq5 = (x31 - r1)^2 + (y31 - r1)^2 - (r1 - r3)^2
eq6 = (x32 - r1)^2 + (y32 - r1)^2 - (r1 - r3)^2
res2 = solve([eq1, eq2, eq5, eq6], (x31, y31, x32, y32))

   4-element Vector{NTuple{4, Sym}}:
    (8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5)
    (8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19))
    (96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5)
    (96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19))

4 組の解が得られるが，2 番目と 3 番目のものが適解である（小円の位置関係が逆になるだけ）。

(8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19)) |> println
(96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5) |> println

   (1358.3457014200953, 2266.2407239349286, 2195.2542985799046, 1638.5592760650716)
   (2195.2542985799046, 1638.5592760650716, 1358.3457014200953, 2266.2407239349286)

2 つの解をまとめると，以下のようになる。

   鈎 = 3000;  股 = 4000;  r3 = 228
   r1 = 1250;  r2 = 300;  x31 = 2195.25;  y31 = 1638.56;  x32 = 1358.35;  y32 = 2266.24
   大円の直径 = 2500

using Plots

function draw0(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x3) = (1250, 300, 456//2, 120*sqrt(19))
   plot()
   circle(0, 0, r1)
   circle(0, r1 - r2, r2, :blue)
   circle(x3, r1 - 2r2 + r3, r3, :green)
   circle(-x3, r1 - 2r2 + r3, r3, :green)
   x = sqrt(r1^2 - (r1 - 2r2)^2)
   segment(-x, r1 - 2r2, x, r1 - 2r2, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1-r2, " r1-r2", :blue)
       point(0, r1-2r2, " r1-2r2", :magenta, delta=-delta)
       point(x3, r1 - 2r2 + r3, "(x3,r1-2r2+r3)", :green, :center, :top, delta=-delta)
       point(r1, 0, "r1 ", :red, :right, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, r3) = (3000, 4000, 456//2)
   (r1, r2) = (1250, 300)
   (x31, y31, x32, y32) = res2[3]
   x2 = r1 + (r1 - r2)*3/5
   y2 = r1 + (r1 - r2)*4/5
   @printf("鈎 = %g;  股 = %g;  r3 = %g\n", 鈎, 股, r3)
   @printf("r1 = %g;  r2 = %g;  x31 = %g;  y31 = %g;  x32 = %g;  y32 = %g\n", r1, r2, x31, y31, x32, y32)
   @printf("大円の直径 = %.8g\n", 2r1)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x31, y31, r3, :green)
   circle(x32, y32, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, r1, " 大円:r1,(r1,r1)", :red, :left, :vcenter)
       point(r1 + (r1 - r2)*3/5, r1 + (r1 - r2)*4/5, " 中円:r2,(x2,y2)", :blue, :left, :vcenter)
       point(x31, y31, " 小円:r3,(x31,y31)", :green, :left, :vcenter)
       point(x32, y32, " 小円:r3,(x32,y32)", :green, :left, :bottom, delta=delta)
       point(股, 0, "股", :black, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;