算額（その466）

算額（その466）

愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)

平田浩一: 伊佐爾波神社の算額にみる江戸末期の和算，愛媛大学教育学部紀要，第60巻，195-206, 2013.
https://www.ed.ehime-u.ac.jp/~kiyou/2013/pdf/20.pdf

円弧内に青，黃，赤，白，黒の 5 個の円がある。青，赤，黒の 3 個の円の直径が与えられたとき，白円の直径を求めよ。

円弧の半径と中心座標を r0, (x0, y0)
青円の半径と中心座標を r1, (x1, y1)
黃円の半径と中心座標を r2, (x2, y2)
赤円の半径と中心座標を r3, (x3, y3); y3 = r3
白円の半径と中心座標を r4, (x4, y4); y4 = r4
黒円の半径と中心座標を r5, (x5, y5); x5 = 0, y5 = r5
とおき，以下の連立方程式を解き，数値解を求める。

include("julia-source.txt")

using SymPy

@syms r1::positive, x1::positive, y1::positive, 
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive,
     r4::positive, x4::positive,
     r5::positive,
     r0::positive, x0::positive, y0::positive;

x5 = 0
y5 = r5
y3 = r3
y4 = r4
eq1 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq3 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq4 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq6 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq7 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = (x0 - x1)^2 + (y0 - y1)^2 - (r0 - r1)^2
eq10 = (x0 - x3)^2 + (y0 - y3)^2 - (r0 - r3)^2
eq11 = (x0 - x5)^2 + (y0 - y5)^2 - (r0 - r5)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11], (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0))

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0) = u
   return [
-(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2,  # eq1
-(r1 + r3)^2 + (-r3 + y1)^2 + (x1 - x3)^2,  # eq2
x1^2 - (r1 + r5)^2 + (-r5 + y1)^2,  # eq3
-(r2 + r3)^2 + (-r3 + y2)^2 + (x2 - x3)^2,  # eq4
(r3 - r4)^2 - (r3 + r4)^2 + (x3 - x4)^2,  # eq5
-(r2 + r4)^2 + (-r4 + y2)^2 + (x2 - x4)^2,  # eq6
x2^2 - (r2 + r5)^2 + (-r5 + y2)^2,  # eq7
x4^2 + (r4 - r5)^2 - (r4 + r5)^2,  # eq8
-(r0 - r1)^2 + (x0 - x1)^2 + (y0 - y1)^2,  # eq9
-(r0 - r3)^2 + (-r3 + y0)^2 + (x0 - x3)^2,  # eq10
x0^2 - (r0 - r5)^2 + (-r5 + y0)^2,  # eq11
   ]
end;

青，赤，黒の 3 個の円の直径がそれぞれ 15, 30, 18 のときの解を求める。

(r1, r3, r5) = (15, 30, 18)
iniv = (r1/31) .* [big"15.0", 45, 6, 21, 25, 62, 10, 25, 57, 37, -15]
iniv = r1 .* [big"7.258", 21.774, 2.903, 10.161, 12.097, 30.0, 4.839, 12.097, 27.581, 17.903, -7.258]
res = nls(H, ini=iniv)
names = ("x1", "y1", "r2", "x2", "y2", "x3", "r4", "x4", "r0", "x0", "y0")
if res[2]
   for (name, x) in zip(names, res[1])
       @printf("%s = %g;  ", name, round(Float64(x), digits=6))
   end
   println()
else
   println("収束していない")
end;
@printf("白円の直径は %g である。\n", 2res[1][7])

   x1 = 17.5989;  y1 = 45.9156;  r2 = 6.54237;  x2 = 23.4663;  y2 = 25.1876;  x3 = 59.6904;  r4 = 9.42818;  x4 = 26.0544;  r0 = 59.7271;  x0 = 40.362;  y0 = 7.41423;  

白円の直径は 18.8564 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0) = res[1]
   @printf("x1 = %g;  y1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  r4 = %g;  x4 = %g;  r0 = %g;  x0 = %g;  y0 = %g\n", x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0)
   @printf("白円の直径は %g である。\n", 2r4)
   plot()
   circle(x1, y1, r1)
   circle(x2, y2, r2)
   circle(x3, r3, r3)
   circle(x4, r4, r4)
   circle(0, r5, r5)
   a = sqrt(r0^2 - y0^2)
   segment(x0 - a, 0, x0 + a, 0, :magenta, lw=0.5)
   θ = atand(y0/a)
   circle(x0, y0, r0, :magenta, beginangle=-θ, endangle=181 + θ)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        point(x1, y1, " 青:r1\n (x1,y1)", :blue, :center, :bottom, delta=delta)
        point(x2, y2, " 黃:r2,(x2,y2)", :black, :left, :vcenter)
        point(x3, r3, " 赤:r3,(x3,r3)", :red, :center, :bottom, delta=delta)
        point(x4, r4, "白:r4,(x4,r4)", :black, :center, :top, delta=-delta)
        point(x5, r5, "黒:r5,(x5,r5)", :black, :center, :top, delta=-delta)
        point(x0, y0, " 円弧:r0,(x0,y0)", :magenta, :left, :bottom, delta=delta)
     hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その465）

算額（その465）

愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)

平田浩一: 伊佐爾波神社の算額にみる江戸末期の和算，愛媛大学教育学部紀要，第60巻，195-206, 2013.
https://www.ed.ehime-u.ac.jp/~kiyou/2013/pdf/20.pdf

正三角形の中に2本の斜線で区切られた領域それぞれに直径の等しい円が内接している。下側の斜線と底辺を左に伸ばし，その2直線と正三角形の左側の斜辺に内接する外接円の直径を求めよ。

正三角形の左隅を原点におき，一辺の長さを a とする。
3個の等円の半径と中心座標を r, (x1, r), (x2, y2), (x3, y3)
外円の半径と中心座標を r0, (x0, r0)
斜線と正三角形の斜辺の交点座標を (x4, √3x4), (x5, (a - x5)√3)
とおき以下の連立方程式の数値解を求める。

include("julia-source.txt")

using SymPy

@syms a::positive, r::positive, x1::positive, x2::negative, y2::positive,
     x3::positive, y3::negative, x4::positive, x5::positive,
     r0::positive, x0::negative;

