裏 RjpWiki

算額（その865）

岩手県平泉町　中尊寺阿弥陀堂（中尊寺地蔵院にて保管） 安政6年(1859)
http://www.wasan.jp/iwate/chusonji3.html

牧下英世：数学史を取り入れた授業実践―算額の教材化と総合的な学習―，2000筑波大学附属駒場論集第40集
https://tsukuba.repo.nii.ac.jp/record/6486/files/10.pdf

外円内に，水平な弦を 1 本と斜線を 2 本描き，区画された領域に大円 1 個と等円 3 個を入れる。外円と大円の直径がわかったとき，等円の直径はいかほどか。

弦と y 軸の交点座標を (0, y)
外円と斜線の交点を (x00, y00), (01, y01)
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
等円の半径と中心座標を r2, (0, R - r2), (x2, y + r2)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     x2::positive, y::positive, x00::positive,
     y00::positice,  x01::positive, y01::positive, d
y = 2r1 - R
eq1 = dist(x01, y01, -x00, y00, x2, y + r2) - r2^2
eq2 = dist(-x01, y01, x00, y00, 0, R - r2) - r2^2
eq3 = dist(-x01, y01, x00, y00, 0, r1 - R) - r1^2
eq4 = x2^2 + (y + r2)^2 - (R - r2)^2
eq5 = (x01^2 + y01^2) - R^2
eq6 = (x00^2 + y00^2) - R^2
eq1 = numerator(apart(eq1, d)) |> simplify
eq2 = numerator(apart(eq2, d)) |> simplify
eq3 = numerator(apart(eq3, d)) |> simplify;

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)
   (r2, x2, x00, y00, x01, y01) = u
   return [
R^2*x00^2 + 2*R^2*x00*x01 + R^2*x01^2 - 4*R*r1*x00^2 - 8*R*r1*x00*x01 - 4*R*r1*x01^2 - 2*R*r2*x00^2 - 4*R*r2*x00*x01 - 2*R*r2*x01^2 + 2*R*x00^2*y01 + 2*R*x00*x01*y00 + 2*R*x00*x01*y01 - 2*R*x00*x2*y00 + 2*R*x00*x2*y01 + 2*R*x01^2*y00 - 2*R*x01*x2*y00 + 2*R*x01*x2*y01 + 4*r1^2*x00^2 + 8*r1^2*x00*x01 + 4*r1^2*x01^2 + 4*r1*r2*x00^2 + 8*r1*r2*x00*x01 + 4*r1*r2*x01^2 - 4*r1*x00^2*y01 - 4*r1*x00*x01*y00 - 4*r1*x00*x01*y01 + 4*r1*x00*x2*y00 - 4*r1*x00*x2*y01 - 4*r1*x01^2*y00 + 4*r1*x01*x2*y00 - 4*r1*x01*x2*y01 - r2^2*y00^2 + 2*r2^2*y00*y01 - r2^2*y01^2 - 2*r2*x00^2*y01 - 2*r2*x00*x01*y00 - 2*r2*x00*x01*y01 + 2*r2*x00*x2*y00 - 2*r2*x00*x2*y01 - 2*r2*x01^2*y00 + 2*r2*x01*x2*y00 - 2*r2*x01*x2*y01 + x00^2*y01^2 + 2*x00*x01*y00*y01 - 2*x00*x2*y00*y01 + 2*x00*x2*y01^2 + x01^2*y00^2 - 2*x01*x2*y00^2 + 2*x01*x2*y00*y01 + x2^2*y00^2 - 2*x2^2*y00*y01 + x2^2*y01^2,  # eq1
R^2*x00^2 + 2*R^2*x00*x01 + R^2*x01^2 - 2*R*r2*x00^2 - 4*R*r2*x00*x01 - 2*R*r2*x01^2 - 2*R*x00^2*y01 - 2*R*x00*x01*y00 - 2*R*x00*x01*y01 - 2*R*x01^2*y00 - r2^2*y00^2 + 2*r2^2*y00*y01 - r2^2*y01^2 + 2*r2*x00^2*y01 + 2*r2*x00*x01*y00 + 2*r2*x00*x01*y01 + 2*r2*x01^2*y00 + x00^2*y01^2 + 2*x00*x01*y00*y01 + x01^2*y00^2,  # eq2
R^2*x00^2 + 2*R^2*x00*x01 + R^2*x01^2 - 2*R*r1*x00^2 - 4*R*r1*x00*x01 - 2*R*r1*x01^2 + 2*R*x00^2*y01 + 2*R*x00*x01*y00 + 2*R*x00*x01*y01 + 2*R*x01^2*y00 - r1^2*y00^2 + 2*r1^2*y00*y01 - r1^2*y01^2 - 2*r1*x00^2*y01 - 2*r1*x00*x01*y00 - 2*r1*x00*x01*y01 - 2*r1*x01^2*y00 + x00^2*y01^2 + 2*x00*x01*y00*y01 + x01^2*y00^2,  # eq3
x2^2 - (R - r2)^2 + (-R + 2*r1 + r2)^2,  # eq4
-R^2 + x01^2 + y01^2,  # eq5
-R^2 + x00^2 + y00^2,  # eq6
   ]
end;
R = 10
r1 = 7
iniv = R .* BigFloat[0.24326219426446974, 0.5286032540384469, 0.9666966139569063, -0.2559250995198634, 0.5325907580969444, 0.8463728991347266]
res = nls(H, ini=iniv)

   ([2.1106684356050365, 4.9901186161311815, 9.935099850667232, -1.137449320748823, 5.01652946256841, 8.65068969222588], true)

外円の直径が 20，大円の直径が 14 のとき，等円の直径は 4.221336871210073 である。

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

   R = 10;  r1 = 7;  y = 4;  r2 = 2.11067;  x2 = 4.99012;  x00 = 9.9351;  y00 = -1.13745;  x01 = 5.01653;  y01 = 8.65069

算額に書かれている術は次のとおりである。注
術曰置外圓径大圓径以除之開平方加一個以除外圓径二段内滅大圓径除得等圓径合問

現代文で書くと，以下のようになろう。
外円の直径を大円の直径で割り，この平方根を求め，1 を加えたもので，外円の直径を 2 倍したものを割り，大円の直径を引くと等円の直径が得られる，これは問に合っている。

式で書くと，外円，大円の直径をそれぞれ「外」，「大」とすれば，2外/(sqrt(外/大) + 1) - 大 で，外円の直径が 20，大円の直径が 14 のとき，等円の直径が 4.22133687121007 であるとわかる。
これは，上に述べた数値解と一致する。

注：牧下は術の前半 17 文字を抜いてるので，意味をなさなくなっている。wasan.jp の算額の写真は不鮮明であるがここに書いたものと相違ない（数値解と一致するので間違いない）。

外 = 20
大 = 14
2外/(sqrt(外/大) + 1) - 大 |> println

   4.22133687121007

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # R = 1
   # r1 = 0.65
   y = 2r1 - R
   (r2, x2, x00, y00, x01, y01) = [28, 60, 110, -32, 60, 97]
   (r2, x2, x00, y00, x01, y01) = res[1]
   println("外円の直径が $(2R)，大円の直径が $(2r1) のとき，等円の直径は $(2r2) である。")
   @printf("R = %g;  r1 = %g;  y = %g;  r2 = %g;  x2 = %g;  x00 = %g;  y00 = %g;  x01 = %g;  y01 = %g\n", R, r1, y, r2, x2, x00, y00, x01, y01)
   plot()
   circle(0, 0, R)
   circle(0, r1 - R, r1, :blue)
   circle(0, R - r2, r2, :green)
   circle2(x2, y + r2, r2, :green)
   x0 = sqrt(R^2 - y^2)
   segment(-x0, y, x0, y, :magenta)
   segment(-x01, y01, x00, y00)
   segment(x01, y01, -x00, y00)
   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 - r2, "等円:r2\n(0,R-r2)", :green, :center, delta=-delta/2)
       point(x2, y + r2, "等円:r2\n(x2,y+r2)", :green, :center, delta=-delta/2)
       point(0, r1 - R, "大円:r1\n(0,r1-R)", :blue, :center, delta=-delta/2)
       point(x00, y00, "(x00,y00)", :red, :right, delta=-delta/2, deltax=-delta/2)
       point(x01, y01, "(x01,y01)", :red, :left, :bottom, delta=delta/2)
   end
end;

算額（その864）

岩手県平泉町　中尊寺阿弥陀堂（中尊寺地蔵院にて保管） 安政6年(1859)
http://www.wasan.jp/iwate/chusonji3.html

牧下英世：数学史を取り入れた授業実践―算額の教材化と総合的な学習―，2000筑波大学附属駒場論集第40集
https://tsukuba.repo.nii.ac.jp/record/6486/files/10.pdf

