算額（その900）

算額（その900）

七四 加須市大字外野 棘脱地蔵堂 明治7年(1874)埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に正三角形と大円，小円をそれぞれ 2 個ずつ入れる。外円の直径が 3 寸のとき小円の直径はいかほどか。

算額（その899）にもう一種類の円を加えたものであるが，SymPy の性能では数式解を求めることができない。

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

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive, y2::positive
eq1 = (R - r1)cosd(Sym(30)) - r1
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = (x2 - R + r1)^2 + y2^2 - (r1 + r2)^2
eq4 = dist2(0, 0, R*cosd(Sym(60)), R*sind(Sym(60)), x2, y2, 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 Float64.(v), r.f_converged
end;

function H(u)
   (r1, r2, x2, y2) = u
   return [
       -r1 + sqrt(3)*(R - r1)/2,  # eq1
       x2^2 + y2^2 - (R - r2)^2,  # eq2
       y2^2 - (r1 + r2)^2 + (-R + r1 + x2)^2,  # eq3
       -r2^2 + 3*x2^2/4 - sqrt(3)*x2*y2/2 + y2^2/4,  # eq4
   ]
end;

R = 3/2
iniv = BigFloat[0.7, 0.24, 0.83, 0.95]
res = nls(H, ini=iniv)

   ([0.6961524227066319, 0.2460188097551155, 0.8278641121314027, 0.9418650844642568], true)

外円の直径が 3 寸のとき，小円の直径は 0.492037619510231 寸である。

その他のパラメータは以下のとおりである。
R = 1.5; r1 = 0.696152;  r2 = 0.246019;  x2 = 0.827864;  y2 = 0.941865;  x = 0.75;  y = 1.29904

算額の「答」，「術」ともに不適切であろう。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 3/2
   (r1, r2, x2, y2) = res[1]
   (x, y) = R .* (cosd(60), sind(60))
   @printf("外円の直径が %g のとき，小円の直径は %.15g である。\n", 2R, 2r2)
   @printf("R = %g; r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x = %g;  y = %g\n", R, r1, r2, x2, y2, x, y)
   plot([x, -x, x, -x, x], [y, y, -y, -y, y], color=:blue, lw=0.5)
   circle(0, 0, R, :green)
   circle2(R - r1, 0, r1, :magenta)
   circle4(x2, y2, r2)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R - r1, 0, "大円:r1,(0,R-r1)", :magenta, :center, delta=-delta)
       point(x2, y2, "小円:r2,(x2,y2)", :black, :center, delta=-delta)
       point(x, y, "(x,y)", :blue, :left, :bottom, delta=delta)
   end
end;

算額（その899）

算額（その899）

七四 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に正三角形と円を 2 個ずつ入れる。外円の直径が 3 寸のとき等円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (R - r, 0)
とおき，以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive
eq1 = (R - r)cosd(Sym(30)) - r
res = solve(eq1, r)[1] |> simplify
res |> println
2res(R => 3/2).evalf() |> N

   R*(-3 + 2*sqrt(3))

   1.3923048454132638

等円の半径は，外円の半径の (2√3 - 3) 倍である。
外円の直径が 3 寸のとき，等円の直径は 3(2√3 - 3) = 1.3923048454132632 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 3/2
   r = R*(2√3 - 3)
   @printf("外円の直径が %g のとき，等円の直径は %.15g である。\n", 2R, 2r)
   (x, y) = R .* (cosd(60), sind(60))
   plot([x, -x, x, -x, x], [y, y, -y, -y, y], color=:blue, lw=0.5)
   circle(0, 0, R, :green)
   circle2(R - r, 0, r)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R - r, 0, "等円:r,(0,R-r)", :black, :left, delta=-delta)
       point(x, y, "(x,y)", :blue, :left, :bottom, delta=delta)
   end
end;

算額（その898）

算額（その898）

七四 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

大円 1 個と，甲円 4 個，乙円 8 個がある。甲円の直径が 1 寸のとき，乙円の直径はいかほどか。

大円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, ((R + r2)/√2, (R + r2)/√2)
とおき以下の連立方程式を解く。

1.　数式解を求める

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive
eq1 = r1^2 + (R - r2 - r1)^2 - (r1 + r2)^2
eq2 = 2((R + r2)/√Sym(2) - r1)^2 - (r1 - r2)^2
res = solve([eq1, eq2], (R, r2))[1];  # 1 of 3