eq1 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x1, r) - r^2
eq2 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x3, y3) - r^2
eq3 = distance(a, 0, a/2, sqrt(Sym(3))a/2, x2, y2) - r^2
eq4 = distance(a, 0, a/2, sqrt(Sym(3))a/2, x3, y3) - r^2
eq5 = distance(x4, sqrt(Sym(3))x4, x5, (a - x5)sqrt(Sym(3)), x2, y2) - r^2
eq6 = distance(x4, sqrt(Sym(3))x4, x5, (a - x5)sqrt(Sym(3)), x3, y3) - r^2
eq7 = distance(x4, sqrt(Sym(3))x4, a, 0, x1, r) - r^2
eq8 = distance(x4, sqrt(Sym(3))x4, a, 0, x2, y2) - r^2
eq9 = distance(x4, sqrt(Sym(3))x4, a, 0, x0, r0) - r0^2  # 外円
eq10 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x0, r0) - r0^2;  # 外円

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r0, x0, r, x1, x2, y2, x3, y3, x4, x5) = u
   return [
-r^2 + (r/4 - sqrt(3)*x1/4)^2 + (-sqrt(3)*r/4 + 3*x1/4)^2,  # eq1
-r^2 + (3*x3/4 - sqrt(3)*y3/4)^2 + (-sqrt(3)*x3/4 + y3/4)^2,  # eq2
-r^2 + (-3*a/4 + 3*x2/4 + sqrt(3)*y2/4)^2 + (-sqrt(3)*a/4 + sqrt(3)*x2/4 + y2/4)^2,  # eq3
-r^2 + (-3*a/4 + 3*x3/4 + sqrt(3)*y3/4)^2 + (-sqrt(3)*a/4 + sqrt(3)*x3/4 + y3/4)^2,  # eq4
-r^2 + (x2 - (3*a^2*x4 - 3*a*x4^2 - 9*a*x4*x5 - sqrt(3)*a*x4*y2 + sqrt(3)*a*x5*y2 + x2*x4^2 - 2*x2*x4*x5 + x2*x5^2 + 6*x4^2*x5 + sqrt(3)*x4^2*y2 + 6*x4*x5^2 - sqrt(3)*x5^2*y2)/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2 + (y2 - (y2*(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2) + (x4 - x5)*(-sqrt(3)*a*x2 + sqrt(3)*a*x4 + sqrt(3)*x2*x4 + sqrt(3)*x2*x5 - 2*sqrt(3)*x4*x5 - x4*y2 + x5*y2))/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2,  # eq5
-r^2 + (x3 - (3*a^2*x4 - 3*a*x4^2 - 9*a*x4*x5 - sqrt(3)*a*x4*y3 + sqrt(3)*a*x5*y3 + x3*x4^2 - 2*x3*x4*x5 + x3*x5^2 + 6*x4^2*x5 + sqrt(3)*x4^2*y3 + 6*x4*x5^2 - sqrt(3)*x5^2*y3)/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2 + (y3 - (y3*(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2) + (x4 - x5)*(-sqrt(3)*a*x3 + sqrt(3)*a*x4 + sqrt(3)*x3*x4 + sqrt(3)*x3*x5 - 2*sqrt(3)*x4*x5 - x4*y3 + x5*y3))/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2,  # eq6
-r^2 + (r - sqrt(3)*x4*(a^2 - a*x1 - a*x4 + sqrt(3)*r*x4 + x1*x4)/(a^2 - 2*a*x4 + 4*x4^2))^2 + (x1 - (x4*(a^2 - 2*a*x4 + 4*x4^2) - (a - x4)*(-a*x1 + a*x4 + sqrt(3)*r*x4 + x1*x4 - 4*x4^2))/(a^2 - 2*a*x4 + 4*x4^2))^2,  # eq7
-r^2 + (x2 - (x4*(a^2 - 2*a*x4 + 4*x4^2) - (a - x4)*(-a*x2 + a*x4 + x2*x4 - 4*x4^2 + sqrt(3)*x4*y2))/(a^2 - 2*a*x4 + 4*x4^2))^2 + (-sqrt(3)*x4*(a^2 - a*x2 - a*x4 + x2*x4 + sqrt(3)*x4*y2)/(a^2 - 2*a*x4 + 4*x4^2) + y2)^2,  # eq8
-r0^2 + (r0 - sqrt(3)*x4*(a^2 - a*x0 - a*x4 + sqrt(3)*r0*x4 + x0*x4)/(a^2 - 2*a*x4 + 4*x4^2))^2 + (x0 - (x4*(a^2 - 2*a*x4 + 4*x4^2) - (a - x4)*(-a*x0 + a*x4 + sqrt(3)*r0*x4 + x0*x4 - 4*x4^2))/(a^2 - 2*a*x4 + 4*x4^2))^2,  # eq9
-r0^2 + (r0/4 - sqrt(3)*x0/4)^2 + (-sqrt(3)*r0/4 + 3*x0/4)^2,  # eq10
   ]
end;

a = 20
iniv = [big"16.0", -10, 10, 19, 47, 21, 35, 42, 14, 52] .* (a/71)
iniv = a .* [big"0.225", -0.141, 0.141, 0.268, 0.662, 0.296, 0.493, 0.592, 0.197, 0.732]
res = nls(H, ini=iniv);
names = ("r0", "x0", "r", "x1", "x2", "y2", "x3", "y3", "x4", "x5")
if res[2]
   for (name, x) in zip(names, res[1])
       @printf("%s = %g;  ", name, round(Float64(x), digits=6))
   end
   println()
else
   println("収束していない")
end;

   r0 = 4.50197;  x0 = -2.59921;  r = 2.96213;  x1 = 5.13055;  x2 = 13.0676;  y2 = 6.08309;  x3 = 10;  y3 = 11.3963;  x4 = 3.86488;  x5 = 14.4977;  

三角形の一辺の長さが 20 寸のとき，外円の直径は 9.00393 寸である。ちなみに，等円の直径は 5.92425 寸である。

using Plots

function abline(x0, y0, slope, minx, maxx, color=:black; lw=0.5)
   f(x) = slope * (x - x0) + y0
   segment(minx, f(minx), maxx, f(maxx), color; lw=lw)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 20
   (r0, x0, r, x1, x2, y2, x3, y3, x4, x5) = res[1]
   @printf("r0 = %g;  x0 = %g;  r = %g;  x1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g;  x4 = %g;  x5 = %g\n", r0, x0, r, x1, x2, y2, x3, y3, x4, x5)
   @printf("外円の直径 = %g;  等円の直径 = %g\n", 2r0, 2r)
   plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:green, lw=0.5)
   circle(x1, r, r, :blue)
   circle(x2, y2, r, :blue)
   circle(x3, y3, r, :blue)
   circle(x0, r0, r0)
   segment(x4, √3x4, x5, (a - x5)√3)
   abline(x4, √3x4, √3x4/(x4 - a), -0.4a, a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r, "(x1,r)", :blue, :center, :top, delta=-delta)
       point(x2, y2, "(x2,y2)", :blue, :center, :top, delta=-delta)
       point(x3, y3, "(x3,y3)", :blue, :center, :top, delta=-delta)
       point(x0, r0, "(x0,r0)", :red, :center, :top, delta=-delta)
       point(x4, √3x4, "  (x4,√3x4)", :black, :left, :vcenter)
       point(x5, √3(a-x5), " (x5,√3(a-x5))", :black, :left, :vcenter)
       point(a/2, √3a/2, " (a/2,√3a/2)", :black, :left, :vcenter)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その464）

算額（その464）

算額（その463）「埼玉県秩父市大宮 秩父神社 明治20年」を解いていて，条件を間違えてしまってできた図が案外きれいだったので記録しておこう。

山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

大円内に甲円 2 個，乙円 4 個，丙円　4 個，丁円 2 個が入っている。甲円の直径が 1 寸のとき，乙円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (0, r4)
として以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::negative,
     r3::positive, x3::negative, r4::positive;
#@syms r0, r1, r2, x2, r3, x3, r4;

eq1 = x2^2 + r2^2 - (r0 - r2)^2
eq2 = x2^2 + (r0 - r1 - r2)^2 - (r1 + r2)^2
eq3 = x3^2 + (r0 - r1 - r3)^2 - (r1 + r3)^2
eq4 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq5 = x3^2 + (r3 - r4)^2 - (r3 + r4)^2
eq6 = 2r4 + 2r1 - r0
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r0, r2, x2, r3, x3, r4))

   1-element Vector{NTuple{6, Sym}}:
    (4*r1*(sqrt(2) + 3)/7, r1*(5 + 4*sqrt(2))/14, -2*r1*sqrt(7*sqrt(2) + 21)/7, 2*r1*(1 + 5*sqrt(2))/49, 2*r1*sqrt(sqrt(2) + 3)*(-4*sqrt(7) + sqrt(14))/49, r1*(-1 + 2*sqrt(2))/7)

乙円の直径は甲円の直径の (5 + 4*sqrt(2))/14 倍である。

すなわち甲円の直径が 1 寸のとき，乙円の直径は 0.7612038749637414 = 7分6厘1毛あまりである。

r0 = 1.2612;  r1 = 0.5;  r2 = 0.380602;  x2 = -0.794104;  r3 = 0.164716;  x3 = -0.293341;  r4 = 0.130602
甲円の直径 = 1;  乙円の直径 = 0.761204

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   (r0, r2, x2, r3, x3, r4) = (4*r1*(sqrt(2) + 3)/7, r1*(5 + 4*sqrt(2))/14, -2*r1*sqrt(7*sqrt(2) + 21)/7, 2*r1*(1 + 5*sqrt(2))/49, 2*r1*sqrt(sqrt(2) + 3)*(-4*sqrt(7) + sqrt(14))/49, r1*(-1 + 2*sqrt(2))/7)
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g\n", r0, r1, r2, x2, r3, x3, r4)
   @printf("甲円の直径 = %g;  乙円の直径 = %g\n", 2r1, 2r2)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r1 - r0, r1, :blue)
   circle4(x2, r2, r2, :orange)
   circle4(x3, r3, r3, :green)
   circle(0, r4, r4)
   circle(0, -r4, r4)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0 - r1, " 甲円:r1\n (0,r0-r1)", :blue, :left, :vcenter)
       point(x2, r2, " 乙円:r2\n (x2,r2)", :orange, :left, :vcenter)
       point(x3, r3, "丙円:r3(r3+r4,0) ", :black, :right, :vcenter)
       point(0, r4, " 丁円:r4(0,r4)", :black, :left, :vcenter)
       point(r0, 0, "r0 ", :red, :right, :bottom, delta=delta/4)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

 

 

