算額（その570）

算額（その570）

群馬の算額 19-1 榛名神社 文化8年(1811)
http://takasakiwasan.web.fc2.com/gunnsann/g019-1.html

菱形の中に，短いほうの対角線の頂点を中心とする 2 つの弧が接しており，容円が 2 個入っている。菱形の対角線の長さが 4 寸，3 寸のとき，容円の直径を求めよ。

対角線の長さを 2a, 2b（a > b）とする。
円弧の半径と中心座標を r1, (0, b)
容円の半径と中心座標を r2, (x2, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

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

   1-element Vector{Tuple{Sym, Sym}}:
    (b*(a^3 + a*b^2 + b^2*sqrt(a^2 + b^2) - b*sqrt(a^4 + 2*a^3*sqrt(a^2 + b^2) + 2*a^2*b^2 + 2*a*b^2*sqrt(a^2 + b^2) + b^4))/(a^2*sqrt(a^2 + b^2)), b*(-a*b - b*sqrt(a^2 + b^2) + sqrt(a^4 + 2*a^3*sqrt(a^2 + b^2) + 2*a^2*b^2 + 2*a*b^2*sqrt(a^2 + b^2) + b^4))/a^2)

菱形の一辺の長さを「斜辺」と置くと，容円の半径の式は以下のようになる。

@syms 斜辺
res[1][1](sqrt(a^2 + b^2) => 斜辺) |> println

   b*(a^3 + a*b^2 + b^2*斜辺 - b*sqrt(a^4 + 2*a^3*斜辺 + 2*a^2*b^2 + 2*a*b^2*斜辺 + b^4))/(a^2*斜辺)

さらに sqrt() のなかみを因子分析すると以下のようになる。

a^4 + 2*a^3*斜辺 + 2*a^2*b^2 + 2*a*b^2*斜辺 + b^4 |> factor |> println

   (a^2 + b^2)*(a^2 + 2*a*斜辺 + b^2)

まとめると
b*(a^3 + a*b^2 + b^2*斜辺 - b*sqrt((a^2 + b^2)*(a^2 + 2*a*斜辺 + b^2)))/(a^2*斜辺)
のようになり，a, b に具体的な値を与えると答えが得られる。

(a, b) = (4, 3) .// 2
斜辺 = sqrt(Sym(a^2 + b^2))
b*(a^3 + a*b^2 + b^2*斜辺 - b*sqrt((a^2 + b^2)*(a^2 + 2*a*斜辺 + b^2)))/(a^2*斜辺) |> println

   87/32 - 9*sqrt(65)/32

容円の半径は (87 - 9√65)/32 寸，直径は (87 - 9√65)/16 ≒ 0.90248 寸である。

(87 - 9√65)/16

   0.9024800165820661

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

   容円の直径 = 0.90248;  a = 2;  b = 1.5;  r1 = 1.5;  r2 = 0.45124;  x2 = 1.24793

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (4, 3) .// 2
   r1 = b
   (r2, x2) = (87/32 - 9*sqrt(65)/32, -81/32 + 15*sqrt(65)/32)
   @printf("容円の直径 = %g;  a = %g;  b = %g;  r1 = %g;  r2 = %g;  x2 = %g\n", 2r2, a, b, r1, r2, x2)
   θ = atand(b/a)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:blue, lw=0.5)
   circle(0, b, r1, beginangle=180+θ, endangle=(360-θ))
   circle(0, -b, r1, beginangle=θ, endangle=(180-θ))
   circle(x2, 0, r2, :green)
   circle(-x2, 0, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
       point(0, b, " b = r1", :blue, :left, :bottom, delta=delta/2)
       point(x2, 0, "等円:r2\n(x2,0)", :green, :center, :bottom, delta=delta)
   end
end;

算額（その569）

算額（その569）

群馬の算額 19-8 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-8.html

長方形の中に直角三角形，等円，小円が入っている。小円の直径が 1 寸，直角三角形の股（直角を挟むに辺のうち長い方）の長さが 6 寸のとき，等円の直径を求めよ。

長方形の長辺，短辺の長さを a, b，直角三角形の頂点の座標を （0, c), (d, b) とする。
等円の半径と中心座標を r1, (r1, r1), (a - r1, b - r1)
小円の半径と中心座標を r2, (r2, b - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, d::positive, r1::positive, r2::positive
x = 6
r2 = 1//2
eq1 = d*r1- b*r2
eq2 = (b - c)*r1 - (a - d)*r2
eq3 = c + a - sqrt(a^2 + c^2) - 2r1
eq4 = (a - d)^2 + b^2 - x^2
eq5 = b + (a - d) - x - 2r1;
res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, c, d, r1))

   1-element Vector{NTuple{5, Sym}}:
    (28/5, 24/5, 33/10, 2, 6/5)