3 組の解が得られるが，最初のものが適解である。
解は虚数解として得られるが，虚部は限りなく 0 に近いので，実部だけを取ればよい（あるいは，絶対値を取る）。

外円の半径 R = 0.975400964516989, 乙円の半径 r2 = 0.115852908334779

res[1](r1 => 1/2).evalf(), res[2](r1 => 1/2).evalf()

   (0.975400964516989 - 2.10363882672262e-33*I, 0.115852908334779 + 0.e-23*I)

乙円の直径は 0.115852908334779*2 = 0.231705816669558 である。

術は「置六個八分七厘壱毛減四個開平方減貳個六分貳厘以除四個」とあるが，(2.62 - sqrt(6.871 - 4))/4 と不思議な数と不思議な演算で 0.23139936260671184 を得て，答で「乙円径二分三厘有奇」としているが，いずれも不適切である。連立方程式で得られる式はとてつもなく長い（しかも，微細ではあるが虚部を含む）。

(2.62 - sqrt(6.871 - 4))/4

   0.23139936260671184

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   (R, r2) = (0.975400964516989, 0.115852908334779)
   @printf("甲円の直径が %g のとき，乙円の直径は %.15g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  R = %g;  r2 = %g\n", r1, R, r2)
   plot()
   circle(0, 0, R, :blue)
   circle4((R + r2)/√2, (R + r2)/√2, r2)
   circle42(0, R - r2, r2)
   circle4(r1, r1, r1, :green)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;

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 Float64.(v), r.f_converged
end;

function H(u)
   (R, r2) = u
   return [
       r1^2 - (r1 + r2)^2 + (R - r1 - r2)^2,  # eq1
       2*(-r1 + sqrt(2)*(R + r2)/2)^2 - (r1 - r2)^2,  # eq2
   ]
end;

r1 = 1/2
iniv = BigFloat[1, 0.5]
res = nls(H, ini=iniv)

   ([0.9754009645169893, 0.11585290833477907], true)

外円の半径 R = 0.975400964516989, 乙円の半径 r2 = 0.115852908334779

 

算額（その897）

算額（その897）

七四 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

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

外円の半径と中心座標を R, (0, 0)
1/3 円の半径は外円の半径と同じで，中心座標は外円の円周上にあり，(x0, y0), (-x0, y0), (0, -R); x0 = R*cosd(Sym(30)), y0 = R*sind(Sym(30))
大円の半径と中心座標を r1, (x1, y1)
小円の半径と中心座標を r2, (0, y2); y2 = y0
と置く。

小円の半径は図をよく見れば，方程式を解くまでもなく以下のように求める。
小円の半径は，外円の半径の (1 - √3/2) 倍である。
外円の直径が 2 寸のとき，小円の直径は 0.2679491924311228 寸である。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, y2::positive

R - R*cosd(Sym(30)) |> simplify |> println
2*(1 - √3/2) |> println

   R*(2 - sqrt(3))/2
   0.2679491924311228

図を描くために大円の半径を求める。中心座標は簡単に求めることができる。

x0 = R*cosd(Sym(30))
y0 = R*sind(Sym(30))
y1 = (R - r1)*sind(Sym(30))
x1 = (R - r1)*cosd(Sym(30))
eq1 = (x0 + x1)^2 + (y0 - y1)^2 - (R + r1)^2
res = solve(eq1, r1)[1]
res |> simplify |> println

   2*R/5

大円の半径は外円の半径の 2/5 倍である。

外円の直径が 2 寸のとき，その他のパラメータは以下のとおりである。
R = 1;  r1 = 0.4;  x1 = 0.519615;  y1 = 0.3;  r2 = 0.133975; y2 = 0.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 2//2
   r1 = 2R/5
   x0 = R*cosd(30)
   y0 = R*sind(30)
   y1 = (R - r1)*sind(30)
   x1 = (R - r1)*cosd(30)
   y2 = R/2
   r2 = R - x0
   @printf("外円の直径が %g のとき，小円の直径は %.15g である。\n", 2R, 2r2)
   @printf("R = %g;  r1 = %g;  x1 = %g;  y1 = %g;  r2 = %g;  y2 = %g\n", R, r1, x1, y1, r2, y2)
   plot()
   circle(0, 0, R, :blue)
   circle(x0, y0, R, beginangle=150, endangle=270)
   circle(-x0, y0, R, beginangle=270, endangle=390)
   circle(0, -R, R, beginangle=30, endangle=150)
   rotate(x1, y1, r1, :green)
   rotate(0, y2, r2, :magenta)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(-x0, y0, "(-x0,y0)")
       point(x1, y1, "大円:r1,(x1,y1)", :green, :center, delta=-delta/2)
       point(0, y2, " 小円:r2,(0,y2)", :magenta, :left, :bottom, delta=delta/2)
   end