算額（その463）

算額（その463）

埼玉県秩父市大宮 秩父神社 明治20年(1887)
山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

九三 秩父市大宮 秩父神社 明治20年(1887)
一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

大円内に甲円 2 個，乙円 4 個，丙円　2 個，丁円 1 個が入っている。甲円の直径が 1 寸のとき，乙円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (r3 + r4, 0)
乙円の半径と中心座標を r4, (0, 0)
として以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive,
     r3::positive, r4::positive;

eq1 = x2^2 + r2^2 - (r0 - r2)^2
eq2 = x2^2 + (r0 - r1 - r2)^2 - (r1 + r2)^2
eq3 = (r3 + r4)^2 + (r0 - r1)^2 - (r1 + r3)^2
eq4 = (x2 - (r3 + r4))^2 + r2^2 - (r2 + r3)^2
eq5 = r4 + 2r1 - r0
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r2, x2, r3, r4))

   1-element Vector{NTuple{5, Sym}}:
    (9*r1/4, 5*r1/8, 3*r1/2, 5*r1/12, r1/4)

乙円の直径は甲円の直径の 5/8 である。
すなわち甲円の直径が 1 寸のとき，乙円の直径は 5/8 = 0.625 = 6分2厘5毛である。

r0 = 1.125;  r1 = 0.5;  r2 = 0.3125;  x2 = 0.75;  r3 = 0.208333;  r4 = 0.125
甲円の直径 = 1;  乙円の直径 = 0.625

using Plots

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

算額（その462）

算額（その462）

埼玉県秩父市大宮 秩父神社 明治20年

山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

正三角形内に半円 3 個，大円 1 個，小円 3 個が入っている。小円の直径が 1 寸のとき，大円の直径はいかほどか。

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

