算額（その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;

算額（その893）

算額（その893）

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

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

大きさの同じ正三角形 2 個に挟まれて甲円と乙円がある。乙円は両方の正三角形の斜辺に外接し，甲円と内接している。

条件式は 1 個だけであり，正三角形の大きさには依存しない。
甲円と乙円の半径を r1, r2 として以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive
eq1 = r2/(2r1 - r2) - 1//2
res = solve(eq1, r2)[1]
res |> println

   2*r1/3

乙円の半径は甲円の半径の 2/3 倍である。
甲円の直径が 3 寸のと，乙円の直径は 2 寸である。

function draw(more=false)
   pyplot(size=(600, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 2
   r1 = 3//2
   r2 = 2r1/3
   plot([2a, a, 0, -a, -2a, 2a], [0, √3a, 0, √3a, 0, 0], color=:blue, lw=0.5)
   circle(0, r1, r1)
   circle(0, 2r1 - r2, r2, :green)
   circle(0, 0, 0.2r1, :blue, beginangle=0, endangle=60)
   circle(0, 0, 0.23r1, :blue, beginangle=0, endangle=60)
   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, r1, "甲円:r1,(0,r1)", :red, :center, :bottom, delta=delta/2)
       point(0, 2r1 - r2, "乙円:r2,(0,2r1-r2)", :green, :center, :bottom, delta=2delta)
       point(0.25r1, 0.2r1, "60°", :blue, :left, :vcenter, mark=false)
   end
end;

算額（その892）

算額（その892）

七十 加須市大字外野 棘脱地蔵堂 明治6年(1873)

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

直角三角形内に全円 1 個，等円 2 個，両方の等円に接する斜線を入れる。全円，等円の直径がそれぞれ 2 寸，1.2 寸のとき，鈎（図のように置いたときの直角三角形の高さ）を求めよ。

図のように置いたときの直角三角形の底辺の長さ，高さを a, b
斜線と底辺の交点座標を (c, 0)
全円の半径と中心座標を r1, (r1, r1)
等円の半径と中心座標を r1, (r2, r2), (x2, r2)
とおき，以下の連立方程式を解く。

1.　数値解を求める

条件を記述するときに r1, r2 の数値を代入した式にすれば，4 元連立方程式は解ける。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, 弦::positive,
     r1::positive, r2::positive, x2::positive
(r1, r2) = (2//2, 12//20)
eq1 = b + c - sqrt(b^2 + c^2) - 2r2
eq2 = dist2(a, 0, 0, b, r1, r1, r1)
eq3 = dist2(a, 0, 0, b, x2, r2, r2)
eq4 = dist2(c, 0, 0, b, x2, r2, r2);
res = solve([eq1, eq2, eq3, eq4], (a, b, c, x2));
res[1]  # 1 of 3

   (13/4, 18/5, 3/2, 19/10)

鈎は 18/5 = 3.6 寸である。「答」では 3.7 寸になっているが，誤記である。
また，一般的には直角三角形の高さ（鈎）は底辺（股）より短いことが多いが，この図形においていは鈎のほうが長い。そして，算額の図は鈎が股より短く書かれているのでミスリーディングである。

2.　数式解を求める

18/5 がなぜ出てくるかは，r1, r2 を記号のままで 4 元連立方程式を解かなければならない。しかし，SymPy では自動的には解けないので，以下のように手動で解く。

方針としては，それぞれの方程式を眺めて，途中の式が複雑にならないものから変数を消去していき，最終的に b, r1, r2 だけを含む式を求める。

using SymPy
@syms a::positive, b::positive, c::positive, 弦::positive,
     r1::positive, r2::positive, x2::positive
eq1 = b + c - sqrt(b^2 + c^2) - 2r2
eq2 = dist2(a, 0, 0, b, r1, r1, r1)
eq3 = dist2(a, 0, 0, b, x2, r2, r2)
eq4 = dist2(c, 0, 0, b, x2, r2, r2);
eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println

   b + c - 2*r2 - sqrt(b^2 + c^2)
   a^2*b^2 - 2*a^2*b*r1 - 2*a*b^2*r1 + 2*a*b*r1^2
   a^2*b^2 - 2*a^2*b*r2 - 2*a*b^2*x2 + 2*a*b*r2*x2 - b^2*r2^2 + b^2*x2^2
   b^2*c^2 - 2*b^2*c*x2 - b^2*r2^2 + b^2*x2^2 - 2*b*c^2*r2 + 2*b*c*r2*x2

eq3, eq4 を対象にして x2 を消去する。

res_x2 = solve(eq4, x2)[2]  # 2 of 2
res_x2 |> println

   (c*(b - r2) + r2*sqrt(b^2 + c^2))/b

eq4 にはなかった sqrt(b^2 + c^2) が出てくる。これは，eq1 に出てくるので，b + c - 2r2 に置き換える。

res_x2 = res_x2(sqrt(b^2 + c^2) => (b + c - 2r2)) |> factor
res_x2 |> println

   -(-b*c - b*r2 + 2*r2^2)/b

eq3, eq4 の x2 に res_x2 を代入する。

eq13 = eq3(x2 => res_x2) |> simplify

   a^2*b^2 - 2*a^2*b*r2 - 2*a*b*(b*c + b*r2 - 2*r2^2) + 2*a*r2*(b*c + b*r2 - 2*r2^2) - b^2*r2^2 + (b*c + b*r2 - 2*r2^2)^2

eq14 = eq4(x2 => res_x2)/2r2^2 |> simplify
eq14 |> println

   b*c - 2*b*r2 - 2*c*r2 + 2*r2^2

eq14 を解いて c を求め，eq13 の c に代入する。代入後の eq13 = 0 は，eq13 の分子が 2 式の積になるので，どちらかが 0 であればよい。

res_c = solve(eq14, c)[1]
res_c |> println

   2*r2*(b - r2)/(b - 2*r2)

eq13(c => res_c) |> simplify |> numerator |> factor |> println

   (a*b - 2*a*r2 - 2*b*r2 + 2*r2^2)*(a*b^3 - 4*a*b^2*r2 + 4*a*b*r2^2 - 4*b^3*r2 + 12*b^2*r2^2 - 16*b*r2^3 + 8*r2^4)

t1 ≠ 0，t2 = 0 である。

t1 = a*b - 2*a*r2 - 2*b*r2 + 2*r2^2
t2 = a*b^3 - 4*a*b^2*r2 + 4*a*b*r2^2 - 4*b^3*r2 + 12*b^2*r2^2 - 16*b*r2^3 + 8*r2^4;

eq2 は a, b, r1 を含むので，a についてとき，これを t2 の a に代入する。

res_a = solve(eq2, a)[1]
res_a |> println

   2*r1*(b - r1)/(b - 2*r1)

t12 = t2(a => res_a) |> simplify
t12 |> println

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

これを b について求めると r1, r2 だけを含む式（解）になる。

res_b = solve(t12, b)[1]  # 1 of 4
res_b |> println

   -2*r2^2/(r1 - 2*r2)

これは「術曰置等円径自之等円径二段全円径差以除得鈎合問」すなわち「等円の直径を二乗したものを全円の直径と等円の直径の二倍との差で割る」
(2r1)^2/(2r1 - 2(2r2)) = 2r1^2/(r1 - 2r2) である。

r1 = 1, r2 = 0.6 を代入すると鈎(b) が得られる。すなわち，鈎は 3.6 寸である。

res_b(r1 => 1, r2 => 0.6) |> println

   3.60000000000000

以下，逆順に代入していくことで，すべての解を得ることができる。

res_a(r1 => 1, r2 => 0.6, b => 3.6) |> println

   3.25000000000000

res_c(r1 => 1, r2 => 0.6, b => 3.6, a => 3.25) |> println

   1.50000000000000

res_x2(r1 => 1, r2 => 0.6, b => 3.6, a => 3.25, c => 1.5) |> println

   1.90000000000000

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (1, 0.6)
   (a, b, c, x2) =  (13/4, 18/5, 3/2, 19/10)
   plot([0, a, 0, 0], [0, 0, b, 0], color=:blue, lw=0.5)
   circle(r1, r1, r1)
   circle(r2, r2, r2, :green)
   circle(x2, r2, r2, :green)
   segment(c, 0, 0, b, :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, b, " b(鈎)", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a(股)", :blue, :left, :bottom, delta=delta/2)
       point(c, 0, "c", :magenta, :left, :bottom, delta=delta/2)
       point(r1, r1, "全円:r1,(r1,r1)", :red, :right, delta=-delta/2)
       point(r2, r2, "等円:r2,(r2,r2)", :green, :center, delta=-delta/2)
       point(x2, r2, "等円:r2,(x2,r2)", :green, :center, delta=-delta/2)
   end
end;

 

算額（その891）

算額（その891）

七十 加須市大字外野 棘脱地蔵堂 明治6年(1873)

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

外円の中に円弧 2 本，甲円，乙円各 2 個を入れる。乙円の径が□□甲円の直径はいかほどか。

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

include("julia-source.txt");

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

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (33*r2/4, 33*r2/8, 5*r2, 3*r2)

甲円の半径は乙円の半径の 33/8 倍である。

算額では乙円の直径は欠損文字であるが，「答曰甲円径三十二寸」とあるので，甲円の直径が 32 寸になるときの乙円の直径を逆算すると 256/33 = 7.75757757575 である。実に半端な数である。答が「甲円径三十三寸」の誤記ならば乙円の直径が 8 寸なので，丸く収まるのだが。
術も妥当なものではない。「術曰置乙円径四十三個以除八個二得甲円径合問」とあるが，これも「置乙円径三十三個以除八個乗甲円径」の誤記なら納得はできる。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 128/33
   (R, r1, x2, y2) = r2 .* (33/4, 33/8, 5, 3)
   @printf("乙円の直径が %g のとき，甲円の直径は %g\n", 2r2, 2r1)
   plot()
   circle(0, 0, R, :blue)
   circle(0, -R, R, beginangle=30, endangle=150)
   circle(0, R, R, beginangle=210, endangle=330)
   circle42(0, r1, r1, :magenta)
   circle4(x2, y2, r2, :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)
       point(r1, 0, "甲円:r1,(r1,0)", :magenta, :center, delta=-delta/2, deltax=6delta)
       point(0, r1, "甲円:r1,(r1,0)", :magenta, :center, :bottom, delta=delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :black, :center, delta=-delta, deltax=-6delta)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その890）

算額（その890）

七十 加須市大字外野 棘脱地蔵堂 明治6年(1873)

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

外円内に円弧 6 個と等円 6 個を入れる。外円の直径が 3 寸のとき，等円の直径はいかほどか。

外円の半径を R とする。
円弧は半径 r1 の円の一部で，図の太線で示す円弧に 2 個の等円が接している。また，等円の中心は正六角形を構成する。
等円の半径を r2 として，以下の連立方程式を解く。

include("julia-source.txt");

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

   (R*(-3 + 4*sqrt(3))/6, R*(-3 + 4*sqrt(3))/13)

円弧がその一部である円の半径と等円の半径は，R(4√3 - 3)/6，R(4√3 - 3)/13 である。
外円の直径が 3 寸のとき，等円の直径は 3(4√3 - 3)/13 = 0.9065084377558866 寸である。

算額では「答曰等円径□分二厘八毛余」とあるが，欠落文字はともかくとして，正しいものではない。
また，「術曰置二十四個開平方一個減以除十三個乗外円径得等円径合問」とあるが，数式は「外円直径 × (√24 - 1)/13」なので，これも 0.899764496669159 になり，正しくもないし，「答」とすら一致しない。

以下では，図を描くために外円と円弧の交点座標 (x0, y0) を求める。

@syms x0, y0
eq3 = x0^2 + (y0 - r1)^2 - r1^2
eq4 = x0^2 + y0^2 - R^2
res2 = solve([eq3, eq4], (x0, y0))
x0 = res2[2][1] |> simplify
y0 = res2[2][2] |> simplify
(x0, y0) |> println

   (2*R*sqrt(28 - 6*sqrt(3))/13, R*(3 + 4*sqrt(3))/13)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 3//2
   (r1, r2) = R .* (4sqrt(3) - 3) .* (1/6, 1/13)
   (x0, y0) = R .* (2sqrt(28 - 6*sqrt(3)), 3 + 4sqrt(3)) ./ 13
   θ = atand(y0 - r1, x0)
   @printf("外円の直径 = %g;  等円の直径 = %g;  円弧の直径 = %g\n", 2R, 2r2, 2r1)
   plot()
   circle(0, 0, R, :blue)
   rotate(R - r2, 0, r2, angle=60)
   #circle((R - r2)/2, (R - r2)√3/2, r2)#, angle=60)
   for i = 0:5
       theta = 90 + 60i
       circle(r1*cosd(theta), r1*sind(theta), r1, :green, beginangle=180 + theta, endangle=270 + theta + θ, lw=i == 0 ? 1.5 : 0.5)
   end
   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, r1)
       point((R - r2)/2, (R - r2)√3/2, "等円:r2\n((R-r2)/2,\n(R-r2)√3/2)", :red, :center, delta=-delta/2)
       point(R - r2, 0, "等円:r2,(0,R-r2)", :red, :center, delta=-delta/2)
       point(0, r1, "r1", :green, :center, :bottom, delta=2delta)
       point(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=delta)
   end
end;

算額（その889）

算額（その889）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

外円内に正三角形と甲円を 1 個ずつ，乙円を 2 個入れる。甲円の直径と正三角形の一辺の長さが共に 1 寸のとき，乙円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive, r1::positive,
     r2::positive, x2::positive, y2::negative
(r1, a) = (1, 1) .// 2
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x2^2 + (y2 - r1 + R)^2 - (r1 + r2)^2
eq3 = (a - x2)^2 + (R - a√Sym(3) - y2)^2 - r2^2
eq4 = a√Sym(3) + 2r1 - 2R
res = solve([eq1, eq2, eq3, eq4], (R, r2, x2, y2));

res[1][1] |> simplify |> println
res[1][2] |> simplify |> println
res[1][3] |> simplify |> println
res[1][4] |> simplify |> println

   sqrt(3)/4 + 1/2
   -19*sqrt(888 + 534*sqrt(3))/2704 - 9/104 + 3*sqrt(296 + 178*sqrt(3))/338 + 25*sqrt(3)/104
   -2*sqrt(888 + 534*sqrt(3))/169 + 3*sqrt(3)/26 + 15/52 + 19*sqrt(296 + 178*sqrt(3))/676
   -sqrt(296 + 178*sqrt(3))/52 - sqrt(3)/8 + sqrt(3)*sqrt(296 + 178*sqrt(3))/208 + 3/8

かなり不明な定数が出てくるが，乙円の半径は
-2*sqrt(888 + 534*sqrt(3))/169 + 3*sqrt(3)/26 + 15/52 + 19*sqrt(296 + 178*sqrt(3))/676
である。

res[1][1].evalf() |> println
res[1][2].evalf() |> println
res[1][3].evalf() |> println
res[1][4].evalf() |> println

   0.933012701892219
   0.248826680729423
   0.675359395597227
   -0.109545416760023

数値としては，乙円の半径は 0.248826680729423（直径は 0.497653361458846）である。

2*0.248826680729423

   0.497653361458846

甲円の直径，正三角形の一辺の長さが共に 1 寸のとき，乙円の直径は 0.497653361458846 寸である。

その他のパラメータは以下のとおりである。
a = 0.5;  r1 = 0.5;  R = 0.933013;  r2 = 0.248827;  x2 = 0.675359;  y2 = -0.109545

算額の答では，「乙円径五分一厘有奇」となっているがそれだとすると，外円をはみ出す（図の赤円）。得られた答えは青円で，ちゃんと外円に接している。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1) = [1, 1] .//2
   (R, r2, x2, y2) = (
       sqrt(3)/4 + 1/2,
       -19*sqrt(888 + 534*sqrt(3))/2704 - 9/104 + 3*sqrt(296 + 178*sqrt(3))/338 + 25*sqrt(3)/104,
       -2*sqrt(888 + 534*sqrt(3))/169 + 3*sqrt(3)/26 + 15/52 + 19*sqrt(296 + 178*sqrt(3))/676,
       -sqrt(296 + 178*sqrt(3))/52 - sqrt(3)/8 + sqrt(3)*sqrt(296 + 178*sqrt(3))/208 + 3/8)
   @printf("甲円の直径，正三角形の一辺の長さが共に 1 寸のとき，乙円の直径は %.15g 寸である。\n", 2r2)
   @printf("a = %g;  r1 = %g;  R = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", a, r1, R, r2, x2, y2)
   plot()
   circle(0, 0, R, :blue)
   plot!([a, 0, -a, a], [R - √3a, R, R - √3a, R - √3a], color=:green, lw=0.25)
   circle(0, r1 - R, r1)
   circle(x2, 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(a, R - √3a, "(a,R-√3a) ", :green, :right, :bottom, delta=delta/2)
       point(0, r1 - R, "甲円:r1,(0,r1-R)", :red, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :magenta, :center, delta=-delta/2)
   end
end;

算額（その888）

算額（その888）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

楕円の中に正方形 2 個と全円 1 個を入れる。正方形の一辺の長さが 1 寸のとき，全円の直径はいかほどか。
ちなみに，以下の図は正方形の一辺の長さが 123.45 のときの図である。

楕円の長半径と短半径，中心座標を a, b, (0, 0) とする。全円の半径は短半径と同じである。
正方形の一辺の長さを c とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,　c::positive
s2 = √Sym(2)
eq1 = (b + s2*c/2)^2/a^2 + (s2*c/2)^2/b^2 - 1
eq2 = b + s2*c - a
res = solve([eq1, eq2], (a, b))[1]

   (c*(1 + sqrt(2)), c)

b すなわち全円の半径は，正方形の一辺の長さ c である。
正方形の一辺の長さが c のとき，全円の直径は 2c である。
正方形の一辺の長さが 1 寸のとき，全円の直径は 2 寸である。

function diamond(a, c)
   plot!([a, a - c/√2, a - √2c, a - c/√2, a], [0, c/√2, 0, -c/√2, 0], color=:green, lw=0.5)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = 123.45
   b = c
   a = b + c√2
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle(0, 0, b, :blue)
   diamond(a, c)
   diamond(-b, c)
   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, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :red, :left, :bottom, delta=delta/2)
       point(b + √2c/2, 0, "b+√2c/2", :green, :center, :bottom, delta=delta/2)
   end
end;

算額（その887）

算額（その887）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

長方形内に菱形 1 個，甲円 4 個，乙円 2 個，丙円 2 個を入れる。長方形の長辺と短辺が 15 寸，14 寸のとき，甲円の直径はいかほどか。

長方形の長辺と短辺を 2a, 2b
甲円の半径と中心座標を r1, (a - r2, b - r1)
乙円の半径と中心座標を r2, (0, b - r2)
丙円の半径と中心座標を r3, (a - r3, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive
@syms a, b, r1, r2, r3
eq1 = (r1 - r3)^2 + (b - r1)^2 - (r1 + r3)^2 |> expand
eq2 = (a - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq3 = dist2(a - 2r3, 0, 0, b - 2r2, a - r1, b - r1, r1);

eq1, eq2 を連立させて r2，r3 を求めると，それぞれは a と r1，b と r2 を含む式になる。

res = solve([eq1, eq2], (r2, r3))
res[r2] |> factor |> println
res[r3] |> factor |> println

   (-a + r1)^2/(4*r1)
   (-b + r1)^2/(4*r1)

r2, r3 を eq3 に代入し，因子分析する。

eq13 = eq3(r2 => (-a + r1)^2/(4*r1), r3 => (-b + r1)^2/(4*r1))
eq13 |> factor

eq13 = 0 となるのは，分子の 4 項のいずれかが 0 になればよい。

eq13 |> factor |> numerator |> println

   (-a*b - a*r1 - b*r1 + r1^2)*(-a*b + a*r1 + b*r1 + r1^2)*(a*b - 3*a*r1 - b*r1 + r1^2)*(a*b - a*r1 - 3*b*r1 + r1^2)

それぞれの項 = 0 を解いて r1 を求め，図形の条件を満たすものを探索すると，
-a*b + a*r1 + b*r1 + r1^2 = 0 を解いた 2 つの解のうちの r1 = -(a + b - sqrt(a^2 + 6*a*b + b^2))/2 が適解である（ほかの 7 個の解はすべて不適）。

res2 = solve(-a*b + a*r1 + b*r1 + r1^2, r1);
res2[2] |> factor |> println
res2[2] |> factor |> display

   -(a + b - sqrt(a^2 + 6*a*b + b^2))/2

r1 = (sqrt(a^2 + 6*a*b + b^2) - a - b)/2 を (r1 - a)^2/4r1，(r1 - b)^2/4r1 に代入すれば r2, r3 が求まる。
数式は複雑で SymPy では簡約化できないが，数値としては簡単に求まる。

res_r2 = res[r2](r1 => (sqrt(a^2 + 6*a*b + b^2) - a - b)/2) |> factor |> display
res_r3 = res[r3](r1 => (sqrt(a^2 + 6*a*b + b^2) - a - b)/2) |> factor |> display

まとめると，r1, r2, r3 は a, b を用いて以下のように計算される。
   r1 = (sqrt(a^2 + 6a*b + b^2) - a - b)/2
   r2 = (a - r1)^2/4r1
   r3 = (b - r1)^2/4r1

長方形の長辺 2a，短辺 2b が 15, 14 のとき，甲円の直径は 6 である。

なお，計算された r1 の値が　2r1 > a または 2r1 > b になると不適切な図になる。

a = 7.5
b = 7
sqrt(a^2 + 6a*b + b^2) - a - b |> println

   6.0

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # (a, b) から (r1, r2, r3) を計算する
   r1 = (sqrt(a^2 + 6a*b + b^2) - a - b)/2
   if 2r1 > a || 2r1 > b
       println("不適切な図になります")
       return
   end
   r2 = (a - r1)^2/4r1
   r3 = (b - r1)^2/4r1
   @printf("長方形の長辺，短辺が %g, %g のとき，甲円の直径は %g である\n", 2a, 2b, 2r1)
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g\n", a, b, r1, r2, r3)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   plot!([0, a - 2r3, 0, 2r3 - a, 0], [2r2 - b, 0, b - 2r2, 0, 2r2 - b], color=:orange, lw=0.5)
   circle4(a - r1, b - r1, r1, :green)
   circle22(0, b - r2, r2)
   circle2(a - r3, 0, 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(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a - r1, b - r1, "甲円:r1,(a-r1,b-r1)", :green, :center, delta=-delta/2)
       point(0, b - r2, "乙円:r2\n(0,b-r2)", :red, :center, delta=-delta/2)
       point(a - r3, 0, "丙円:r3\n(a-r3,0)", :magenta, :center, delta=-delta/2)
   end
end;

   長方形の長辺，短辺が 24.69, 20.222 のとき，甲円の直径は 9.22291 である
   a = 12.345;  b = 10.111;  r1 = 4.61145;  r2 = 3.24235;  r3 = 1.63967

算額（その886）

算額（その886）

六十四 加須市不動岡 総願寺 慶応二丙寅(1866)

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

外円の中に 2 本の斜線を引き大円 2 個，小円 2 個を入れる。大円の直径が 1 寸のとき，小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, -r1), (0, r1)
小円の半径と中心座標を r2, (x2, y2); y2 ≠ 0 小円の中心は x 軸上にはない
外円と斜線の交点座標を (x, y); y < 0

include("julia-source.txt");

using SymPy
@syms d, R::positive, r1::positive,
     r2::positive, x2::positive, y2, x::positive, y::negative
R = 2r1
eq1 = dist2(0, R, x, y, 0, -r1, r1)
eq2 = dist2(0, R, x, y, x2, y2, r2)
eq3 = x2^2 + y2^2 - (R - r2)^2
eq4 = x2^2 + (r1 - y2)^2 - (r1 + r2)^2
eq5 = R^2 - (x^2 + y^2)|> expand;
solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, x, y))[1]  # 1 of 8

   (16*r1/25, 24*sqrt(2)*r1/25, 2*r1/25, 8*sqrt(2)*r1/9, -14*r1/9)

8 組の解が得られるが，最初のものが適解である。

小円の半径は大円の半径の 16/25 倍である。
大円の直径が 1 寸のとき，小円の直径は 16/25 = 0.64 寸 = 6 分 5 厘である。

その他のパラメータは以下のとおりである。
r1 = 0.5;  R = 1;  r2 = 0.32;  x2 = 0.678823;  y2 = 0.04; x = 0.628539;  y = -0.777778

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   # (r2, x2, y2, x, y) = (16*r1/25, 24*sqrt(2)*r1/25, 2*r1/25, 8*sqrt(2)*r1/9, -14*r1/9)
   (r2, x2, y2, x, y) = r1 .* (16/25, 24√2/25, 2/25, 8√2/9, -14/9)
   R = 2r1
   @printf("大円の直径が %g のとき，小円の直径は %g である\n", 2r1, 2r2)
   @printf("r1 = %g;  R = %g;  r2 = %g;  x2 = %g;  y2 = %g; x = %g;  y = %g\n", r1, R, r2, x2, y2, x, y)
   plot()
   circle(0, 0, R, :blue)
   circle22(0, r1, r1)
   circle2(x2, y2, r2, :green)
   segment(0, R, x, y, :magenta)
   segment(0, R, -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, " R", :blue, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, "大円:r1\n(0,r1)", :red, :center, delta=-delta/2)
       point(0, -r1, "大円:r1\n(0,-r1)", :red, :center, delta=-delta/2)
       point(0, 0, "", :blue)
       point(x2, y2, "小円:r2\n(x2,y2)", :green, :center, delta=-delta)
   end
end;