end;

算額（その896）

算額（その896）

七四 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に，正三角形と，大きさの異なる 4 種類の円が入っている。それぞれの円の直径を求めよ。

注：問，答，術共に欠損文字があり，更に円の名前も不明瞭で，大きさの順にもなっていないので，大きい順に甲，乙，丙，丁とし，それぞれの直径を外円の直径に対する比で求める。当然ながら，ある円の直径を特定の値としたとき別の円の直径を求めるのは容易なことである。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 0); r1 = R/2
乙円の半径と中心座標を r2, (x2, r1 - r2)
丙円の半径と中心座標を r3, (x3, r1 + r3)
丁円の半径と中心座標を r4, (0, r1 + r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive, r4::positive
r1 = R/2
eq1 = (R - r1 - r4) - 2r4
eq2 = x2^2 + (r1 - r2)^2 - (R - r2)^2
eq3 = x3^2 + (r1 + r3)^2 - (R - r3)^2
eq4 = dist2(R√Sym(3)/2, -R/2, 0, R, x2, r1 - r2, r2)
eq5 = dist2(R√Sym(3)/2, -R/2, 0, R, x3, r1 + r3, r3)
res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, r3, x3, r4))[1];  # 1 of 4

4 組の解が得られるが，最初のものが適解である。それぞれは以下のように簡約化できる。

for i in 1:5
   res[i] |> simplify |> println
end

   R*(-1 + sqrt(3))/3
   R*(6 - sqrt(3))/6
   R*(-5 + 3*sqrt(3))
   3*R*(2 - sqrt(3))/2
   R/6

円についてだけいえば，大きい順に外円の直径を a とすれば，それぞれの直径は
甲円 a/2
乙円 a(√3 - 1)/3
丙円 a(3√3 - 5)
丁円 a/6
である。

外円の直径が 12 寸のとき，丁円の直径は 2 寸である。

その他のパラメータは以下のとおりである。

R = 6;  r1 = 3;  r2 = 1.4641;  x2 = 4.26795;  r3 = 1.17691;  x3 = 2.41154;  r4 = 1

この結果から見ると，算額の欠落文字を含む「只云□圓径□寸問丙圓径」に対して「丙円□□□分七厘有奇」は，「只云丁圓径一寸問丙圓径」に対して「丙円径一寸一分七厘有奇」なのであろう。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 3
   r1 = R/2
   (r2, x2, r3, x3, r4) = (
       R*(-1 + sqrt(3))/3,
       R*(6 - sqrt(3))/6,
       R*(-5 + 3*sqrt(3)),
       3*R*(2 - sqrt(3))/2,
       R/6)
   @printf("外円の直径が %g 寸のとき，丁円の直径は %g 寸\n", 2R, 2r4)
   @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g\n", R, r1, r2, x2, r3, x3, r4)
   plot([√3R/2, 0, -√3R/2, √3R/2], [-R/2, R, -R/2, -R/2], color=:green, lw=0.5)
   circle(0, 0, R, :blue)
   circle(0, 0, r1, :red)
   circle2(x2, r1 - r2, r2, :magenta)
   circle2(x3, r1 + r3, r3, :brown)
   circle(0, r1 + r4, r4, :orange)
   segment(-√3R/2, r1, √3R/2, r1, :black)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, 0, "甲円:r1,(0,0)", :red, :center, delta=-delta/2)
       point(x2, r1 - r2, "乙円:r2\n(x2,r1-r2)", :magenta, :center, delta=-delta/2)
       point(x3, r1 + r3, "丙円:r3\n(x3,r1+r3)", :brown, :center, delta=-delta/2)
       point(0, r1 + r4, "丁円:r4\n(0,r1+r4)", :black, :center, delta=-delta/2)
       point(0, r1, " r1", :red, :center, delta=-delta/2)
       point(0, -r1, " -r1", :red, :center, delta=-delta/2)
   end