なお，eq1 を解いて r0 を求めれば，順次 eq2, eq3 を解いてもよい。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, r2::positive, a::positive
eq1 = distance(0, sqrt(Sym(3))a, -a, 0, a/2, sqrt(Sym(3))a/2) - r0^2
eq2 = (a/2)^2 + (sqrt(Sym(3))a/2 - a/sqrt(Sym(3)))^2 - (r0 - r1)^2
eq3 = (a/2)^2 + (sqrt(Sym(3))a/2 - r2)^2 - (r0 + r2)^2;

res = solve([eq1, eq2, eq3], (r0, r1, r2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(3)*a/2, sqrt(3)*a/6, sqrt(3)*a/24)
    (sqrt(3)*a/2, 5*sqrt(3)*a/6, sqrt(3)*a/24)

解は二通り得られるが，最初のものが適解である。
また，大円の半径は小円半径の4倍なので，小円の直径が 1 寸のとき，大円の直径は 4 寸である。

res[1][2]/res[1][3] |> println

   4

ちなみに，半円の半径は小円の半径の 12 倍である。

res[1][1]/res[1][3] |> println

   12

また，a = 4√3 のとき 小円の直径が 1 寸になる。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 4√3
   (r0, r1, r2) = (sqrt(3)*a/2, sqrt(3)*a/6, sqrt(3)*a/24)
   @printf("r0 = %g;  r1 = %g;  r2 = %g\n", r0, r1, r2)
   @printf("r1/r2 = %g\n", r1/r2)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   circle(0, 0, r0, beginangle=0, endangle=180)
   circle(a/2, √3a/2, r0, beginangle=120, endangle=300)
   circle(-a/2, √3a/2, r0, beginangle=240, endangle=420)
   circle(0, a/√3, r1, :blue)
   circle(0, r2, r2, :green)
   circle(0, r2, r2, :green)
   circle(a/2 - √3r2/2, √3a/2 - r2/2, r2, :green)
   circle(-a/2 + √3r2/2, √3a/2 - r2/2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a/2, √3a/2, " 半円:r0,(a/2,√3a/2)", :red, :left, :vcenter)
       point(0, r2, "  小円:r2,(0,r2)", :green, :left, :bottom, delta=2delta)
       point(a/2 - r2√3/2, √3a/2 - r2/2, "     小円:r2\n      (a/2-r2√3/2,√3a/2-r2/2)", :green, :left, :top, delta=-delta)
       point(0, a/√3, " 大円:r1\n (0,a/√3)", :black, :left, :bottom, delta=delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(0, √3a, " √3a", :black, :left, :vcenter)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

算額（その461）

算額（その461）

埼玉県秩父市大宮 秩父神社 明治20年(1887)

山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

外円内に大円，中円，小円と弦が図のように配置されている。中円の直径が 3 寸のとき，小円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (0, r1), (0, -r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (x31, 2r3), (x32, 0)
弦が y 軸と交わる y 座標は r3
として以下の連立方程式を解く。

なお，r2 が既知として方程式を解くと SymPy ではなぜか解が求まらない。
r2 を未知して，他の未知数も r0 により表せるとすれば解は求まる。
中円の直径が 3 寸のときの小円の直径は比例関係で求めることができる。

include("julia-source.txt")

using SymPy

@syms a::positive, r0::positive, r1::positive,
     r2::positive, x2::negative, y2::negative, r3::positive, x31::negative, x32::negative

r1 = r0//2
eq1 = x31^2 + (r1 - 2r3)^2 - (r1 + r3)^2 |> expand
eq2 = x31^2 + 4r3^2 - (r0 - r3)^2 |> expand
eq3 = x32^2 + r1^2 - (r1 + r3)^2 |> expand
eq4 = x2^2 + (r1 + y2)^2 - (r1 + r2)^2 |> expand
eq5 = (x2 - x32)^2 + y2^2 - (r2 + r3)^2 |> expand
eq6 = x2^2 + y2^2 - (r0 - r2)^2 |> expand;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r2, x2, y2, r3, x31, x32))

   1-element Vector{NTuple{6, Sym}}:
    (3*r0/14, -2*sqrt(6)*r0/7, -5*r0/14, r0/5, -2*sqrt(3)*r0/5, -sqrt(6)*r0/5)

r3/r2 = (1/5)/(3/14) = 14/15 なので， r3 = (14/15)r2 である。
r2 = 3/2 ならば r3 = 7/5 = 1.4。
中円の直径が 3 寸ならば，小円の直径は2寸8分である。

   r2 = 0.214286; x2 = -0.699854; y2 = -0.357143; r3 = 0.2; x31 = -0.69282; x32 = -0.489898
   r0 = 1 のとき，中円の直径 = 0.428571;  小円の直径 = 0.4
   r0 = 1 のとき，中円の直径 = 3//7;  小円の直径 = 2//5

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 1
   r1 = r0//2
   (r2, x2, y2, r3, x31, x32) = (3*r0//14, -2*sqrt(6)*r0//7, -5*r0/14, r0//5, -2*sqrt(3)*r0//5, -sqrt(6)*r0//5)
   @printf("r2 = %g; x2 = %g; y2 = %g; r3 = %g; x31 = %g; x32 = %g\n", r2, x2, y2, r3, x31, x32)
   @printf("r0 = %g のとき，中円の直径 = %g;  小円の直径 = %g\n", r0, 2r2, 2r3)
   println("r0 = $r0 のとき，中円の直径 = $(2r2);  小円の直径 = $(2r3)")
   plot()
   circle(0, 0, r0, :black)
   circle(0, r1, r1)
   circle(0, -r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(x31, 2r3, r3, :orange)
   circle(-x31, 2r3, r3, :orange)
   circle(x32, 0, r3, :orange)
   circle(-x32, 0, r3, :orange)
   segment(-sqrt(r0^2 - r3^2), r3, sqrt(r0^2 - r3^2), r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1, " 大円:r1,(0,r1)", :red, :left, :vcenter)
       point(x2, y2, "中円:r2\n(x2,y2)", :blue, :center, :bottom, delta=delta)
       point(x31, 2r3, "小円:r3\n(x31,2r3)", :black, :center, :top, delta=-delta)
       point(x32, 0, "小円:r3\n(x32,0)", :black, :center, :top, delta=-delta)
       point(0, r3, " r3", :green, :left, :top)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その460）

算額（その460）

埼玉県鴻巣市 薬師堂 明治23年
山口正義(2016): やまぶき2, 第36号
https://yamabukiwasan.sakura.ne.jp/ymbk36.pdf

一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

外円の中に正方形，大円，中円，小円，弦がある。小円の直径が 1 寸のとき，中円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, 2r2 - r0 + r1)
中円の半径と中心座標を r2, (0, 3r2 - r0), (0, r2 - r0)
小円の半径と中心座標を r3, (x3, 2r2 - r0 - r3)
として以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, x1::positive, r2::positive, r3::positive, x3::positive;

eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x1^2 + (2r2 - r0 + r1)^2 - (r0 - r1)^2
eq3 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = x3^2 + (2r2 - r0 - r3)^2 - (r0 - r3)^2
eq5 = distance(0, 4r2 - r0, r0 - 2r2, 2r2, x1, 2r2 - r0 + r1) - r1^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, x1, r2, x3))

   2-element Vector{NTuple{5, Sym}}:
    (25*r3/4, 5*r3/2, 5*sqrt(2)*r3/2, 5*r3/4, sqrt(5)*r3)
    (r3*(59 - 30*sqrt(2))/4, -4*sqrt(2)*r3 + 13*r3/2, -5*r3 + 3*sqrt(2)*r3/2, r3*(2*sqrt(2) + 7)/4, r3*sqrt(2*sqrt(2) + 7))