等円の半径は 6/5 寸，直径は 2.4 寸である。

その他のパラメータは以下の通りである。
a = 5.6;  b = 4.8;  c = 3.3;  d = 2;  r1 = 1.2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, c, d, r1) = res[1]  # (28/5, 24/5, 33/10, 2, 6/5)
   @printf("a = %g;  b = %g;  c = %g;  d = %g;  r1 = %g\n", a, b, c, d, r1)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:blue, lw=0.5)
   circle(r1, r1, r1)
   circle(a - r1, b - r1, r1)
   circle(r2, b - r2, r2, :magenta)
   plot!([0, a, d, 0], [c, 0, b, c], color=:green, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       # hline!([0], color=:black, lw=0.5)
       # vline!([0], color=:black, lw=0.5)
       point(d, b, "(d,b)", :blue, :center, :bottom, delta=delta/2)
       point(0, c, "  c", :blue, :left, :vcenter)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(a - r1, b - r1, "等円:r1,(a-r1,b-r1)", :red, :center, :top, delta=-delta/2)
       point(r1, r1, "等円:r1,(r1,r1)", :red, :center, :top, delta=-delta/2)
       point(r2, b - r2, "小円:r2\n(r2,b-r2)", :magenta, :center, :bottom, delta=delta/2)
   end
end;

算額（その568）

算額（その568）

群馬の算額 19-5 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-5.html

長方形内に斜線を引き，大円，小円を入れる。長方形の長辺の長さが 7 寸のとき，短辺の長さと小円の直径を求めよ。

長方形の短辺と長辺の長さを a, b とする。斜線と y 軸の交点座標を (0, c) とする。
大円の半径と中心座標を a/2, (a/2, b - a/2)
小円の半径と中心座標を r2, (r2, r2), (a - r2, y2)
とおき，以下の連立方程式を解く。
eq3, eq4 は何通りか考えられる。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r2::positive, y2::positive, c::positive
b = 7
eq1 = (a - r2 - a/2)^2 + (b - a/2 - y2)^2 - (a/2 + r2)^2
eq2 = a + c - sqrt(a^2 + c^2) - 2r2
eq3 = distance(a, 0, 0, c, r2, r2) - r2^2
eq3 = a*c - 2*a*r2 - 2*c*r2 + 2*r2^2
eq3 = (b - a/2 - c)/(a/2) - r2/y2
eq4 = distance(a, 0, 0, c, a - r2, y2) - r2^2;
eq4 = -a*r2^2 + a*y2^2 - 2*c*r2*y2;
eq4 = (a/2)/(b - a/2) - r2/y2;
# res = solve([eq1, eq2, eq3, eq4], (a, r2, y2, c))

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, r2, y2, c) = u
   return [
       (a/2 - r2)^2 - (a/2 + r2)^2 + (-a/2 - y2 + 7)^2,  # eq1
       a + c - 2*r2 - sqrt(a^2 + c^2),  # eq2
       -r2/y2 + 2*(-a/2 - c + 7)/a,  # eq3
       a/(2*(7 - a/2)) - r2/y2,  # eq4
   ]
end;

iniv = BigFloat[5.98, 0.76, 1.01, 1.78]
res = nls(H, ini=iniv)
   (BigFloat[5.999999999999999999999999999999999999999999999999999999999999999863119186219052, 0.7500000000000000000000000000000000000000000000000000000000000002629111561500711, 0.9999999999999999999999999999999999999999999999999999999999999991291513094263072, 1.75000000000000000000000000000000000000000000000000000000000000039296728020577], true)