end;

算額（その895）

算額（その895）

七二 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

仁円 1 個，義円 2 個，禮円 4 個と 2 本の斜線がある。義円の直径が 3 寸のとき，禮円の直径はいかほどか。

仁円の半径と中心座標を R, (0, 0)
義円の半径と中心座標を r1, (0, R + r1), (0, R - r1)
禮円の半径と中心座標を r2, (x21, y21), (x22, y21)
斜線と人円の交点座標を (x, y)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x::positive, y::negative, R::positive, r1::positive, r2::positive,
     x21::positive, y21::positive, x22::positive, y22::positive
eq1 = dist2(0, R + 2r1, x, y,0, R - r1, r1)
eq2 = dist2(0, R + 2r1, x, y, x21, y21, r2)
eq3 = dist2(0, R + 2r1, x, y, x22, y22, r2)
eq4 = x21^2 + y21^2 - (R - r2)^2
eq5 = x22^2 + (y22 - R - r1)^2 - (r1 - r2)^2
eq6 = y21/x21 - x/(R + 2r1 - y)
eq7 = (y22 - R - r1)/x22 - x/(R + 2r1 - y)
eq8 = x^2 + y^2 - R^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 Float64.(v), r.f_converged
end;

function H(u)
   (R, r2, x21, y21, x22, y22, x, y) = u
   return [
       -R^2*r1^2 - 4*R*r1^3 + 2*R*r1^2*y - 4*r1^4 + 4*r1^3*y + 8*r1^2*x^2 - r1^2*y^2,  # eq1
       -R^2*r2^2 + R^2*x^2 - 2*R^2*x*x21 + R^2*x21^2 - 4*R*r1*r2^2 + 4*R*r1*x^2 - 8*R*r1*x*x21 + 4*R*r1*x21^2 + 2*R*r2^2*y - 2*R*x^2*y21 + 2*R*x*x21*y + 2*R*x*x21*y21 - 2*R*x21^2*y - 4*r1^2*r2^2 + 4*r1^2*x^2 - 8*r1^2*x*x21 + 4*r1^2*x21^2 + 4*r1*r2^2*y - 4*r1*x^2*y21 + 4*r1*x*x21*y + 4*r1*x*x21*y21 - 4*r1*x21^2*y - r2^2*x^2 - r2^2*y^2 + x^2*y21^2 - 2*x*x21*y*y21 + x21^2*y^2,  # eq2
       -R^2*r2^2 + R^2*x^2 - 2*R^2*x*x22 + R^2*x22^2 - 4*R*r1*r2^2 + 4*R*r1*x^2 - 8*R*r1*x*x22 + 4*R*r1*x22^2 + 2*R*r2^2*y - 2*R*x^2*y22 + 2*R*x*x22*y + 2*R*x*x22*y22 - 2*R*x22^2*y - 4*r1^2*r2^2 + 4*r1^2*x^2 - 8*r1^2*x*x22 + 4*r1^2*x22^2 + 4*r1*r2^2*y - 4*r1*x^2*y22 + 4*r1*x*x22*y + 4*r1*x*x22*y22 - 4*r1*x22^2*y - r2^2*x^2 - r2^2*y^2 + x^2*y22^2 - 2*x*x22*y*y22 + x22^2*y^2,  # eq3
       x21^2 + y21^2 - (R - r2)^2,  # eq4
       x22^2 - (r1 - r2)^2 + (-R - r1 + y22)^2,  # eq5
       -x/(R + 2*r1 - y) + y21/x21,  # eq6
       -x/(R + 2*r1 - y) + (-R - r1 + y22)/x22,  # eq7
       -R^2 + x^2 + y^2,  # eq8
   ]
end;

r1 = 6/2
iniv = BigFloat[6.1, 1.1, 4.8, 1.6, 1.8, 9.6, 5.2, -2.8]
res = nls(H, ini=iniv)

   ([6.0, 1.0, 4.714045207910317, 1.6666666666666667, 1.8856180831641267, 9.666666666666666, 5.261948151328113, -2.883036880224506], true)