2 通りの解が得られるが，最初のものが適解である。

   r0 = 3.125;  r1 = 1.25;  x1 = 1.76777;  r2 = 0.625;  x3 = 1.11803
   中円の直径 = 1.25

中円の半径は 5*r3/4 で，小円の半径 r3 の 5/4 倍である。
r3 = 1/2寸のとき，小円の半径は 5/8寸，直径は 5/4寸 = 1寸2分5厘である。

using Plots

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

算額（その459）

算額（その459）

長野県長野市 美和神社 文化10年(1813)

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

上下 2 個の長方形の上段には大円，下段には小円を並べたところ，上下の長方形の横の長さの差が 5 寸になった。大円と小円の個数を求めよ。

算額では大円径，小円径を 91, 75 として解を求めているが，一般解を求めればよい。
どちらの長方形が大きいかの制約条件はない。

術は「剰一術（盈一術），現代解法（中村）は「法による除算」を用いているが，以下に示すコンピュータプログラムで全探索してもたかが知れている。

function prog(r1=91, r2=75; n=100)
   for i in 1:n
       for j in 1:n
           if abs(i*r1 - j*r2) == 5
               println("大円の個数 = $i, 小円の個数 = $j, $(i*r1) - $(j*r2)")
           end
       end
   end
end;

prog()

   大円の個数 = 5, 小円の個数 = 6, 455 - 450
   大円の個数 = 70, 小円の個数 = 85, 6370 - 6375
   大円の個数 = 80, 小円の個数 = 97, 7280 - 7275

円の個数の上限を 100 個までにして探索して複数の答えを得る。
どちらの長方形が大きいかの制約条件はないので，大円 5 個，小円 6 個が最小解である。
術，中村は大円 70 個，小円 85 個としているが，小円を入れる長方形の方が大きい最小の解である。

大円径，小円径が 67, 47 のとき，以下のような解が求まる。

prog(67, 47)

   大円の個数 = 12, 小円の個数 = 17, 804 - 799
   大円の個数 = 35, 小円の個数 = 50, 2345 - 2350
   大円の個数 = 59, 小円の個数 = 84, 3953 - 3948

大円径，小円径が 1009, 997 のときは探索範囲を 1000 個まで広げて，以下のような解が求まる。

prog(1009, 997, n=1000)

   大円の個数 = 415, 小円の個数 = 420, 418735 - 418740
   大円の個数 = 582, 小円の個数 = 589, 587238 - 587233

include("julia-source.txt")

using Plots

function draw(r1, r2, n1, n2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (r1, r2) ./ 2
   (w1, w2) = (n1*2r1, n2*r2)
   plot([0, n1*2r1, n1*2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   for i in 1:n1
       circlef((2i - 1)r1, r1, r1, :cyan)
   end
   plot!([0, n2*2r2, n2*2r2, 0, 0], 2r1 .+ [0, 0, 2r2, 2r2, 0], color=:black, lw=0.5)
   for i in 1:n2
       circlef((2i - 1)r2, r2 + 2r1, r2, :red)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       annotate!(r1, r1, text(@sprintf("%i寸 × %i個 = %i寸", 2r1, n1, n1*2r1), :left, :vcenter, 10))
       annotate!(r1, 2r1 + r2, text(@sprintf("%i寸 × %i個 = %i寸", 2r2, n2, n2*2r2), :left, :vcenter, 10))
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その458）

算額（その458）

長野県長野市 久保寺観音堂 享和3年(1803)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

円に正三角形が内接している。
正三角形の一辺と外積（円と正三角形の間の面積）の和を A，矢（弦の中点と弧の中点を結ぶ線分の長さ）と円周の和を B とするとき，円の面積を求めよ。

最初に，簡単な解を求める方法について書く。条件 A は不要で，かつ，条件 B はとてつもなく簡単なものである。

図において，「矢」は直径の 1/4 である。
したがって，円の半径を r とすれば，「矢と円周の和を B とする」は 「r/2 + 2r*π = B」 である。
半径は，この方程式を r について解けば得られる。

include("julia-source.txt")

using SymPy
@syms r, B
solve(r/2 + 2r*PI - B, r)[1] |> println

   2*B/(1 + 4*pi)

円の半径が 2B/(1 + 4π) なので，面積は π*(2B/(1 + 4π))^2 である。

B = 8.139822368615503
r = 2B/(1 + 4π)
(r, π*(2B/(1 + 4π))^2)

   (1.2, 4.523893421169302)

---

まともにやるととてつもなく面倒なことになる。

円の半径を r，正三角形の一辺の長さを 2a とする。
円，正三角形の面積をそれぞれ s1, s2 とおく。

include("julia-source.txt")

using SymPy

@syms r::positive, a::positive, s1::positive, s2::positive, 矢::positive, A::positive, B::positive
s1 = PI*r^2
a = r*sqrt(Sym(3))/2
s2 = a^2*sqrt(Sym(3))
矢 = r - r/2;

A, B を定義する。

eq1 = 2a + (s1 - s2) - A
eq2 = 矢 + PI*2r - B;

差を取って，r について解くのが常套手段ではあるが，eq2 を見るとこれを r について解けばよいことがわかる。

solve(eq2, r) |> println

   Sym[2*B/(1 + 4*pi)]

ということで，冒頭の簡単な方法へ戻る。

なお，eq1 は不要であるし，差を取ったのでは条件 A が残り，面倒な式になる。無理にやると解は二通り得られるが，どちらが適解かは単純にはいえない。以下に示すように，どちらを取るかは r  に依存するのである。これでは自分の尻尾を飲み込もうとしている蛇と同じである。

eq3 = eq1 - eq2
res = solve(eq3, r);

res[1] |> println

   (-sqrt(-12*sqrt(3)*A + 16*pi*A - 16*pi*B + 12*sqrt(3)*B - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2) - 4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3))

res[2] |> println

   (sqrt(-12*sqrt(3)*A + 16*pi*A - 16*pi*B + 12*sqrt(3)*B - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2) - 4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3))

円の半径を与えて，A, B を計算し，それを eq3 の解 res により円の半径を再計算する関数 rev を書いてみる。
r > (-4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3) の場合には最初のもの，そうでない場合には2番めのものが適解である。

function rev(r)  # 半径を与えて検算する
   s1 = pi*r^2
   a = r*sqrt(3)/2
   s2 = a^2*sqrt(3)
   矢 = r - r/2
   A0 = 2a + (s1 - s2)
   B0 = 矢 + pi*2r
   println(A0, ", ", B0)
   term1 = sqrt(-12*sqrt(3)*A0 + 16*pi*A0 - 16*pi*B0 + 12*sqrt(3)*B0 - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2)
   term2 = - 4*pi - 1 + 2*sqrt(3)
   ans1 = (-term1 + term2)/(-4*pi + 3*sqrt(3))
   ans2 = ( term1 + term2)/(-4*pi + 3*sqrt(3))
   return (r > (-4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3)) ? ans1 : ans2)