短辺の長さは 6 寸，小円の直径は 1.5 寸である。
ちなみに，y2 = 1 寸, c = 1.75 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 7
   (a, r2, y2, c) = res[1]
   @printf("a = %g;  r2 = %g;  y2 = %g;  c = %g\n", a, r2, y2, c)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:blue, lw=0.5)
   circle(a/2, b - a/2, a/2)
   circle(r2, r2, r2, :magenta)
   circle(a - r2, y2, r2, :magenta)
   segment(0, c, a, 0)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a/2, b - a/2, "大円:a/2,(a/2,b-a/2)", :red, :center, delta=-delta/2)
       point(r2, r2, " 小円:r2\n (r2,r2)", :magenta, :center, delta=-delta/2)
       point(a - r2, y2, " 小円:r2\n (a-r2,y2)", :magenta, :center, delta=-delta/2)
       point(0, c, " c", :blue, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
   end
end;

 

算額（その567）

算額（その567）

群馬の算額 29-3 山名八幡宮 文化11年
http://takasakiwasan.web.fc2.com/gunnsann/g029-3.html

楕円の中に長径に沿って等円を描く。両端の円は楕円の頂点で楕円に接している。楕円の長径，短径（差渡し径）が 15 寸，5 寸のとき，等円の個数の最小値はいくつか。

長径，短径が a, b の楕円の曲率円（長径の端の 1 点で内接する半径が最大の円）の半径は，b^2/a である。

等円の個数の最小値は「楕円の長径/曲率円の半径」である。
a = 15/2, b = 5/12 のとき 9 個である。

a = 15//2
b = 5//2
println("曲率円の半径 = $(b^2/a)")
a/(b^2/a)

   曲率円の半径 = 5//6

   9//1

   a = 7.5;  b = 2.5;  曲率円の半径 = 0.833333

include("julia-source.txt");

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (15, 5) .// 2
   r = b^2/a
   @printf("a = %g;  b = %g;  曲率円の半径 = %g\n", a, b, b^2/a)
   plot()
   ellipse(0, 0, a, b, color=:red)
   for x in -a + r:2r:a - r
       circle(x, 0, r, :blue)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, b, " b", :red, :left, :bottom, delta=delta)
       point(a, 0, " a", :red, :left, :bottom, delta=delta)
       point(r - a, 0, "r-a", :blue, :center, :bottom, delta=delta)
       point(a - r, 0, "a-r", :blue, :center, :bottom, delta=delta)
   end
end;

算額（その566）

算額（その566）

群馬の算額 56-4改 八幡宮 天保5年
http://takasakiwasan.web.fc2.com/gunnsann/g056-4kai.html

正方形の中に 2 本の斜線を引き，大小 5 個の円を入れる。直角三角形の鈎（直角を挟む 2 辺のうち，短い方の辺）の長さが 1 寸のとき，大円の直径はいかほどか。

正方形の中心を原点とし，一辺の長さを 2a とする。
斜線と正方形の一辺が交差する座標を (b, -a)
大円の半径と中心座標を r1, (x1, a - r1)
小円の半径と中心座標を r2, (a - r2, r2 - a)
として，以下の連立方程式を解く

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, k:: positive,
     r1::positive, x1::positive, r2::positive

eq1 = (a - b) - k
eq2 = x1^2 + (a - r1)^2 - (r1 - r2)^2
eq3 = 2a + k - sqrt(4a^2 + k^2) - 2r2
eq4 = a*(a + b) - r1*sqrt(4a^2 + k^2)  # 中央の平行四辺形の面積を二通りの表現で
eq5 = distance(b, -a, a, a, x1, a - r1) - r1^2|>simplify|>numerator;
eq5 = a^3 - a^2*r1 - 2*a^2*x1 + a*b*r1 - a*r1^2 + a*r1*x1 + a*x1^2 - b*r1*x1;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, r1, x1, r2))

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, r1, x1, r2) = u
   return [
a - b - k,  # eq1
x1^2 + (a - r1)^2 - (r1 - r2)^2,  # eq2
2*a + k - 2*r2 - sqrt(4*a^2 + k^2),  # eq3
a*(a + b) - r1*sqrt(4*a^2 + k^2),  # eq4
a^3 - a^2*r1 - 2*a^2*x1 + a*b*r1 - a*r1^2 + a*r1*x1 + a*x1^2 - b*r1*x1,  # eq5
   ]
end;

k = 1
iniv = BigFloat[1.95, 0.945, 1.3, 0.236, 0.532]
res = nls(H, ini=iniv)

  (BigFloat[1.550697252562826473616653568524307307723291771022322587402262886914361142019489, 0.5506972525628264736166535685243073077232917710223225874022628869143611420202746, 0.9999999999999999999999999999999999999999999999667428781642976015715212111819139, 0.1775643993863705460811874612232267136362404546140429319571589691343659605198328, 0.4213839097383412741278163330894292515538227879203389202145579779541173446366729], true)

大円の直径は「鈎」の 2 倍である。