正方形内に正三角形と斜線，甲円と乙円が入っている。甲円と乙円の直径の差がわかっているとき，正方形の一辺の長さを求めよ。

正方形の一辺の長さを 2a
甲円の半径と中心座標を r1, (2a - r1, y1)
乙円の半径と中心座標を r2, (r2, y2); r1 > r2
直径の差を K; r1 = r2 + K/2
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, y1::positive,
     r2::positive, y2::positive, K::positive, d
s3 = √Sym(3)
r1 = r2 + K/2
eq1 = dist(2a, 0, a, s3*a, 2a - r1, y1) - r1^2
eq2 = dist(0, 0, a, s3*a, r2, y2) - r2^2
eq3 = dist(2a, 2a, a, s3*a, 2a - r1, y1) - r1^2
eq4 = dist(2a, 2a, a, s3*a, r2, y2) - r2^2
eq1 = numerator(apart(eq1, d)) |> simplify
eq2 = numerator(apart(eq2, d)) |> simplify
eq3 = numerator(apart(eq3, d)) |> simplify
eq4 = numerator(apart(eq4, d)) |> simplify;
# res = solve([eq1, eq2, eq3, eq4], (a, y1, r2, y2))

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")

   -K^2/16 - K*r2/4 - sqrt(3)*K*y1/4 - r2^2/4 - sqrt(3)*r2*y1/2 + y1^2/4,  # eq1
   -r2^2/4 - sqrt(3)*r2*y2/2 + y2^2/4,  # eq2
   -K^2/8 - sqrt(3)*K^2/16 - K*a/2 - K*r2/2 - sqrt(3)*K*r2/4 + K*y1/4 + sqrt(3)*a^2 + 2*a^2 - a*r2 - 2*a*y1 - sqrt(3)*a*y1 - r2^2/2 - sqrt(3)*r2^2/4 + r2*y1/2 + sqrt(3)*y1^2/4 + y1^2/2,  # eq3
   2*a^2 - a*r2 + sqrt(3)*a*r2 - sqrt(3)*a*y2 - a*y2 - r2^2/2 - sqrt(3)*r2^2/4 - r2*y2/2 + sqrt(3)*y2^2/4 + y2^2/2,  # eq4

eq2, eq1 を解いて，r2, y1 を求める。

res1 = solve(eq2, r2)[1]  # r2
res1 |> println

   y2*(2 - sqrt(3))

res2 = solve(eq1, y1)[2] |> factor # y1 2番目が適解
res2 |> println

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

eq3, eq4 に r2, y1 を代入して，変数消去する。

eq13 = eq3(y1 => res2, r2 => res1) |> simplify
eq13 |> println

   sqrt(3)*K^2 + 7*K^2/4 - 4*K*a - 2*sqrt(3)*K*a + sqrt(3)*K*y2 + 2*K*y2 + sqrt(3)*a^2 + 2*a^2 - 4*a*y2 + y2^2

eq14 = eq4(y1 => res2, r2 => res1) |> simplify
eq14 |> println

   2*a^2 - 6*a*y2 + 2*sqrt(3)*a*y2 - y2^2 + sqrt(3)*y2^2

eq14 から y2 を求めて，それを eq13 に代入して変数消去する。

res3 = solve(eq14, y2)[2] |> simplify # y2 2番目が適解
res3 |> println

   a*(-sqrt(42 - 24*sqrt(3)) - sqrt(14 - 8*sqrt(3)) + 2*sqrt(3))/2

eq23 = eq13(y2 => res3) |> simplify
eq23 |> println

   sqrt(3)*K^2 + 7*K^2/4 - K*a - 3*K*a*sqrt(42 - 24*sqrt(3))/2 - 5*K*a*sqrt(14 - 8*sqrt(3))/2 - 4*sqrt(3)*a^2 - a^2*sqrt(14 - 8*sqrt(3)) + a^2*sqrt(42 - 24*sqrt(3)) + 7*a^2

eq23 を a について求める。K だけを含む式になる。

res99 = solve(eq23, a)[1] |> simplify  # a 1番目が適解
res99 |> println

   K*(5*sqrt(14 - 8*sqrt(3)) + 3*sqrt(42 - 24*sqrt(3)) + 2 + 2*sqrt(sqrt(3) + 2))/(4*(-4*sqrt(3) - sqrt(2)*sqrt(7 - 4*sqrt(3)) + sqrt(6)*sqrt(7 - 4*sqrt(3)) + 7))

res99 には K が含まれているので simplify できないが，K にかかる係数は単に数値演算なので，簡約化できる。a = K*(4√3 + 7)/2 である。

@syms d
res999 = apart(res99/K, d)*K |> simplify  # a
res999 |> println

   K*(4*sqrt(3) + 7)/2

K に具体的な数値を代入すれば a の値が求まる。

res_a = res999(K => 5042).evalf()  # a
res_a |> println

   35113.0003435246

この後，逆順に代入していくことにより各変数の値を求めることができる。

res_y2 = res3(K => 5042, a => res_a).evalf()  # y2
res_y2 |> println

   42641.6741619770

res_r2 = res1(K => 5042, y2 => res_y2).evalf()  # r2
res_r2 |> println

   11425.8021556128

res_y1 = res2(K => 5042, r2 => res_r2).evalf()  # y1
res_y1 |> println

   52050.1742478581

甲円と乙円の直径の差が 5042 のとき，正方形の一辺の長さは 35113.0003435246 である。

その他のパラメータは以下のとおりである。
K = 5042; a = 35113;  r1 = 13946.802;  y1 = 52050.174;  r2 = 11425.802;  y2 = 42641.674

以下により，一挙に数値解を求めることもできる。

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)
   (a, y1, r2, y2) = u
   return [
       -(K/2 + r2)^2 + (y1 - sqrt(3)*(sqrt(3)*a*y1 - a*(-K/2 - r2))/(4*a))^2 + (-K/2 - r2 - (-sqrt(3)*a*y1/2 + a*(-K/2 - r2)/2)/(2*a))^2,  # eq1
       -r2^2 + (r2 - (a*r2/2 + sqrt(3)*a*y2/2)/(2*a))^2 + (y2 - sqrt(3)*(a*r2 + sqrt(3)*a*y2)/(4*a))^2,  # eq2
       -(K/2 + r2)^2 + (-K/2 + a*(-a*(-K/2 - r2) + (-2*a + y1)*(-2*a + sqrt(3)*a))/(a^2 + (-2*a + sqrt(3)*a)^2) - r2)^2 + (-2*a + y1 - (-2*a + sqrt(3)*a)*(-a*(-K/2 - r2) + (-2*a + y1)*(-2*a + sqrt(3)*a))/(a^2 + (-2*a + sqrt(3)*a)^2))^2,  # eq3
       -r2^2 + (-2*a + y2 - (-2*a + sqrt(3)*a)*(-a*(-2*a + r2) + (-2*a + y2)*(-2*a + sqrt(3)*a))/(a^2 + (-2*a + sqrt(3)*a)^2))^2 + (-2*a + a*(-a*(-2*a + r2) + (-2*a + y2)*(-2*a + sqrt(3)*a))/(a^2 + (-2*a + sqrt(3)*a)^2) + r2)^2,  # eq4
   ]