end;

半径が 1.2 の場合

rev(1.2)

   4.731739518077568, 8.139822368615503

   1.1999999999999988

半径が 3.45 の場合

rev(3.45)

   27.906580792648718, 23.401989309769576

   3.4499999999999993

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1.2
   a = r*√3/2
   b = -r/2
   c = 3r*(2 - √3)/2
   d = r*(11 - 6*√3)/2
   plot([a, 0, -a, a], [-r/2, r, -r/2, -r/2], color=:blue, lw=0.5)
   circle(0, 0, r)
   segment(0, -r/2, 0, -r, :green, lw=2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r, 0, "r ", :red, :right, :bottom, delta=delta/2)
       point(a, -r/2, "(a,-r/2) ", :blue, :right, :bottom, delta=delta/2)
       point(0, -r/2, " -r/2", :blue, :left, :bottom, delta=delta/2)
       point(0, -3r/4, " 矢", :green, mark=false)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その457）

算額（その457）

長野県長野市 久保寺観音堂 享和3年(1803)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

円に正三角形，正三角形に正方形が内接している。円の直径，正三角形の辺，正方形の辺のそれぞれの 3 乗の和を A，円の直径，正方形の辺の 3 乗の差を B とするとき，円の直径を A, B で表わせ。

include("julia-source.txt")

using SymPy

円の半径，正三角形の辺，正方形の辺をそれぞれ r, 2a, 2c とおく。
正方形の下辺と上辺がy軸と交差するy座標を b, d とする。
a = r√3/2, b = -r/2 から，c, d を求める。

@syms a, b, c, d, r
a = r*√Sym(3)/2
b = -r/2
eq1 = d - b - 2c
eq2 = 2c/(a - c) - sqrt(Sym(3))
solve([eq1, eq2], (c, d))

   Dict{Any, Any} with 2 entries:
     d => -3*sqrt(3)*r + 11*r/2
     c => -3*sqrt(3)*r/2 + 3*r

r, a, c と A, B の関係式を記述する。

@syms r1::positive, r2::positive,
     r3::positive, x3::positive,
     A::positive, B::positive;
a = r*√Sym(3)/2
c = 3r*(2 - √Sym(3))/2
eq1 = (2r)^3 + (2a)^3 + (2c)^3 - A
eq2 = (2r)^3 - (2c)^3 - B;

eq1 + eq2 としてA, B, r の関係式を得る。

eq3 = eq1 + eq2
eq3 |> println

   -A - B + 3*sqrt(3)*r^3 + 16*r^3

eq3 を r について解く（r を A, B で表す式を求める）。

solve(eq3, r)[1] |> factor |> println

   (A + B)^(1/3)/(3*sqrt(3) + 16)^(1/3)

(中身)^(1/3) の形なので，3 乗根の中身を簡約化する。

@syms x, A::positive, B::positive;
apart((A + B)/(3sqrt(Sym(3)) + 16), x) |> simplify |>factor |> println

   -(-16 + 3*sqrt(3))*(A + B)/229

求める関係式は (16 - 3√3)*(A + B)/229 の 3 乗根である。
求まるのは半径なので，2倍して直径にする関数 f を定義する。

f(A, B) = 2(A + B)^(1/3)/(3*sqrt(3) + 16)^(1/3);

検算してみる。

r = 1; a = r*√3/2; c = 3r*(2 - √3)/2 のとき，
A = 13.715575357311328; B = 7.480577065395304
である。

r = 1
a = r*√3/2
c = 3r*(2 - √3)/2
A = (2r)^3 + (2a)^3 + (2c)^3
B = (2r)^3 - (2c)^3
(A, B) |> println

   (13.715575357311328, 7.480577065395304)

計算された A, B から f() によって円の直径を求める。
確かに 2.0 になる。

f(13.715575357311328, 7.480577065395304)

   2.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1
   a = r*√3/2
   b = -r/2
   c = 3r*(2 - √3)/2
   d = r*(11 - 6*√3)/2
   plot([a, 0, -a, a], [-r/2, r, -r/2, -r/2], color=:blue, lw=0.5)
   plot!([c, c, -c, -c, c], [-r/2, d, d, -r/2, -r/2], color=:green, lw=0.5)
   circle(0, 0, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r, 0, "r ", :red, :right, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(0, d, " d", :green, :left, :bottom, delta=delta/2)
       point(c, 0, "c ", :green, :right, :bottom, delta=delta/2)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta/2)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その456）