鈎 = k = 1 のときのパラメータの値。
a = 1.5507;  b = 0.550697;  r1 = 1;  x1 = 0.177564;  r2 =0.421384
鈎 = 1;  大円の直径 = 2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r1, x1, r2) = res[1]
   k = a - b
   @printf("a = %g;  b = %g;  r1 = %g;  x1 = %g;  r2 =%g\n", a, b, r1, x1, r2)
   @printf("鈎 = %g;  大円の直径 = %g\n", k, 2r1)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle(x1, a - r1, r1)
   circle(-x1, r1 - a, r1)
   circle(0, 0, r2, :green)
   circle(a - r2, r2 - a, r2, :green)
   circle(r2 - a, a - r2, r2, :green)
   segment(b, -a, a, a, :magenta)
   segment(-a, -a, -b, a, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, a - r1, "大円:r1,(x1,a-r1)", :red, :center, :bottom, delta=delta)
       point(a - r2, r2 - a, "小円:r2\n(a-r2,r2-a)", :green, :center, delta=-delta)
       point(a, 0, "a ", :blue, :right)
       point(0, -a, "-a ", :blue, :right, :bottom, delta=delta)
   end
end;

算額（その565）

算額（その565）

群馬の算額 51-3 稲荷神社 文政11年
http://takasakiwasan.web.fc2.com/gunnsann/g051-3.html

直角三角形内に全円と正方形が内接している。この直角三角形の弦を一辺とし，それと直角をなす辺の長さを前述の正方形の一辺の長さとする直角三角形を作る。この直角三角形に内接する円の直径が 5 寸，全円の直径が 7 寸のとき，もとの直角三角形の三辺（鈎，股，弦）の和を求めよ。

元の直角三角形の三辺を「鈎」，「股」，「弦」とする。
全円の半径と中心座標を r1, (r1, r1)
正方形の一辺の長さを a
上の直角三角形に内接する円の半径と中心座標を r2, (x2, y2)，上の頂点を (x, y)
として，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     鈎2::positive, 股2::positive, 弦2::positive,
     r1::positive, r2::positive, a::positive,
     x::positive, y::positive,
     x2::positive, y2::positive

鈎2 = a  # sqrt(x^2 +(y - 鈎)^2)
x = sqrt(a^2 - (y - 鈎)^2)
弦2 = sqrt((股 - x)^2 + y^2)
股2 = 弦

eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎 + 股 - 弦 - 2r1
eq3 = (鈎 - a)/a - 鈎/股
eq4 = 鈎2^2 + 股2^2 - 弦2^2
eq5 = 鈎2 + 股2 - 弦2 - 2r2
eq6 = distance(0, 鈎, 股, 0, x2, y2) - r2^2 |> numerator;
eq7 = distance(x, y, 股, 0, x2, y2) - r2^2 |> numerator;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (鈎, 股, 弦, a, y, x2, 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, y, x2, y2) = u
   return [
       -弦^2 + 股^2 + 鈎^2,  # eq1
       -2*r1 - 弦 + 股 + 鈎,  # eq2
       -鈎/股 + (-a + 鈎)/a,  # eq3
       a^2 - y^2 + 弦^2 - (股 - sqrt(a^2 - (y - 鈎)^2))^2,  # eq4
       a - 2*r2 + 弦 - sqrt(y^2 + (股 - sqrt(a^2 - (y - 鈎)^2))^2),  # eq5
       -r2^2 + (x2 - 股*(x2*股 - y2*鈎 + 鈎^2)/(股^2 + 鈎^2))^2 + (y2 - 鈎*(-x2*股 + y2*鈎 + 股^2)/(股^2 + 鈎^2))^2,  # eq6
       -r2^2 + (x2 - (x2*(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2) - y*(x2*y - y*股 + y2*股 - y2*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2)))/(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2))^2 + (-y*(-x2*股 + x2*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) + y*y2 + 股^2 - 股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2))/(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2) + y2)^2,  # eq7
   ]
end;

(r1, r2) = (7, 5) .// 2
iniv = BigFloat[9.41379, 16.89655, 19.7931, 6.27586, 15.2069, 3.45, 10.55]
res = nls(H, ini=iniv)

  (BigFloat[10.49999999999999999999999999999999999999999129919479365547439836695548116931581, 14.00000000000000000000000000000000000000001100960233399594221790231656113690163, 17.5000000000000000000000000000000000000000023087971276514166162692720423047526, 5.99999999999999999999999999999999999999999995377737351386752851972765331489346, 15.29999999999999999999999999999999999999999475386066794051341170073747364389192, 3.500000000000000000000000000000000000000000057743557772276350570950054497150227, 10.99999999999999999999999999999999999999999392286748897320175788883023235161591], true)

元の直角三角形の鈎，股，弦の和は 42 寸である。