禮円の半径は義円の半径の 1/3 である。
義円の直径が 6 寸のとき，禮円の直径は 2 寸である。

ちなみに，仁円の直径は義円の直径の 2 倍の 12 寸である。その他のパラメータは以下のとおりである。
R = 6;  r2 = 1;  x21 = 4.71405;  y21 = 1.66667;  x22 = 1.88562;  y22 = 9.66667;  x = 5.26195;  y = -2.88304

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 6/2
   (R, r2, x21, y21, x22, y22, x, y) = res[1]
   @printf("義円の直径が %g 寸のとき，禮円の直径は %g 寸\n", 2r1, 2r2)
   @printf("R = %g;  r2 = %g;  x21 = %g;  y21 = %g;  x22 = %g;  y22 = %g;  x = %g;  y = %g\n", R, r2, x21, y21, x22, y22, x, y)
   plot()
   circle(0, 0, R, :blue)
   circle(0, R - r1, r1)
   circle(0, R + r1, r1)
   circle2(x22, y22, r2, :green)
   circle2(x21, y21, r2, :green)
   segment(0, R + 2r1, x, y, :magenta)
   segment(0, R + 2r1, -x, y, :magenta)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R + 2r1, "R+2r1", :red, :center, :bottom, delta=delta/2)
       point(x, y, "(x,y) ", :blue, :right, :vcenter)
       point(0, 0, "仁円:R,(0,0)", :blue, :center, delta=-delta)
       point(0, R + r1, "義円:r1,(0,R+r1)", :red, :center, delta=-1.5delta)
       point(0, R - r1, "義円:r1,(0,R-r1)", :red, :center, delta=-delta)
       point(x21, y21, "禮円:r2,(x21,y21) ", :green, :right, :vcenter)
       point(x22, y22, " 禮円:r2,(x22,y22) ", :green, :left, :bottom)
   end
end;

算額（その894）

算額（その894）

七二 加須市大字外野 棘脱地蔵堂 明治7年(1874)

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円内に 8 個の円を入れる。外円の直径が 55 寸のとき，甲円の直径はいかほどか。

注：外円内には 8 個の円が入っているとしか言っていないが，上下の 2 個（緑，甲円と呼んでいるもの）と左右の 4 個（マゼンタ，多分乙円だろうが）は大きさが違う。青円には名前がついていないが，大円と呼んでおく。

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

include("julia-source.txt");

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

   (R*(-1 + sqrt(5))/2, R*(3 - sqrt(5))/2, R*(-1 + sqrt(5))/4, R*(-1/2 + sqrt(5)/2))

甲円（上下の小さい方の 2️ 円）の半径 r2 は，外円の半径 R の (3 - √5)/2 倍である。
外円の直径が 55 寸のとき，甲円の直径は 55*(3 - √5)/2 = 21.00813061875578 寸である。

算額では「答曰甲円貳拾四寸有奇」とあるが，「貳拾四寸」は「貳拾一寸」の誤記である。
術は「置五个開平方以減三余乗半」つまり，「5 の平方根から 3 を引いて，0.5 を掛ける (3 - √5)/2」といっているので正しい。

ちなみに，大円，乙円の直径は 33.99186938124422 寸，16.99593469062211 寸である。

外円の直径が 832040 寸のとき，甲円の直径は 317811.00000053743 である。

function draw(more=false)
   pyplot(size=(600, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 55//2
   (r1, r2, r3, x3) = (R*(-1 + sqrt(5))/2, R*(3 - sqrt(5))/2, R*(-1 + sqrt(5))/4, R*(-1/2 + sqrt(5)/2))
   @printf("大円の直径 = %g; 甲円の直径 = %g;  乙円の直径 = %g\n", 2r1, 2r2, 2r3)
   plot()
   circle(0, 0, R)
   circle22(0, R - r1, r1, :blue)
   circle22(0, R - r2, r2, :green)
   circle4(x3, r3, r3, :magenta)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - r1, "大円:r1,(0,R-r1)", :blue, :center, delta=-delta/2)
       point(0, R - r2, "甲円:r2,(0,R-r2)", :green, :center, delta=-delta/2)
       point(x3, r3, "乙円:r3\n(x3,r3)", :magenta, :left, delta=-delta/2)
   end
end;