算額（その456）

(20) 兵庫県青垣町遠坂字後岶 熊野神社 明治18年(1885)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

丹波新聞(2023.03.25): これわかる?　小学生が「和算」挑戦　138年前の難問　神社に掲げられた「算額」
https://tanba.jp/2023/03/%E3%81%93%E3%82%8C%E3%82%8F%E3%81%8B%E3%82%8B%EF%BC%9F%E3%80%80%E5%B0%8F%E5%AD%A6%E7%94%9F%E3%81%8C%E3%80%8C%E5%92%8C%E7%AE%97%E3%80%8D%E6%8C%91%E6%88%A6%E3%80%80138%E5%B9%B4%E5%89%8D%E3%81%AE/

キーワード：円4個，直線上

直線の上に天円，地円，人円がある。それぞれは互いに外接している。地円，人円の直径がそれぞれ　2 寸，3 寸のとき，天円の直径を求めよ。

算額（その82）と同じであった。ただし，こちらのほうが古い。

天円の半径と中心座標を r1, (r1 + 2r2)
地円の半径と中心座標を r2, (0, r2)
人円の半径と中心座標を r3, (x3, r3)
とし，連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r1::positive, r2::positive,
     r3::positive, x3::positive;

eq1 = x3^2 + (r1 + 2r2 - r3)^2 - (r1 + r3)^2
eq2 = x3^2 + (r3 - r2)^2 - (r2 + r3)^2;
res = solve([eq1, eq2], (r1, x3))

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

   r1 = 2;  x3 = 2.44949
   天円の直径 = 4

天円の半径は，地円の半径の2乗を人円と地円の半径の差で割ることにより得られる。

r2 = 2/2 寸, r3 = 3/2 寸のとき，天円の半径は r1 = 2 寸。つまり，直径は 4 寸。

r2 = 2/2
r3 = 3/2
println(-r2^2/(r2 - r3))

   2.0

using Plots

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

算額（その455）

算額（その455）

長野県中野市 田上観音堂 文化6年(1809)

中村信弥「改訂増補 長野県の算額」(p.89)
http://www.wasan.jp/zoho/zoho.html
県内の算額1

鈎股弦内に2個の正方形と，大円，中円，小円が1個ずつ入っている。大円，中円，小円の直径がそれぞれ 280,210,168 寸のとき，弦の長さを求めよ。

結論から先にいうと，小円の直径が 168 寸では，算額に示されているような図が得られない。
そこで，小円の直径は未知ということで解を求める。

補助として，図のように勾股弦に内接する半径 r の円を加える。
大円，中円，小円の半径を r1, r2, r3 として，以下の連立方程式を解く。
なお，一度にすべての解を求めるのは無理なので，鈎，股，弦，r, x2 を先に求め，次いで r3, x1, y1 を求める。

include("julia-source.txt")

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     r::positive,
     x1::positive, x2::positive, y1::positive,
     r1::positive, r2::positive, r3::positive;

(r1, r2) = (280, 210) .// 2
y2 = x2

eq0 = 鈎 + 股 - 弦 - 2r
eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = r*y2 - 鈎*r1
eq3 = r*x2 - 股*r2
eq4 = y2/(股 - x2) - 鈎/股;
res1 = solve([eq0, eq1, eq2, eq3, eq4], (鈎, 股, 弦, r, x2))

   1-element Vector{NTuple{5, Sym}}:
    (735, 980, 1225, 245, 420)

#@syms 鈎::positive, 股::positive, 弦::positive,
     r::positive,
     x1::positive, x2::positive, y1::positive,
     r1::positive, r2::positive, r3::positive;

#(r1, r2) = (280, 210) .// 2
#y2 = x2
(鈎, 股, 弦, r, x2) = (735, 980, 1225, 245, 420)

#eq5 = distance(0, y1, x1, 0, r3, r3) - r3^2
eq5 = sqrt(x1^2 + y1^2)/弦 - r3/r

eq6 = distance(0, 鈎, 股, 0, x1, 0) - (x1^2 + y1^2);
eq7 = y1/x1 - 鈎/股