その他のパラメータは以下の通り。
鈎 = 10.5;  股 = 14;  弦 = 17.5;  a = 6;  x = 3.6;  y = 15.3; , x2 = 3.5;  y2 =11

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (7, 5) .// 2
   (鈎, 股, 弦, a, y, x2, y2) = Float64.(res[1])
   x = sqrt(a^2 - (y -鈎)^2)
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  a = %g;  x = %g;  y = %g; , x2 = %g;  y2 =%g\n", 鈎, 股, 弦, a, x, y, x2, y2)
   @printf("元の直角三角形の鈎，股，弦の和 = %g\n", 鈎 + 股 + 弦)
   plot([0, 0, 股, x, 0, 股], [鈎, 0, 0, y, 鈎, 0], color=:blue, lw=0.5)
   plot!([0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, y2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       # hline!([0], color=:black, lw=0.5)
       # vline!([0], color=:black, lw=0.5)
       point(x, y, " (x,y)", :blue, :left, :vcenter)
       point(0, 鈎, " 鈎", :blue, :left, :vcenter, delta=delta/2)
       point(股, 0, "股", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
       point(0, a, " a", :magenta, :left, :bottom, delta=delta/2)
       point(r1, r1, "全円:r1,(r1,r1)", :red, :center, :top, delta=-delta/2)
       point(x2, y2, "円:r2,(x2,y2)", :green, :center, :top, delta=-delta/2)
   end
end;

算額（その564）

算額（その564）

一 群馬県高崎市石原町 清水寺 享保20年(1735)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

群馬の算額 1 清水寺
http://takasakiwasan.web.fc2.com/yattemimasennka1.html

キーワード：円6個，外円

外円内に大円 2 個，中円 1 個，小円 2 個が入っている。大円，中円，小円の直径が与えられたとき，外円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, y1)
中円の半径と中心座標を r2, (0, y2)
小円の半径と中心座標を r3, (r3, y3)
とおき，以下の連立方程式を解く。

流石に，未知数が 4 個とはいえ，他の 3 変数を変数のまま連立方程式を解くのは無理なようで，r1, r2, r3 を数値で与えれば解ける。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, y1::negative, r2::positive, y2::positive, r3::positive, y3::positive
(r1, r2, r3) = (3, 2, 1)
eq1 = r1^2 + y1^2 - (R - r1)^2
eq2 = r3^2 + y3^2 - (R - r3)^2
eq3 = r1^2 + (y2 - y1)^2 - (r1 + r2)^2
eq4 = r3^2 + (y3 - y2)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3, eq4], (R, y1, y2, y3))

   1-element Vector{NTuple{4, Sym}}:
    (22/7 + 16*sqrt(2)/7, -8/7 - 2*sqrt(2)/7, 20/7 - 2*sqrt(2)/7, 12*sqrt(2)/7 + 20/7)

大円，中円，小円の直径がそれぞれ 6, 4, 2 のとき，外円の直径は (22 + 16√2)/7 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, y1, y2, y3) = (22/7 + 16*sqrt(2)/7, -8/7 - 2*sqrt(2)/7, 20/7 - 2*sqrt(2)/7, 12*sqrt(2)/7 + 20/7)
   plot()
   circle(0, 0, R)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(0, y2, r2, :green)
   circle(r3, y3, r3, :magenta)
   circle(-r3, y3, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r1, y1, "大円:r1,(r1,y1)", :blue, :center, delta=-delta/2)
       point(0, y2, " 中円:r2\n (0,y2)", :green, :left, :vcenter)
       point(r3, y3, " 小円:r3,(r3,y3)", :magenta, :left, :vcenter)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
   end
end;

 

算額（その563）

算額（その563）

群馬の算額 20 旧福寿庵
http://takasakiwasan.web.fc2.com/gunnsann/g020.html

長方形内に 2 本の斜線，上円 1 個，下円　2 個が入っている。長方形の縦×横が 336 平方寸，上円と下円の直径の積が 63 平方寸のとき，上円の直径はいかほどか。

長方形の短辺，長辺の長さを a, b とおく。
上円の半径と中心座標を r1, (a/2, b - r1)
下円の半径と中心座標を r2, (r2, r2)
とおき，以下の連立方程式を解く。
簡単な連立方程式に見えるが SymPy では解けないので数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, a::positive, b::positive