end;
K = 5042
iniv = BigFloat[35, 52, 11, 42].*1000
res = nls(H, ini=iniv)

   ([35113.000343524545, 52050.17424785807, 11425.802155612791, 42641.67416197693], true)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   s3 = √3
   K = 5042
   a = K*(4s3 + 7)/2
   y2 = a*(2s3 - sqrt(42 - 24s3) - sqrt(14 - 8s3))/2
   r2 = y2*(2 - s3)
   y1 = (s3 + 2)*(K + 2*r2)/2
   r1 = r2 + K/2
   @printf("K = %.8g; a = %.8g;  r1 = %.8g;  y1 = %.8g;  r2 = %.8g;  y2 = %.8g\n", K, a, r1, y1, r2, y2)
   plot([0, 2a, 2a, 0, 0], [0, 0, 2a, 2a, 0], color=:black, lw=0.5)
   plot!([0, 2a, a, 0], [0, 0, √3a, 0], color=:green, lw=0.5)
   circle(2a - r1, y1, r1)
   circle(r2, y2, r2, :blue)
   segment(0, 2a - 2(2a - √3a), 2a, 2a, :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(2a - r1, y1, "甲円:r1,(2a-r1,y1)", :red, :center, delta=-delta/2)
       point(r2, y2, "乙円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(a, √3a, "(a,√3a)", :green, :right, :bottom, delta=delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(2a, 0, " 2a", :black, :left, :bottom, delta=delta/2)
       plot!(xlims=(0, 75000))
   end
end;

算額（その863）

三十三　岩手県一関市舞川相川 菅原神社 嘉永3年(1850)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正方形内に四分円 2 個，半円 1 個，大円 1 個，小円 1 個を入れる。小円の直径が 17 寸のとき，大円の直径はいかほどか。

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

include("julia-source.txt");

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

   2-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (9*r2/4, 9*r2/17, 45*r2/68, 5*r2/4, 3*r2)
    (33*r2/4, 33*r2/17, 165*r2/68, 13*r2/4, 3*r2)

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

大円の半径は小円の半径の 33/17 倍である。
すなわち，小円の直径が 17 寸のとき，大円の直径は 33 寸である。

その他のパラメータは以下のとおりである。
r2 = 8.5;  r0 = 70.125;  r1 = 16.5;  x1 = 20.625;  x2 = 27.625;  y2 = 25.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 17/2
   (r0, r1, x1, x2, y2) = r2 .* (33/4, 33/17, 165/68, 13/4, 3)
   @printf("小円の直径が %g のとき，大円の直径は %g\n", 2r2, 2r1)
   @printf("r2 = %g;  r0 = %g;  r1 = %g;  x1 = %g;  x2 = %g;  y2 = %g\n", r2, r0, r1, x1, x2, y2)
   plot([0, r0, r0, 0, 0], [0, 0, r0, r0, 0], color=:black, lw=0.5)
   circle(0, 0, r0, beginangle=0, endangle=90)
   circle(r0, r0, r0, beginangle=180, endangle=270)
   circle(r0, r0/2, r0/2, :magenta, beginangle=90, endangle=270)
   circle(x1, r0 - x1, r1, :green)
   circle(x2, y2, r2, :orange)
   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(r0, 0, " r0", :black, :left, :bottom, delta=delta/2)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(r0, r0/2, "(r0,r0/2)", :black, :right, :bottom, delta=delta)
       point(x1, r0 - x1, "大円:r1,(x1,r0-x1)", :green, :center, delta=-delta/2)
       point(x2, y2, "小円:r2,(x2,y2)", :orange, :center, delta=-delta/2)
   end
end;

算額（その862）

二十六　岩手県一関市萩荘 赤萩観音寺前額

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

全円内に大円 1 個，中円 2 個，小円 4 個と菱形が入っている。小円の直径を知って，全円の直径を求めよ。注：「中円 4 個，小円 2 個」と書いているが明らかに誤記であろう。

大円の直径は小円の直径の (1 + √2) 倍である。これを踏まえると，大円と中円と菱形が全円の中に入っているとき，全円の直径を求めよという問題になる。それは，算額（その445）https://blog.goo.ne.jp/r-de-r/e/9f35b3c1b37302b50c61355f43956949 と同じ問題であるが，そこに示したように，描画パラメータは 5 個で，既知のパラメータはそのうちの 2 個，条件は 3 つなので，方程式は解くことができる。しかし，本問では既知の条件は1つしかないので，方程式は解けない。方程式を解くためには既知のパラメータをもう 1 つ見つける必要がある。とはいっても，もっとも単純に「大円の直径が全円の半径に等しい」と仮定する以外の妥当な方策はない。

大円の直径は小円の直径の (1 + √2) 倍なので，全円の直径は小円の直径の 2(1 + √2) ≒ 4.82842712474619 倍である。

術では「全円の直径は小円の直径の 5 倍」としているが，疑問である。

算額の答えはここまでである。

以下では，図を描くために，中円のパラメータを求める。

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

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

   r3 = 0.5;  r1 = 1.20711;  R = 2.41421;  r2 = 0.804738;  x2 = 1.60948;  y2 = 0

include("julia-source.txt");

using SymPy
@syms b::positive, R::positive, r1::positive,
     r2::positive, x2::positive, y2,
     r3::positive, d
r3 = 1//2
r1 = (1 + √Sym(2))r3
R = 2r1
b = R - r1
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq3 = dist(0, 0, sqrt(R^2 - (b - R)^2), b - R, x2, y2) - r2^2
# eq4 = 2(r1 - r3)^2 - 4r3^2
res = solve([eq1, eq2, eq3], (r2, x2, y2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (1/3 + sqrt(2)/3, sqrt(8*sqrt(2)/9 + 4/3), 0)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   r1 = (1 + √2)r3
   R = 2r1
   (r2, x2, y2) = (1/3 + sqrt(2)/3, sqrt(8*sqrt(2)/9 + 4/3), 0)
   @printf("r3 = %g;  r1 = %g;  R = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", r3, r1, R, r2, x2, y2)
   b = R - r1
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r1, r1)
   circle2(r1 - r3, R - r1, r3, :magenta)
   circle(0, R - r3, r3, :magenta)
   circle(0, R - 2r1 + r3, r3, :magenta)
   circle(x2, y2, r2, :blue)
   x = sqrt(R^2 - (b - R)^2)
   plot!([0, x, 0, -x, 0], [-R, b - R, 2b -R, b - R, -R], color=:orange, lw=0.5)
   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", :green, :left, :bottom, delta=delta/2)
       point(0, R - r1, "大円:r1,(0,R-r1)", :red, :center, delta=-delta)
       point(x2, y2, "中円:r2\n(x2,y2)", :blue, :center, delta=-delta/2)
       point(r1 - r3, R - r1, "小円:r3\n(r1-r3,R-r1)", :magenta, :center, :bottom, delta=delta)
       point(0, b - R, " b", :orange, :left, :vcenter)
   end
end;

算額（その861）

二十四　岩手県一関市萩荘字八幡 達古袋八幡神社 弘化3年

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

大円の一部である円弧内に正三角形2個と小円が入っている。正三角形の一辺の長さと小円の直径がわかっているときに大円の直径を求める術は如何に。

問では，「正三角形の一辺の長さと小円の直径がわかっているときに」とあるが，一般的にこの2つを任意に指定すると，図形を構成することは不可能である。正三角形の一辺の長さ，小円の直径，大円の直径の 3 変数のうち，1 つだけを指定して残りの 2 つのパラメータを決定するというのが妥当であろう。

よって，以下のように定めて連立方程式を解く。
大円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (0, R -r)
正三角形の一辺の長さを a とする。

eq1 は小円の中心と正三角形の一辺までの距離が小円の半径に等しいこと，
eq2 は正三角形の頂点が大円の円周上にあることを意味している。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, a::positive, x::positive
x = sqrt(R^2 - (R - 2r)^2)
eq1 = dist(x - a, R - 2r, x - a/2, R - 2r + a√Sym(3)/2, 0, R - r) - r^2
eq2 = (x - a/2)^2 + (R - 2r + a√Sym(3)/2)^2 - R^2
res = solve([eq1, eq2], (R, a));

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

大円の半径は，小円の半径の 7/3，正三角形の一辺の長さは小円の半径の √3 倍である。

「術」は不適切である。

res[2]

   (7*r/3, sqrt(3)*r)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1//2
   (R, a) = (7*r/3, sqrt(3)*r)
   x = sqrt(R^2 - (R - 2r)^2)
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r, r)
   plot!([x - a, x - a/2, x], R - 2r .+ [0, √3a/2, 0], color=:blue, lw=0.5)
   plot!(-[x - a, x - a/2, x], R - 2r .+ [0, √3a/2, 0], color=:blue, lw=0.5)
   segment(-x, R - 2r, x, R - 2r, :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", :green, :left, :bottom, delta=delta/2)
       point(0, R, " R", :green, :left, :bottom, delta=delta/2)
       point(0, R - 2r, "R-2r", :black, :center, delta=-delta/2)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(sqrt(R^2 - (R - 2r)^2), R - 2r, "(√(R^2-(R-2r)^2),R-2r)", :black, :right, delta=-delta/2)
   end
end;

算額（その860）

二十四　岩手県一関市萩荘字八幡 達古袋八幡神社 弘化3年

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

水平線の上に菱形と小円 3 個を内接する大円とそれ左側にあり大円と外接する小円がある。小円の直径が 1 寸のとき，左側の小円と大円と直線に挟まれる面積（黒積）を求めよ。

大円とそれに内接する小円 3 個，菱形 1 個は，算額（その851）にあるとおりで，小円の直径を r とすると，大円の直径は 3r である。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive, y::positive,
   x2::negative
y = r
eq1 = x^2 + y^2 - (R - r)^2
eq2 = x^2 + (R - r - y)^2 - 4r^2
eq3 = x2^2 + (r - R)^2 - (R + r)^2　
res = solve([eq1, eq2, eq3], (R, x, x2))[1]

   (3*r, sqrt(3)*r, -2*sqrt(3)*r)

小円，大円と直線により区分される面積は，台形の面積から小円の面積の 1/3 と大円の面積の 1/6 を引いたものである。
S = 2√3*r - 11π*r^2/6

@syms r, R
R = 3r
s3 = √Sym(3)
S = (r + R)*s3(R - r)/2 - PI*r^2/3 - PI*R^2/6
S |> println

   -11*pi*r^2/6 + 2*sqrt(3)*r

r = 1/2 のとき，S = 0.2921541746735554 である。

S(r => 1/2) |> N

   0.2921541746735554

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (R, x, x2) = r .* (3, √3, -2√3)
   y = r
   b = R - r
   y0 = -r
   x0 = sqrt(R^2 - y0^2)
   plot([x2, 0, 0, x2, x2], [-R, -R, 0, r - R, -R], seriestype=:shape, fillcolor=:gray30)
   circlef(0, 0, R, :orange)
   circlef(0, R - r, r, :cyan)
   circlef(x, y, r, :cyan)
   circlef(-x, y, r, :cyan)
   plot!([x0, 0, -x0, 0, x0], [-r, R - 2r, -r, -R, b - R], color=:red, seriestype=:shape, lw=0.5, alpha=0.5)
   circlef(x2, r - R, r, :cyan)
   segment(x2, -R, 0, -R, :cyan)
   plot!([x2, 0, 0, x2, x2], [-R, -R, 0, r - R, -R], color=:black, lw=1)
   circle(0, 0, 0.2, :black, beginangle=210, endangle=270)
   circle(x2, r - R, 0.2, :black, beginangle=270, endangle=390)
   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, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :vcenter)
       point(0, R - r, "小円:r,(0,R-r)", :green, :center, delta=-delta/2)
       point(x, r, "小円:r,(x,r)", :green, :center, delta=-delta/2)
       point(x2, r - R, "小円:r \n(x2,r-R) ", :green, :right, :vcenter)
       point(0, -r, " -r", :white, :left, :vcenter)
       point(2√2r, -r, "(2√2r,-r)  ", :white, :right, :vcenter)
       point(0, -0.25, "60°", :black, :right, mark=false)
       point(x2 + 0.2, r - R, "120°", :black, :left, mark=false)
   end
end;

算額（その859）

二十四　岩手県一関市萩荘字八幡 達古袋八幡神社 弘化3年

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

長方形内に大円が内接する不等辺三角形を入れ，大きさの同じ小円を左側に複数個，右側に 1 個入れる。小円の直径がわかっているとき，左側の小円の個数と大円の直径を求めよ。

なお，前もって述べるが，術では「(左側の小円の個数÷2 + 1)×小円径=大円径」とあるが，この式は一般には成り立たない。

また，左側の小円の個数が 1 個の場合は「算額（その857）」と同じになるが，その場合も，算額（その857）で述べた通り，示された条件だけでは不足であり，たとえば長方形の長辺の長さが既知でなければ一意な解はない。さらに，左側の小円の個数も未知として方程式を解くことは適切とは思えない。小円の個数も与えたうえで大円の直径を求めよとするのが妥当であろう。

長方形の長辺と短辺の長さを a, b
不等辺三角形の頂点の座標を (c, b)
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (r2, b - (2n - 1)*r2), (a - r2, b - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, n::integer,
     r1::positive, x1::positive, r2::positive, d
l1 = sqrt(b^2 + c^2)
l2 = sqrt(b^2 + (a - c)^2)
S1 = (a + l1 + l2)*r1/2
S2 = (b + c*(2n - 1) + l1)*r2/2
S3 = ((a - c) + b + l2)*r2/2
eq1 = (S2 + S3) - S1
eq2 = numerator(apart(dist(0, 0, c, b, x1, r1) - r1^2, d))
eq3 = numerator(apart(dist(c, b, a, 0, x1, r1) - r1^2, d))
eq4 = numerator(apart(dist(c, b, a, 0, a - r2, b - r2) - r2^2, d))
eq5 = (a - c) + b - l2 - 2r2;

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

左側の小円が 1 個の場合に，算額（その857）と同じになる。
算額（その857）は大円の直径から小円の直径を求めるものであるが，長方形の長辺が 2，大円の直径が 1 のとき，小円の直径が 2/3 である。
本算額では，以下のように小円の直径から大円の直径を求めることになる。

a = 2
n = 1
r2 = 1/3
iniv = BigFloat[2, 1.5, 0.6, 1.2]
res = nls(H, ini=iniv)

   ([1.3333333333333333, 1.0, 0.49999999999999994, 1.0], true)

長方形の長辺が 2，小円の直径が 2/3，個数が 1 のとき，大円の直径は 1 である。

「術」の「(左側の小円の個数÷2 + 1)×小円径=大円径」は一般には成り立たない。

a = 4;  b = 1.5;  c = 2;  n = 1;  r1 = 0.666667;  x1 = 2
小円の直径が 1，個数が 1，長方形の長辺が 4 のとき，大円の直径は 1.33333 である。
術による大円の直径 = 1.5

a = 4;  b = 2.70711;  c = 2.70711;  n = 2;  r1 = 1;  x1 = 2.41421
小円の直径が 1，個数が 2，長方形の長辺が 4 のとき，大円の直径は 2 である。
術による大円の直径 = 2

a = 4;  b = 4.10074;  c = 2.83875;  n = 3;  r1 = 1.23801;  x1 = 2.36272
小円の直径が 1，個数が 3，長方形の長辺が 4 のとき，大円の直径は 2.47602 である。
術による大円の直径 = 2.5

a = 4;  b = 5.53944;  c = 2.88985;  n = 4;  r1 = 1.39379;  x1 = 2.29917
小円の直径が 1，個数が 4，長方形の長辺が 4 のとき，大円の直径は 2.78757 である。
術による大円の直径 = 3

a = 4;  b = 7;  c = 2.91667;  n = 5;  r1 = 1.5;  x1 = 2.25
小円の直径が 1，個数が 5，長方形の長辺が 4 のとき，大円の直径は 3 である。
術による大円の直径 = 3.5

a = 20;  b = 15.5568;  c = 18.9657;  n = 15;  r1 = 5.17517;  x1 = 14.4693
小円の直径が 1，個数が 15，長方形の長辺が 20 のとき，大円の直径は 10.3503 である。
術による大円の直径 = 8.5

a = 20
n = 15
r2 = 1/2
iniv = BigFloat[2, 1.5, 0.6, 1.2].*10
res = nls(H, ini=iniv)
function draw(more=false)
   pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, c, r1, x1) = res[1]
   str = @sprintf("2r2=%g, n=%g, a=%g\n --> 2r1=%g", 2r2, n, a, 2r1)
   @printf("a = %g;  b = %g;  c = %g;  n = %g;  r1 = %g;  x1 = %g\n", a, b, c, n, r1, x1)
   @printf("小円の直径が %g，個数が %g，長方形の長辺が %g のとき，大円の直径は %g である。\n", 2r2, n, a, 2r1)
   @printf("術による大円の直径 = %g\n", (n/2 + 1)*2r2)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:blue, lw=0.5)
   circle(x1, r1, r1)
   circle(a - r2, b - r2, r2, :green)
   for i = 1:n
       circle(r2, b - (2i - 1)*r2, r2, :orange)
   end
   plot!([0, c, a], [0, b, 0], color=:magenta, lw=0.5)
   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(x1, 0.1r1, str, :black, :center, :bottom, delta=delta/2, mark=false)
       point(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(r2, b - (2n - 1)r2, "小円:r2,(r2,b-(2n-1)r2)", :black, :left, :bottom, delta=delta/2)
       point(a - r2, b - r2,"小円:r2,(a-r2,b-r2)", :black, :right, delta=-delta/2)
   end
end;

算額（その858）

二十二　岩手県一関市瑞山 駒形根神社 明治41年(1908)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正三角形内に甲円，乙円，丙円を各 1 個，丁円を 2 個入れる。甲円の直径がわかっているときに丁円の直径を求めよ。

注：引用元の図は誤解を招く不正確なものである。以下のように丙円，丁円は外接している。「術」に述べられている解になるためには，このような図でなければならない。

正三角形の一辺の長さを 2a
甲円の半径と中心座標を r1, (0, 2r3 + 2r2 + r1)
乙円の半径と中心座標を r2, (0, 2r3 + r2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, r2::positive, r3::positive,
     r4::positive, x4::positive, d
s3 = √Sym(3)
eq1 = 2r1 + r1 + 2r2 + 2r3 - s3*a
eq2 = (s3*a - (2r3 + r2)) - 2r2
eq3 = dist(a, 0, 0, s3*a, x4, r4) - r4^2 |> (x -> apart(x, d)) |> numerator
eq4 = x4^2 + ((2r3 + r2) - r4)^2 - (r2 + r4)^2
eq5 = x4^2 + (r4 - r3)^2 - (r3 + r4)^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r2, r3, r4, x4));

2 組の解が得られるが，2 番目のものが適解である。そのうちの 4 番目が丁円の半径。

res[2] |> println
res[2][4] |> println

   (r1*(-6*sqrt(265 + 156*sqrt(3)) + sqrt(795 + 468*sqrt(3)) + 132 + 110*sqrt(3))/33, 3*r1, r1*(-2*sqrt(795 + 468*sqrt(3)) + 11 + sqrt(265 + 156*sqrt(3)) + 44*sqrt(3))/22, r1*(-10*sqrt(265 + 156*sqrt(3)) - 2*sqrt(795 + 468*sqrt(3)) + 187 + 110*sqrt(3))/33, -3*sqrt(3)*r1*(-20*sqrt(265/324 + 13*sqrt(3)/27) + 70/9 + 20*sqrt(3)/3)/10)

   r1*(-10*sqrt(265 + 156*sqrt(3)) - 2*sqrt(795 + 468*sqrt(3)) + 187 + 110*sqrt(3))/33

甲円の半径に対する丁円の半径の倍数は，以下のように簡約化できる。

res1 = res[2][4]/r1
res1 |> println

   -10*sqrt(265 + 156*sqrt(3))/33 - 2*sqrt(795 + 468*sqrt(3))/33 + 17/3 + 10*sqrt(3)/3

res2 = (res1.args[2] + res1.args[3])^2 |> simplify |> sqrt |> simplify
res2 |> println

   2*sqrt(100 + 58*sqrt(3))/3

res3 = res1.args[1] + res1.args[4] - res2
res3 |> println

   -2*sqrt(100 + 58*sqrt(3))/3 + 17/3 + 10*sqrt(3)/3

甲円の半径に対する丁円の半径の倍数は約 2 である。正確に 2 ではないことに驚く。まさに奇跡。

res3 |> N

   2.001267660707633

術には以下の式が書かれている。見かけは簡単であるが，上述の結果と同じである。

天 = √3 + 1
(sqrt(8天 + 5) - 天)^2/3

   2.001267660707633

function draw(more=false)
   pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 0.5
   s3 = √3
   s = sqrt(795 + 468s3)
   t = sqrt(265 + 156s3)
   (a, r2, r3, r4, x4) = r1 .* (
       (s - 6t + 132 + 110s3)/33,
       3,
       (t - 2s + 11 +44s3)/22,
       (187 + 110s3 - 10t - 2s)/33,
       -3s3*(70/9 + 20s3/3 - 20sqrt(265/324 + 13s3/27))/10)
   @printf("甲円の直径が %g のとき，丁円の直径は %g\n", 2r1, 2r4)
   @printf("r1 = %g;  a = %g;  r2 = %g;  r3 = %g; r4 = %g;  x4 = %g\n", r1, a, r2, r3, r4, x4)
   plot([a, 0, -a, a], [0, s3*a, 0, 0], color=:blue, lw=0.5)
   circle(0, s3*a - 2r1, r1)
   circle(0, s3*a - 3r1 - r2, r2, :green)
   circle(0, r3, r3, :orange)
   circle2(x4, r4, r4, :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, s3*a - 2r1, " 甲円:r1,(0,√3*a-2r1)", :red, :left, :vcenter)
       point(0, s3*a - 3r1 - r2, "乙円:r2,(0,√3*a-3r1-r2)", :green, :center, delta=-delta/2)
       point(0, r3, "丙円;r3\n(0,r3)", :orange, :center, delta=-delta/2)
       point(x4, r4, "丁円;r4,(x4,r4)", :magenta, :center, delta=-delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
   end
end;

算額（その857）

二十二　岩手県一関市瑞山 駒形根神社 明治41年(1908)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

長方形の中に 2 本の弦と大円 1 個，小円 2 個を入れる。大円の直径が 1 寸のとき，小円の直径を求める術を述べよ。

長方形の長辺と短辺を 2a, b
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (a - r2, b - r2)
とおき，以下の連立方程式を解く。
なお，分析結果からいうと，「大円の直径が 1 寸」という条件だけでは不足である。r1 が一定でも，a の値により小円の直径は変わる。つまり，この図形を決めるためには，大円の直径と長方形の長辺の両方が必要である。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive, d::positive
eq1 = dist(a, 0, 0, b, 0, r1) - r1^2
eq1 = numerator(apart(eq1, d))
eq2 = dist(a, 0, 0, b, a - r2, b - r2) - r2^2
eq2 = numerator(apart(eq2, d))
res = solve([eq1, eq2], (b, r2))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (2*a^2*r1/((a - r1)*(a + r1)), a^2/(a - r1))
    (2*a^2*r1/((a - r1)*(a + r1)), a*r1/(a + r1))

2 組の解が得られるが，2 番目のものが適解である。
大円の直径が 1 寸で，長方形の長辺の長さが 2 寸のとき，小円の直径は 2/3 寸である。

r1 = 1/2  # 2r1 = 1
a = 1     # 2a = 2
2(a*r1/(a + r1))

   0.6666666666666666

図に示すように，r1 が一定でも a が大きくなるにつれて r2 も大きくなり，0.5 に収束する。つまり，r2 の最大値は r1 で，小円の直径は大円の直径に近づいていく。

f = res[2][2](r1 => 1//2)
f |> println

   a/(2*(a + 1/2))

pyplot(size=(500, 300), showaxis=true, grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(f, xlims=(0.6, 100), xlabel="a", ylabel="r2")

算額には「術」は書かれておらず，「答曰小円径六分〇九毛有奇」とのみある。山村は「図面は誤り。実証図はもっと大きい。」と書いているが，大円の直径に近づくとは気づいていない。

function draw(more=false)
   pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1//2
   a = 1
   (b, r2) = (2*a^2*r1/(a^2 - r1^2), a*r1/(a + r1))
   @printf("小円の直径 = %g\n", 2r2)
   @printf("r1 = %g;  a = %g;  b = %g;  r2 = %g\n", r1, a, b, r2)
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:blue)
   circle(0, r1, r1)
   circle2(a - r2, b - r2, r2, :green)
   plot!([-a, 0, a], [0, b, 0], color=:orange, lw=0.5)
   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(0, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta/2)
       point(a - r2, b - r2, "小円:r2\n(a-r2,b-r2)", :green, :center, delta=-delta/2)
       point(a, b - 2.1r2, @sprintf("r1=1/2,a=%g ", a), :black, :right, mark=false)
   end
end;

算額（その856）

二十二　岩手県一関市瑞山 駒形根神社 明治41年(1908)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

外円内に大円 2 個，小円 6 個が入っている。外円の直径が 10 寸のとき，小円の直径を得る術を問う。

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

include("julia-source.txt");

using SymPy

@syms R, r1, r2, x2, y2
r1 = R/2
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = (r1 - x2)^2 + y2^2 - (r1 + r2)^2
eq3 = x2^2 + (R - r2 - y2)^2 - 4r2^2;
res = solve([eq1, eq2, eq3], (r2, x2, y2));

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

(res[2][1](R=>5).evalf(),
res[2][2](R=>5).evalf(),
res[2][3](R=>5).evalf())

   (1.01331143350161, 1.96006569949517, 3.47157430856830)

それぞれは非常に長い式であるが，R が与えられれば小円の直径を求める「術」はある。以下のように計算すればよい。

# res[2][1] を求める関数を定義する

function 直径(R)
   f = ∛2
   h = f^2
   m = √69
   e = 69m + 997
   c = ∛e
   b = c^2
   a = 3h*b
   d = 330f + 37*c + a
   g = √d
   i = ∜d
   j = e^(5/6)
   k = d^(3/4)
   s = √e
   q = e^(1/6)
   r = (138m + 1994)^(2/3)
   n = sqrt(-a*g - 330f*g + 430*s + 74*c*g)
   o = sqrt(-207h*b*g - 22770f*g + 29670s + 5106c*g)
   p = sqrt(22770f + 2553c + 207h*b)
   R*(-17750573f*b*k*n/1242 - 154556461h*c*k*n/1242 - 2602054165k*n/1242 - 214763h*c*k*o/18 - 3923195k*o/18 - 21809f*b*k*o/18 - 4907453r*g/2 - 1625993r*p/6 - 341870903*c*g/9 - 11855327c*p/3 - 171811097f*g - 48431557f*p/3 + 658961h*m*s*i*n/9 + 28361659q*i*o/18 + 83669f*j*i*o/9 + 494602117h*s*i*n/621 + 18810834773q*i*n/1242 + 64451243f*j*i*n/621 + 8457376h*m*j + 488797196sqrt(4761*m + 68793)/3 + 1660341296f*m*q + 14095396844s/9 + 45213761744f*q/3 + 268114864h*j/3)/(19029096h*m*j + 366597897*m*s + 3735767916f*m*q + 3523849211s + 33910321308f*q + 201086148h*j)
end;
直径(10/2)

   1.01331143350161

算額では「術」はなく「答」も天下り式に「小円二寸」とある。
山村の記述では，「外円径÷5 = 小円径」としているが，こじつけである。

何度も経験することであるが，この文書に収められている算額および山村の記述には問題が多い。

function draw(more=false)
   pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 10/2
   r1 = R/2
   (r2, x2, y2) = (1.01331143350161, 1.96006569949517, 3.47157430856830)
   plot()
   circle(0, 0, R)
   circle2(r1, 0, r1, :blue)
   circle22(0, R - r2, r2, :magenta)
   circle4(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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(r1, 0, "大円:r1,(r1,0)", :blue, :center, delta=-delta/2)
       point(0, R - r2, "小円:r2\n(0,R-r2)", :magenta, :center, delta=-delta/2)
       point(x2, y2, "小円:r2\n(x2,y2)", :magenta, :center, delta=-delta/2)
   end
end;

算額（その855）

三十　岩手県一関市山ノ目 配志和神社 嘉永5年(1848)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正方形内に直径が 1 寸の 4 個の等円を入れる。赤で示した部分の面積を求めよ。

正方形の一辺の長さを 2a，等円の半径を r とおくと，正方形の一辺の長さは等円の直径の 1 + √2 倍である。

include("julia-source.txt");

using SymPy

@syms a, r
eq = (2r)^2 - 2(a - r)^2
a = solve(eq, a)[2]
a |> println

   r*(1 + sqrt(2))

赤積は，点 (0, 0), (a, 0), (a, a - r), (a - r, a - r) を結んでできる台形の面積から，緑で示した等円の面積の 45/360 と，赤で示した等円の面積の 135/360 を引いたものである。

赤積 = 2((a + r)*(a - r)/2 - pi*r^2*(1/8 + 135/360)) = (a^2 - r^2) - pi*r^2

r = 1/2 のとき，a = 1 + √2 なので，赤積 = 0.4217086177890992 である。

r = 1/2
a = r*(1 + sqrt(2))
(a^2 - r^2) - pi*r^2

   0.4217086177890992

また，緑の等円の（上側の）方程式 f(x) と 赤の等円の（下側）の方程式 g(x) について積分することで，赤積を求めることもできる。
当然ながら両者は一致する。

@syms r, a, x, f, g
r = 1//2
a = r*(1 + √Sym(2))

f = sqrt(r^2 - x^2)
S1 = integrate(f, (x, r/√(Sym(2)), r));
S1 |> println

g = -sqrt(r^2 - (x - (a - r))^2) + (a - r)
S2 = integrate(g, (x, r/√(Sym(2)), a));
S2 |> println

S = S2 - S1
S |> println

2S.evalf() |> println

   -1/16 + pi/32
   -5/16 - 3*pi/32 + sqrt(2)*(1/2 + sqrt(2)/2)/2
   -pi/8 - 1/4 + sqrt(2)*(1/2 + sqrt(2)/2)/2
   0.421708617789099

念の為，GeoGebra で図形を描き，面積を多角形で近似すると，ほぼ同じ数値になることを確かめた。

答曰も，術曰も赤積は 0.884 有奇としているが，明白な誤りである。
0.7854 は円積率で，0.7854 * 4 は円周率に極めて近い数値なのであるが，ほかは何を計算しているのだろうか。

(sqrt(0.75) + 1 - 0.7854 * 1.25) |> println
sqrt(0.75) + 1 |> println
0.7854 * 4 |> println
0.7854 * 1.25 |> println
1.866 - 0.98175 |> println

   0.8842754037844386
   1.8660254037844386
   3.1416
   0.98175
   0.8842500000000001

function draw(more=false)
   pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1//2
   a = r*(1 + sqrt(2))
   赤積 = 2((a + r)*(a - r)/2 - pi*r^2*(1/8 + 135/360))
   println("赤積 = $赤積")
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle4(a - r, a - r, r)
   circle(0, 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(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(a - r, a - r, "等円:r,(a-r,a-r)", :red, :center, delta=-delta/2)
       point(r/√2, r/√2, " (r/√2,r/√2)", :red, :left, :vcenter)
       point(r, 0, " r", :red, :left, delta=-delta/2)
   end
end;

算額（その854）

十八　岩手県平泉町 弁慶堂（現在は地蔵堂にて保管） 安政6年(1859)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正方形内に斜線を 2 本引き（右下頂点を通る甲斜，左下と右上の頂点を通る対角線の乙斜），正方形に内接する円（予円と書かれている）と乙円，丙円，丁円，戊円の 4 円を入れる。正方形の一辺の長さが与えられたとき，「甲斜，乙・丙・丁・戊・己円の直径の六和」を求めよ。

まず，「己円」が定義されていないという大問題がある。4 円の名前が乙円から始まっているのも不自然であるが，正方形に内接する円は「予円」と名付けられている。予と己を誤記するとも思えない。（山村の誤記か，もともとの誤記かもわからないが）

問は（術に示すように），「甲斜の引き方によらず，六和は正方形の一辺の長さに定数を掛けたものになる」ことを意図しているのであろうが，とりあえずは正方形の一辺の長さと甲斜の位置を与えて結果がどうなるかを見てみよう。

問を解く前に述べておくが，術は，「二個開平方以滅二個開平方加三個乗方面」とあり，山村は，
六和 = (sqrt(2 - √2) + 3)×方 = 4.8477×方
としているが，そもそも (sqrt(2 - √2) + 3) = 3.7653668647301792 である（これは山村の誤記であろう）。

以上の諸問題があるが，六和を求める以前に，図形を描くことができるかどうかを検討する。

正方形の一辺の長さを a, 甲斜と正方形の一辺との交点座標を (0, b)
乙円，丙円，丁円，戊円の半径と中心座標を，r2,(x2, y2) ... r5,(x5,y5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

dist2(x1, y1, x2, y2, x3, y3, r) = numerator(apart(dist(x1, y1, x2, y2, x3, y3) - r^2, d))
@syms d, a, b, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5
# (X0, Y0) = intersection(0, 0, a, a, a, 0, 0, b)
# (X1, Y1) = intersection(x2, y2, x5, y5, x3, y3, x4, y4)
eq1 = (x2 - a/2)^2 + (y2 - a/2)^2 - (a/2 - r2)^2 |> expand
eq2 = (x3 - a/2)^2 + (y3 - a/2)^2 - (a/2 - r3)^2 |> expand
eq3 = (x4 - a/2)^2 + (y4 - a/2)^2 - (a/2 - r4)^2 |> expand
eq4 = (x5 - a/2)^2 + (y5 - a/2)^2 - (a/2 - r5)^2 |> expand
eq5 = dist2(0, 0, a, a, x2, y2, r2)
eq6 = dist2(0, 0, a, a, x3, y3, r3)
eq7 = dist2(0, 0, a, a, x4, y4, r4)
eq8 = dist2(0, 0, a, a, x5, y5, r5)
eq9 = dist2(a, 0, 0, b, x2, y2, r2)
eq10 = dist2(a, 0, 0, b, x3, y3, r3)
eq11 = dist2(a, 0, 0, b, x4, y4, r4)
eq12 = dist2(a, 0, 0, b, x5, y5, r5)
# eq13 = X1 - X0 |> simplify |> numerator
# eq14 = Y1 - Y0 |> simplify |> numerator

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)
   (r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
a^2/4 + a*r2 - a*x2 - a*y2 - r2^2 + x2^2 + y2^2,  # eq1
a^2/4 + a*r3 - a*x3 - a*y3 - r3^2 + x3^2 + y3^2,  # eq2
a^2/4 + a*r4 - a*x4 - a*y4 - r4^2 + x4^2 + y4^2,  # eq3
a^2/4 + a*r5 - a*x5 - a*y5 - r5^2 + x5^2 + y5^2,  # eq4
-r2^2 + x2^2/2 - x2*y2 + y2^2/2,  # eq5
-r3^2 + x3^2/2 - x3*y3 + y3^2/2,  # eq6
-r4^2 + x4^2/2 - x4*y4 + y4^2/2,  # eq7
-r5^2 + x5^2/2 - x5*y5 + y5^2/2,  # eq8
a^2*b^2 - 2*a^2*b*y2 - a^2*r2^2 + a^2*y2^2 - 2*a*b^2*x2 + 2*a*b*x2*y2 - b^2*r2^2 + b^2*x2^2,  # eq9
a^2*b^2 - 2*a^2*b*y3 - a^2*r3^2 + a^2*y3^2 - 2*a*b^2*x3 + 2*a*b*x3*y3 - b^2*r3^2 + b^2*x3^2,  # eq10
a^2*b^2 - 2*a^2*b*y4 - a^2*r4^2 + a^2*y4^2 - 2*a*b^2*x4 + 2*a*b*x4*y4 - b^2*r4^2 + b^2*x4^2,  # eq11
a^2*b^2 - 2*a^2*b*y5 - a^2*r5^2 + a^2*y5^2 - 2*a*b^2*x5 + 2*a*b*x5*y5 - b^2*r5^2 + b^2*x5^2,  # eq12
   ]
end;

a = 100
b = 50
iniv = BigFloat[24.8, 82, 47,   26.7, 39, 77,   18.2, 45, 19,   15.7, 18, 40]
res = nls(H, ini=iniv)
println(res)

   ([24.042227875260373, 73.9205817737793, 39.91973704292161, 24.664395829520114, 28.463273173025932, 63.34399626287168, 13.777280329715465, 36.05369934249495, 16.569682647595275, 10.478025010198175, 15.644697884880152, 30.4628629611869], true)

「己円」が「予円」であるとしても，b を変えれば「六和」も変わる。
したがって，この問題は「無術」と思われる。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 100
   b = 50
   (r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = res[1]
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   circle(a/2, a/2, a/2, :brown)
   circle(x2, y2, r2)
   circle(x3, y3, r3, :green)
   circle(x4, y4, r4, :magenta)
   circle(x5, y5, r5, :orange)
   segment(0, 0, a, a, :green)
   segment(0, b, a, 0, :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(x2, y2, "乙円:r2,(x2,y2)", :red, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3,(x3,y3)", :magenta, :center, delta=-delta/2)
       point(x4, y4, "丁円:r4,(x4,y4)", :orange, :center, delta=-delta/2)
       point(x5, y5, "戊円:r5\n(x5,y5)", :green, :center, delta=-delta/2)
       point(a/2, a/2, " 予円", :brown, :left, :vcenter)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その853）

十七　岩手県平泉町 中尊寺金色堂（現存しない） 文政8年(1825)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

直角三角形（鈎股弦）内に隔斜を 2 本引き，区画された領域に等円を 2 個入れる。鈎，股，弦は，「5(鈎 + 股 + 弦) - 37(股 - 鈎) = 1」という条件式が成り立つとき（「鈎^2 + 股^2 = 弦^2」は当然），等円の直径はいかほどか。

鈎股弦内に隔斜 2 本と等円 2 個を入れるという問は他にもある。算額（その139），算額（その280）

この算額では，鈎股弦の条件の出し方が独特である。

もし大学入試の数学の問題だと整数方程式を如何に効率的に解くかということになるかもしれないが，プログラムで全探索するほうが手っ取り早い。
念の為， 鈎，股 が 10000 まで探索すると，解は 28, 45, 53 の一組しかない。

include("julia-source.txt");
#=
mx = 10000
for 鈎 = 1:mx, 股 = 鈎 + 1:mx
   弦 = sqrt(鈎^2 + 股^2)
   if 5(鈎 + 股 + 弦) - 37(股 - 鈎) == 1
       @printf("%g;  %g;  %g\n", 鈎, 股, 弦)
   end
end
=#

鈎 = 28, 股 = 45 円の半径を r，それぞれの中心座標を (r, y), (x, r) とおき，直線までの距離に関する方程式を立て，解く。
SymPy の solve() では能力不足で解けないので，nlsolve() で数値解を求めるのが簡単であるが，変数を手動で消去していって r を求める手順を示してみよう。

まずは，式を簡単にするために，円の中心から隔斜までの距離の式を通覧する。

using SymPy

dist2(x1, y1, x2, y2, x3, y3, r) = numerator(apart(dist(x1, y1, x2, y2, x3, y3) - r^2, d))
@syms 鈎::positive, 股::positive, 弦::positive,
     r::positive, x::positive, y::positive, x0::positive, y0::positive,
     d;
eq1 = dist2(0, 鈎, x0, y0, r, y, r)  # 円1の中心から BD への距離
eq2 = dist2(0, 鈎, x0, y0, x, r, r)  # 円2の中心から BD への距離
eq3 = dist2(0, 鈎, 股,  0, x, r, r)  # 円2の中心から AB への距離
eq4 = dist2(0, 0,  x0, y0, r, y, r)  # 円1の中心から OC への距離
eq5 = dist2(0, 0,  x0, y0, x, r, r); # 円2の中心から OC への距離

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   -r^2*x0^2 - 2*r*x0*y*y0 + 2*r*x0*y*鈎 + 2*r*x0*y0*鈎 - 2*r*x0*鈎^2 + x0^2*y^2 - 2*x0^2*y*鈎 + x0^2*鈎^2,  # eq1
   -r^2*y0^2 + 2*r^2*y0*鈎 - r^2*鈎^2 - 2*r*x*x0*y0 + 2*r*x*x0*鈎 - 2*r*x0^2*鈎 + x^2*y0^2 - 2*x^2*y0*鈎 + x^2*鈎^2 + 2*x*x0*y0*鈎 - 2*x*x0*鈎^2 + x0^2*鈎^2,  # eq2
   -r^2*鈎^2 + 2*r*x*股*鈎 - 2*r*股^2*鈎 + x^2*鈎^2 - 2*x*股*鈎^2 + 股^2*鈎^2,  # eq3
   -r^2*x0^2 - 2*r*x0*y*y0 + x0^2*y^2,  # eq4
   -r^2*y0^2 - 2*r*x*x0*y0 + x^2*y0^2,  # eq5

最初に，eq5 を x について解く。

res1 = solve(eq5, x0)[1]
res1 |> println

   y0*(-r^2 + x^2)/(2*r*x)

eq1 〜 eq4 に x0 = y0*(-r^2 + x^2)/(2*r*x) を代入する。結果の方程式は eq11 〜 eq14

eq11 = numerator(apart(eq1(x0 => res1), d))/y0/(x - r)/(x + r) |> factor;

eq12 = numerator(apart(eq2(x0 => res1), d))/鈎/(x - r)/(x + r) |> factor;

eq13 = numerator(apart(eq3(x0 => res1), d))/鈎 |> factor;

eq14 = numerator(apart(eq4(x0 => res1), d))/y0^2/(x - r)/(x + r) |> factor;

次に eq14 を x について解く。

res2 = solve(eq14, x)[1]
res2 |>  println

   r*(r + y)/(-r + y)

eq11 〜 eq14 に x = r*(r + y)/(-r + y) を代入する。

eq21 = numerator(apart(eq11(x => res2), d))/(4r^3*鈎) |> factor;

eq22 = numerator(apart(eq12(x => res2), d))/(4r^3) |> factor;

eq23 = numerator(apart(eq13(x => res2), d)) |> factor;

次に eq21 を y0 について解く。

res3 = solve(eq21, y0)[1]
res3 |> println

   (-r^2*y + r^2*鈎 + y^3 - y^2*鈎)/(r^2 + y^2 - y*鈎)

eq21 〜 eq23 に y0 = (-r^2*y + r^2*鈎 + y^3 - y^2*鈎)/(r^2 + y^2 - y*鈎) を代入する。

eq32 = numerator(apart(eq22(y0 => res3), d))/(y*(r + y)^2) |> factor
eq32 |> println

   (-r^2 - 2*r*y + 2*r*鈎 + y^2 - y*鈎)*(-r^3 + r^2*y - r^2*鈎 - r*y^2 + r*鈎^2 + y^3 - y^2*鈎)

eq33 = numerator(apart(eq23(y0 => res3), d)) |> factor;

eq32 は 式1*式2 の形で，-r^2 - 2*r*y + 2*r*鈎 + y^2 - y*鈎 = 0 を解けばよい。

eq32_1 = (-r^2 - 2*r*y + 2*r*鈎 + y^2 - y*鈎);

eq32_1 を y について解く。

res4 = solve(eq32_1, y)[1]
res4 |> println

   r + 鈎/2 - sqrt(8*r^2 - 4*r*鈎 + 鈎^2)/2

eq33 に y = r + 鈎/2 - sqrt(8*r^2 - 4*r*鈎 + 鈎^2)/2 を代入する。

eq42 = eq33(y => res4) |> simplify
eq42 |> factor |> println

   (2*r*股 + 2*r*鈎 - 股*鈎)*(4*r^3 - 4*r^2*股 + 2*r^2*鈎 - 2*r^2*sqrt(8*r^2 - 4*r*鈎 + 鈎^2) + 2*r*股*鈎 - 股*鈎^2 + 股*鈎*sqrt(8*r^2 - 4*r*鈎 + 鈎^2))/2

sqrt(8*r^2 - 4*r*鈎 + 鈎^2) でまとめ，二乗してルートを外す。

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r::positive, x::positive, y::positive, x0::positive, y0::positive

eq01 = (4*r^3 - 4*r^2*股 + 2*r^2*鈎  + 2*r*股*鈎 - 股*鈎^2)^2 - (股*鈎 - 2*r^2)^2*(8*r^2 - 4*r*鈎 + 鈎^2) |> factor
eq01 |> println

   -4*r^2*(4*r^4 + 8*r^3*股 - 8*r^3*鈎 - 4*r^2*股^2 - 8*r^2*股*鈎 + 4*r*股^2*鈎 + 4*r*股*鈎^2 - 股^2*鈎^2)

eq02 = eq01/(-4r^2)
eq02 |> println

   4*r^4 + 8*r^3*股 - 8*r^3*鈎 - 4*r^2*股^2 - 8*r^2*股*鈎 + 4*r*股^2*鈎 + 4*r*股*鈎^2 - 股^2*鈎^2

鈎，股に具体値を代入してプロットしてみる。

eq03 = eq02(鈎 => Sym(28), 股 => Sym(45))

   4*r^4 + 8*r^3*股 - 8*r^3*鈎 - 4*r^2*股^2 - 8*r^2*股*鈎 + 4*r*股^2*鈎 + 4*r*股*鈎^2 - 股^2*鈎^2

using Plots
pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho");
plot(eq03, xlims=(4, 32))
hline!([0])

3 個の実数解が得られるが r = 6 が適解である。

solve(eq03, r)

   3-element Vector{Sym{PyCall.PyObject}}:
                   6
                  30
    -35 + 7⋅√70

あとは，交代代入を繰り返し，ここまでで変数消去されていた変数の値を確定していけばよい。

(鈎, 股) = (28, 45)
r = 6
y = r + 鈎/2 - sqrt(8*r^2 - 4*r*鈎 + 鈎^2)/2
y0 = (-r^2*y + r^2*鈎 + y^3 - y^2*鈎)/(r^2 + y^2 - y*鈎)
x = r*(r + y)/(-r + y)
x0 = y0*(-r^2 + x^2)/(2*r*x)
(r, y, y0, x, x0)

   (6, 10.0, 8.0, 24.0, 15.0)

NLsolve() で数値解を求めるほうが，簡単である。

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number</span>
       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, x, y, x0, y0) = u
   return [
       -r^2*x0^2 - 2*r*x0*y*y0 + 2*r*x0*y*鈎 + 2*r*x0*y0*鈎 - 2*r*x0*鈎^2 + x0^2*y^2 - 2*x0^2*y*鈎 + x0^2*鈎^2,  # eq1
       -r^2*y0^2 + 2*r^2*y0*鈎 - r^2*鈎^2 - 2*r*x*x0*y0 + 2*r*x*x0*鈎 - 2*r*x0^2*鈎 + x^2*y0^2 - 2*x^2*y0*鈎 + x^2*鈎^2 + 2*x*x0*y0*鈎 - 2*x*x0*鈎^2 + x0^2*鈎^2,  # eq2
       -r^2*鈎^2 + 2*r*x*股*鈎 - 2*r*股^2*鈎 + x^2*鈎^2 - 2*x*股*鈎^2 + 股^2*鈎^2,  # eq3
       -r^2*x0^2 - 2*r*x0*y*y0 + x0^2*y^2,  # eq4
       -r^2*y0^2 - 2*r*x*x0*y0 + x^2*y0^2,  # eq5
   ]
end;

(鈎, 股, 弦) = (28, 45, 53)
iniv = BigFloat[5.0, 20, 11, 12, 10]
res = nls(H, ini=iniv)

   ([6.0, 24.0, 10.0, 15.0, 8.0], true)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦) = (28, 45, 53)
   (r, x, y, x0, y0) = res[1]
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(r, y, r)
   circle(x, r, r)
   (X1, Y1) = intersection(0, 0, x0, y0, 0, 鈎, 股, 0)
   segment(0, 0, X1, Y1, :green)
   (X2, Y2) = intersection(0, 鈎, x0, y0, 0, 0, 股, 0)
   segment(0, 鈎, X2, Y2, :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(x, r, "等円:r,(x,r)", :red, :center, delta=-delta/2)
       point(r, y, "等円:r,(r,y)", :red, :center, delta=-delta/2)
       point(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=4delta, deltax=-4delta)
       point(股, 0, " 股", :blue, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その852）

十六　岩手県前沢町赤生津箱根 金刀比羅神社 明治43年(1910)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

外円の中に正方形，さらにその中に 7 個の円が入っている。小円の直径が 50 寸のとき，外円の直径はいくつか。

算額（その77）を 90 度回転し，外円を加えたものに過ぎない。求めるものも，外円の直径になっている。

外円の半径と中心座標を R, (0, 0)
正方形の一辺の長さを 2a = √2R
大円の半径と中心座標を r0, (r1, r0), (-r1, r0); r0 = 2r1
中円の半径と中心座標を r1, (0, 0), (r0, r1)
小円の半径と中心座標を r2, (0, a - r2)
とおき，以下の方程式を解く。

include("julia-source.txt");

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

   5*sqrt(2)*r2

外円の直径は 小円の直径の 5√2 倍である。
小円の直径が 50 寸のとき，外円の直径は 353.5533905932738 になる。
√2 ≒ 1.414 で近似すると 353.5 である。

しかし，（判読不能箇所が多いが）問の末尾近くの一部分には「▢▢只▢五十二有寸外円径問幾何」とあり，52というのは小円の径を言っているのだろうと推測される。小円の径が 52 なら，外円の径は 353.5 にはならない。
術は「置二個開平方乗小円径五段得外円径」と正しいことをいっているのだが。

なお，正方形の長さは小円の直径の 3 倍，中円の直径は小円の直径の 5/3 倍である。

50*5*1.414, 50*5√2

   (353.5, 353.5533905932738)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 50//2
   R = 5√2r2
   a = R/√2
   r1 = a/3
   @printf("外円の直径 = %g;  正方形の一辺の長さ = %g;  中円の直径 = %g;  小円の直径 = %g\n", 2R, 2a, 2r1, 2r2)
   plot(a .* [-1, 1, 1, -1, -1], a .* [-1, -1, 1, 1, -1], color=:black, lw=0.5)
   circle(0, 0, R)
   circle(0, 0, r1, :green)
   circle2(a - r1, 0, r1, :green)
   circle2(r1, 0, 2r1, :blue)
   circle22(0, a - r2, 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, 0, " O", :green, :center, delta=-delta/2)
       point(r1, 0, " 大円:r1\n (r1,0)", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :green, :center, delta=-delta/2)
       point(2r1, 0, "中円:r1\n(2r1,0)", :green, :left, deltax=-2delta, delta=-delta/2)
       point(0, a - r2, " 小円:r2,(0,a-r2)", :black, :left, :vcenter)
       point(a, 0, " a ", :magenta, :left, :vcenter)
       point(R, 0, " R ", :red, :left, :vcenter)
   end
end;

算額（その851）

九　岩手県水沢市佐倉河 胆沢城八幡神社 弘化2年(1845)

山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

大円内に互いに内接・外接する小円 3 個と，菱形を入れる。小円の直径が与えられたとき，黒積を求めよ。
図では，小円，菱形はどれにも接していないが，そんな訳はないだろう。

常識なのかもしれないが，どれが黒積か明瞭ではない。大円の面積から小円および菱形の面積を除いた部分（図では緑色の部分）の面積かもしれないが，術は理解できない。

とりあえず，必要なパラメータを求め，図だけを描く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive, y::positive, d
y = r
eq1 = x^2 + y^2 - (R - r)^2
eq2 = x^2 + (R - r - y)^2 - 4r^2
res = solve([eq1, eq2], (R, x))[1]

   (3*r, sqrt(3)*r)

小円の直径が 1 のとき，大円の直径は 3，菱形の対角線は 4√2 と 2 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (R, x) = r .* (3, √3)
   y = r
   b = R - r
   y0 = -r
   x0 = sqrt(R^2 - y0^2)
   plot([x0, 0, -x0, 0, x0], [-r, R - 2r, -r, -R, b - R], color=:blue, lw=0.5)
   circle(0, 0, R)
   circle(0, R - r, r, :green)
   circle2(x, y, r, :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(0, R, " R", :red, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :vcenter)
       point(0, R - r, "小円:r,(0,R-r)", :green, :center, delta=-delta/2)
       point(x, r, "小円:r,(x,r)", :green, :center, delta=-delta/2)
       point(0, -r, " -r", :blue, :left, :vcenter)
       point(2√2r, -r, "(2√2r,-r)  ", :blue, :right, :vcenter)
   end
end;