res2 = solve([eq5, eq6, eq7], (r3, x1, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (79.4594594594595, 317.837837837838, 238.378378378378)

res = solve([eq0, eq1, eq2, eq3, eq4, eq5, eq6, eq7], (鈎, 股, 弦, r, x2, r3, x1, y1))

   1-element Vector{NTuple{8, Sym}}:
    (735.000000000000, 980.000000000000, 1225.00000000000, 245.000000000000, 420.000000000000, 79.4594594594595, 317.837837837838, 238.378378378378)

   鈎 = 735;  股 = 980;  弦 = 1225;  x1 = 317.838;  x2 = 420;  y1 = 238.378
   r3 = 79.4595;  小円の直径 = 158.919
   傾斜している正方形の一辺の長さ = 397.297
   正立している正方形の一辺の長さ = 420

弦の長さは 1225 寸である。算額の答，術，および中村による現代解も 1295 寸を得ている。この解は，小円の直径が 168 寸であるとして導き出されている。しかし，最初に述べたように直径が 168 寸では正方形が弦と接しない。一般的に，算額の解法では目的とする一つのパラメータ（値）だけを求めるので，他のパラメータがどうなるかについては無頓着なので，矛盾に気づかないのであろう。
いずれにしろ，小円の直径は 2✕79.4594594594595 = 158.918918918919 である。そうすれば，図形の相対位置関係は算額の問にある図と同等になる。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (280, 210) .// 2
   (鈎, 股, 弦, r, x2) = res1[1]  # (735, 980, 1225, 245, 420)
   (r3, x1, y1) = res2[1]  # (79.4594594594595, 317.837837837838, 238.378378378378)
   #(鈎, 股, 弦, r, x2, r3, x1, y1) = res[1]
   y2 = x2
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  x1 = %g;  x2 = %g;  y1 = %g\n", 鈎, 股, 弦, x1, x2, y1)
   @printf("r3 = %g;  小円の直径 = %g\n", r3, 2r3)
   @printf("傾斜している正方形の一辺の長さ = %g\n", sqrt(distance(0, 鈎, 股, 0, 0, y1)))
   @printf("正立している正方形の一辺の長さ = %g\n", x2)
   plot([0, 股, 0, 0],[0, 0, 鈎, 0], color=:gray, lw=0.5)
   circle(r, r, r, :gray)
   circle(x2 + r1, r1, r1, :blue)
   circle(r2, y2 + r2, r2, :red)
   circle(r3, r3, r3, :green)
   segment(0, y2, x2, y2)
   segment(x2, 0, x2, y2)
   l = sqrt(x1^2 + y1^2)
   plot!([x1, l*鈎/弦 + x1, l*鈎/弦, 0, x1], [0, l*股/弦, l*股/弦 + y1, y1, 0], color=:orange, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       point(股, 0, " 股", :black , :left, :bottom, delta=delta/2)
       point(x2 + r1, r1, " 大円:r1,(x2+r1,r1)", :black, :left, :bottom, delta=delta)
       point(r2, y2 + r2, "        中円:r2,(r2,y2+r2)", :black, :center, :bottom, delta=5delta)
       point(r3, r3, " 小円:r3,(r3,r3)", :black, :left, :bottom, delta=delta)
       point(x1, 0, " x1", :black, :left, :bottom, delta=delta/2)
       point(x2, 0, " x2", :black, :left, :bottom, delta=delta/2)
       point(0, y1, " y1", :black, :left, :bottom, delta=delta/2)
       point(0, y2, " y2", :black, :left, :bottom, delta=delta/2)
       point(l*鈎/弦, l*股/弦 + y1, " (l*鈎/弦,l*股/弦+y1)", :black, :left, :vcenter)
       point(l*鈎/弦 + x1, l*股/弦, "  (l*鈎/弦+x1,l*股/弦)\n   l=sqrt(x1^2+y1^2)", :black, :left, :vcenter)

       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

モンティ・ホール問題のシミュレーション -- 2

R で書いたシミュレーションを Julia に翻訳した。

https://blog.goo.ne.jp/r-de-r/e/9dbe50704d057552798b9cfbeb0cf982

using Random
using Statistics

# モンティ・ホール問題のシミュレーション

function MontyHallProblem(ChangeAnswer=true)
    door = shuffle(["景品", "ヤギ", "ヤギ"])  # セット
    Player = rand(1:3)[1]  # どれかを選ぶ
    Monty = setdiff(setdiff(1:3, Player), indexin(["景品"], door))[1]  # 必ずヤギのドアを開ける
    if ChangeAnswer  # 最初の答を変えるなら
        Player = setdiff(setdiff(1:3, Player), Monty)[1]
    end
    return door[Player] == "景品"  # 当たったかどうか結果を返す
end

最初の回答を変える場合

@time mean(MontyHallProblem(true) for _ in 1:10000000) |> println

    0.6664597
     12.091644 seconds (320.48 M allocations: 28.194 GiB, 12.49% gc time, 1.82% compilation time)

最初の回答を変えない場合

@time mean(MontyHallProblem(false) for _ in 1:10000000) |> println

    0.3332185
      8.299969 seconds (200.03 M allocations: 18.032 GiB, 12.26% gc time, 0.16% compilation time)

 

Julia: モンティ・ホール問題のシミュレーション

using Random

# モンティ・ホール問題のシミュレーション

function monty_hall_simulation(num_trials::Int)
    stay_wins = 0 # ドアを変えない場合の勝利回数
    switch_wins = 0 # ドアを変えた場合の勝利回数
    
    for _ in 1:num_trials
        # 3つのドアと1つの正解ドアをランダムに配置
        doors = Random.shuffle([false, false, true])
        
        # プレイヤーが最初にランダムに1つのドアを選択
        player_choice = rand(1:3)
        
        # ホストが正解ドア以外かつプレイヤーが選んでいないドアを開ける
        host_choices = [i for i in 1:3 if i != player_choice && !doors[i]]
        host_open = rand(host_choices)
        
        # プレイヤーがドアを変えない場合
        stay_choice = player_choice
        
        # プレイヤーがドアを変える場合
        switch_choice = [i for i in 1:3 if i != player_choice && i != host_open][1]
        
        # 勝敗を判定
        stay_wins += doors[stay_choice]
        switch_wins += doors[switch_choice]
    end
    
    stay_win_percentage = stay_wins / num_trials * 100
    switch_win_percentage = switch_wins / num_trials * 100
    
    println("ドアを変えない場合の勝率: $stay_win_percentage%")
    println("ドアを変えた場合の勝率: $switch_win_percentage%")
end

# シミュレーションを実行
monty_hall_simulation(10000000) # 10000000回試行

    ドアを変えない場合の勝率: 33.33347%
    ドアを変えた場合の勝率: 66.66653%

 

算額（その454）

算額（その454）

数学史研究，通巻 186 号，2005年7月〜9月
http://www.wasan.jp/sugakusipdf/186.pdf

中村正教(1937):「伊額二十二」昭和12年
算額（その298）https://blog.goo.ne.jp/r-de-r/s/%E3%81%9D%E3%81%AE298 より

半円内に二個の等円が入っている。等円の直径を斜の長さで表わせ。

半円の半径と中心座標を r0, (0, 0)
等円の半径と中心座標を r1, (r1, r1)
斜の長さを a
とおく。

a = √2r0 である。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, a::positive;

a = sqrt(Sym(2))r0
eq = 2r1^2 - (r0 - r1)^2
res = solve(eq, r1)[1]
res |> println

   r0*(-1 + sqrt(2))

r0 = a/√2 を代入する。

res(r0 => a/sqrt(Sym(2))) |> simplify |> println

   a*(2 - sqrt(2))/2

すなわち，等円の直径は (2 - √2)a である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   r0 = a/√2
   r1 = (2 - √2)a/2
   plot([-r0, 0, r0], [0, r0, 0], color=:gray, lw=0.5)
   circle(0, 0, r0, beginangle=0, endangle=180)
   circle(r1, r1, r1, :blue)
   circle(-r1, r1, r1, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, r1, "等円:r1,(r1,r1)  ", :blue, :center, :top, delta=-2delta)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;