弦 = sqrt(b^2 + a^2/4)
eq1 = r1*弦 - a*(b - r1)/2
eq2 = b + a/2 - 弦 - 2r2
eq3 = a*b - 336
eq4 = 2r1*2r2 - 63;
# res = solve([eq1, eq2, eq3, eq4], (r1, r2, a, b))

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, a, b) = u
   return [
       -a*(b - r1)/2 + r1*sqrt(a^2/4 + b^2),  # eq1
       a/2 + b - 2*r2 - sqrt(a^2/4 + b^2),  # eq2
       a*b - 336,  # eq3
       4*r1*r2 - 63,  # eq4
   ]
end;

iniv = BigFloat[5, 3.1, 15, 25]
res = nls(H, ini=iniv)

   (BigFloat[5.250000000000000000000000000000000000000000009939436168066665434919743445334633, 2.999999999999999999999999999999999999999999989678292454450465893181229233646876, 13.99999999999999999999999999999999999999999983811297245359217105296244127297891, 24.00000000000000000000000000000000000000000008495905738224684384339759581898053], true)

上円の半径は 5.25 寸，直径は 10.5 寸である。

その他のパラメータは以下の通り。
上円の直径 = 10.5;  r1 = 5.25;  r2 = 3;  a = 14;  b = 24

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

 

 

算額（その562）

算額（その562）

群馬の算額 104-8 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-8.html

直角三角形と全円の隙間に 3 個の等円を入れる。全円の直径が与えられたとき，等円の直径を求めよ。

全円の半径と中心座標を r1, (0, 0)
等円の半径と中心座標を r2, (1 - 3r2, r2 - r1), (r1 - r2, 0), (x, y)
とおき，以下の方程式を解く。

図を描くのでない限り，斜辺に接する等円の中心座標を求める必要はない。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive

eq = (r1 - 3r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve(eq, r2)
res[1] |> println

   r1/9

等円の半径は全円の半径の 1/9 である。
全円の半径が 1/2 のとき，小円の半径は 1/18，直径は 1/9 である。

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

算額（その561）

算額（その561）

算額（その452）と重複

群馬の算額 104-4 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-4.html

正方形の中に大円，中円，小円を入れる。正方形の一辺の長さが与えられたとき，小円の直径を求めよ。

正方形の一辺の長さを a とする。
大円の半径と中心座標を r1, (r1, 0)
中円の半径と中心座標を r2, (0, 0), (a/2 - r2, 0)
小円の半径と中心座標を r3, (0, a/2 - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, r2::positive, r3::positive

r2 = a/6
r1 = 2r2
eq = r2^2 + (a/2 - r3)^2 - (r1 + r3)^2
res = solve(eq, r3)
res[1] |> println

   a/10

小円の半径は正方形の一辺の長さの 1/10 である。
正方形の一辺の長さが 1 のとき，小円の半径は 0.1，直径は 0.2 である。

その他のパラメータは以下の通りである。
小円の直径 = 0.2;  a = 1;  r1 = 0.333333;  r2 = 0.166667;  r3 = 0.1

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

算額（その560）

算額（その560）

群馬の算額 104-7改 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-7kai.html

正三角形に内接する大円，大円に外接し頂点を通る中円，中円と社編が作る弓形の中の小円がある。大円の直径が与えられたとき，小円の直径を求めよ。

大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (0, 2r1 + r2)
小円の半径と中心座標を r3, (x, y)
とおき，以下の連立方程式を解く。
なお，eq1 と (eq2, eq3) は独立なので，まず eq1 から r3 を求め，(x, y) は (eq2, eq3) の連立方程式を解かなくても三角比から求めることもできる。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive

r2 = r1/2
eq1 = (r2 - 2r3)/r2 - r1/2r1
eq2 = x^2 + (y - 2r1 - r2)^2 - (r2 - r3)^2
eq3 = (y - 2r1 - r2)/x - 1/sqrt(Sym(3))
res = solve([eq1, eq2, eq3], (r3, x, y))
res[1] |> println

   (r1/8, 3*sqrt(3)*r1/16, 43*r1/16)

小円の半径は大円の半径の 1/8 である。
大円の直径が 1 のとき，小円の直径は 0.125 である。

その他のパラメータは以下の通りである。
小円の直径 = 0.125;  r1 = 0.5;  r2 = 0.25;  r3 = 0.0625;  x = 0.16238;  y = 1.34375

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   r2 = r1/2
   (r3, x, y) = r1 .* (1/8, 3*sqrt(3)/16, 43/16)
   @printf("小円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x = %g;  y = %g\n", 2r3, r1, r2, r3, x, y)
   plot([√3r1, 0, -√3r1, √3r1], [0, 3r1, 0, 0], color=:blue, lw=0.5)
   circle(0, r1, r1)
   circle(0, 2r1 + r2, r2, :green)
   circle(x, y, r3, :magenta)
   circle(-x, y, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1, " 大円:r1,(0,r1)", :red, :left, :vcenter)
       point(0, 2r1 + r2, " 中円:r2,(0,2r1+r2)", :green, :left, :vcenter)
       point(x, y, "   小円:r3,(x,y)", :magenta, :left, :vcenter)
       point(√3r1, 0, "√3r1  ", :blue, :right, :bottom, delta=delta/2)
       point(0, 2r1, " 2r1", :green, :left, :top, delta=-delta)
       point(0, 3r1, " 3r1", :green, :left, :bottom)
   end
end;

 

算額（その559）

算額（その559）

群馬の算額 104-3 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-3.html

2 個の外円が交わり，大円，小円が入っている。斜線は大円と小円の共通接線である。外円の直径が与えられたとき，小円の直径を求めよ。

外円の半径と中心座標を R, (r1, 0)
大円の半径と中心座標を r1, (0, 0); r1 = R/2
小円の半径と中心座標を r2, (r1 + r2, 0)
とおき，以下の連立方程式を解く。
eq1 と eq2 は独立で，それぞれ r2, x を求めることができる。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, x::positive, y::positive

r1 = R/2
y = sqrt(R^2 - (x - r1)^2)
eq1 = r1/3r1 - r2/(3r1 - r1 - r2)
eq2 = y/sqrt((R + r1 + x)^2 + y^2) - r1/(R + r1)
solve([eq1, eq2], (r2, x))

   1-element Vector{Tuple{Sym, Sym}}:
    (R/4, R*(5 + 4*sqrt(10))/18)

小円の半径は外円の半径の 1/4 である。
外円の直径が 1 のとき，小円の直径は 0.25 である。

その他のパラメータは以下の通りである。
小円の直径 = 0.25;  R = 0.5;  r1 = 0.25;  r2 = 0.125;  x = 0.490253;  y = 0.438496

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1/2
   (r1, r2, x) = R .* (1/2, 1/4, (5 + 4*sqrt(10))/18)
   y = sqrt(R^2 - (x - r1)^2)
   @printf("小円の直径 = %g;  R = %g;  r1 = %g;  r2 = %g;  x = %g;  y = %g\n",
       2r2, R, r1, r2, x, y)
   plot()
   circle(r1, 0, R)
   circle(-r1, 0, R)
   circle(0, 0, r1, :green)
   circle(r1 + r2, 0, r2, :blue)
   circle(-r1 - r2, 0, r2, :blue)
   plot!([x, -R - r1, x], [y, 0, -y], color=:magenta, lw=0.5)
   plot!([-x, R + r1, -x], [y, 0, -y], color=:magenta, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r1, 0, "r1 ", :green, :right, :bottom, delta=delta/2)
       point(r1 + r2, 0, "r1+r2", :blue, :center, :bottom, delta=delta/2)
       point(r1 + R, 0, "r1+R ", :black, :right, :bottom, delta=delta/2)
       point(x, y, " (x,y)", :magenta, :left, :bottom)
   end
end;

 

 

算額（その558）

算額（その558）

群馬の算額 104-6 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-6.html

外円，大円，中円，小円，等円を図のように配置する。外円の直径が与えられたとき，小円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, y), (r1, y) ; r1 = 2r2
中円の半径と中心座標を r2, (0, R - r2), (r2, y)
小円の半径と中心座標を r3, (0, R - 2r3 - r3), (0, 3r3 - R), (0, r3 - R)
等円の半径と中心座標を r4, (x4, y4)
として，以下の連立方程式を解く。

等円のパラメータは他のパラメータと独立なので，最初に 小円と中円の半径を求める

include("julia-source.txt");

using SymPy

@syms R, r1, y::negative, r2::positive, r3, r4, x4, y4

y = r3 - r2
r1 = 2r2
eq4 = r1^2 + (R - 2r2 - r3 - y)^2 - (r1 + r3)^2
eq5 = r1^2 + (r1 - r3)^2 - (r1 + r3)^2
solve([eq4, eq5], (r2, r3))

   2-element Vector{Tuple{Sym, Sym}}:
    (2*R/7, R/7)
    (2*R, R)

最初のものが適解である。
中円と小円の半径は，外円の半径の 2/7, 1/7 である。
算額の問に対する答えはここまでで完了する。

図を描くために，次に r2, r3 が既知として 等円の半径と中心座標を求める。

@syms R, r1, y::negative, r2::positive, r3, r4, x4, y4

(r2, r3) = (2R/7, R/7)
y = r3 - r2
r1 = 2r2
eq1 = (x4 - r1)^2 + (y4 - y)^2 - (r1 + r4)^2
eq2 = x4^2 + y4^2 - (R - r4)^2
eq3 = x4^2 + (R - r2 - y4)^2 - (r2 + r4)^2
solve([eq1, eq2, eq3], (r4, x4, y4))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (5*R*(39 - sqrt(273))/546, 8*sqrt(273)*R/273, R*(3*sqrt(273)/182 + 5/14))
    (5*R*(sqrt(273) + 39)/546, -8*sqrt(273)*R/273, R*(5/14 - 3*sqrt(273)/182))

最初のものが適解である。

他のパラメータの数値解は以下の通りである。
小円の直径 = 0.285714;  R = 1;  r1 = 0.571429;  r2 = 0.285714;  r3 = 0.142857;  r4 = 0.50845;  x4 = -0.484182;  y4 = 0.0847905

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   (r2, r3) = R .* (2/7, 1/7)
   (r4, x4, y4) = R .*(5(sqrt(273) + 39)/546, -8*sqrt(273)/273, (5/14 - 3*sqrt(273)/182))
   y = r3 - r2
   r1 = 2r2
   @printf("小円の直径 = %g;  R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  x4 = %g;  y4 = %g\n",
       2r3, R, r1, r2, r3, r4, x4, y4)
   plot()
   circle(0, 0, R, :orange)
   circle(0, y, r1)
   circle(r1, y, r1)
   circle(-r1, y, r1)
   circle(0, R - r2, r2, :blue)
   circle(r2, y, r2, :blue)
   circle(-r2, y, r2, :blue)
   circle(0, R - 2r2 - r3, r3, :green)
   circle(0, 3r3 - R, r3, :green)
   circle(0, r3 - R, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(R, 0, " R", :black, :left, :bottom, delta=delta/2)
       point(0, 3r3 - R, " 小円:r3,(0,3r3-R)", :green, :left, :vcenter)
       point(0, r3 - R, " 小円:r3,(0,r3-R)", :green, :left, :vcenter)
       point(0, R - 2r2 - r3, " 小円:r3,(0,R-2r2-r3)", :green, :left, :vcenter)
       point(0, R - r2, " 中円:r2,(0,R-r2)", :blue, :left, :vcenter)
       point(r2, y, "中円:r2\n(r2,y)", :blue, :center, :top, delta=-delta)
       point(r1, y, " 大円:r1\n (r1,y)", :red, :left, :vcenter)
       point(0, y, "y ", :blue, :right, :vcenter)
   end
end;

算額（その557）

算額（その557）

群馬の算額 104-2 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-2.html

3 個の正三角形と大円，小円が図のように配置されている。正三角形の一辺の長さが与えられたとき，小円の直径を求めよ。

正三角形の一辺の長さを a とする。
大円の半径と中心座標を r1, (0, √3a/2); r1 = a/2
小円の半径と中心座標を r2, (0, √3a/2 - r1 + r2)
とおき，方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1, r2, a

r1 = a/2
eq = r2/(r1 - r2) - 1//2
solve(eq, r2)[1] |> println

   a/6

小円の半径は正三角形一辺の長さの 1/6 である。
小円の直径は正三角形一辺の長さの 1/3 である。

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

算額（その556）

算額（その556）

群馬の算額 104-1 長谷寺 文久元年
http://takasakiwasan.web.fc2.com/gunnsann/g104-1.html

外円内に弦を引き，区画された領域に大円 2 個，小円 3 個を入れる。大円の直径が与えられたとき，小円の直径を求めよ。

大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (0, r2), (x2, r2)
斜めの弦と円の交点座標を (x, y)
とおき，方程式を解く。

include("julia-source.txt");

using SymPy

@syms R, r1, r2, x2, x

R = 2r1
eq1 = (R - r2)/(2R - r1) - r2/r1
eq2 = x2^2 + r2^2 - (R - r2)^2
eq3 = x/sqrt(x^2 + (2r1 + sqrt(R^2 - x^2))^2) - r2/(R - r2)
solve([eq1, eq2, eq3], (r2, x2, x))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r1/2, -sqrt(2)*r1, 8*sqrt(2)*r1/9)
    (r1/2, sqrt(2)*r1, 8*sqrt(2)*r1/9)

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

小円の半径は大円の半径の 1/2 である。
また，外円の半径は大円の半径の 2 倍である。

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

   小円の直径 = 1;  r1 = 1;  r2 = 0.5;  x2 = 1.41421;  x = 1.25708;  y = 1.55556

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