算額（その653）

算額（その653）

群馬県前橋市河原浜町 大胡神社 大正4年(1915)
http://www.wasan.jp/gunma/ohgo.html

この算額の問・答・術はほとんど判読できない。図はかろうじて見ることができるので，以下のようなものでもあろうかと推測して問を建ててみた。

外円内に互いに外接する 7 個の等円がある。更に外側の等円に外接し外円に内接する 7 個の小円がある。
外円の直径が与えられたとき，等円と小円の直径を求めよ。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

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

   2-element Vector{Tuple{Sym, Sym}}:
    (-R/6 + sqrt(3)*R/3 + (1/2 - 2*sqrt(3)/3)*(10*sqrt(3)*R/39 + 9*R/13 - 4*sqrt(6)*R*sqrt(11894 - 6867*sqrt(3))/(3*(625 - 361*sqrt(3)))), 10*sqrt(3)*R/39 + 9*R/13 - 4*sqrt(6)*R*sqrt(11894 - 6867*sqrt(3))/(3*(625 - 361*sqrt(3))))
    (-R/6 + sqrt(3)*R/3 + (1/2 - 2*sqrt(3)/3)*(4*sqrt(6)*R*sqrt(11894 - 6867*sqrt(3))/(3*(625 - 361*sqrt(3))) + 10*sqrt(3)*R/39 + 9*R/13), 4*sqrt(6)*R*sqrt(11894 - 6867*sqrt(3))/(3*(625 - 361*sqrt(3))) + 10*sqrt(3)*R/39 + 9*R/13)

2 組の解が得られるが，2 番目のものが適解である。
表示される式は長いが簡約化すれば R/3，(5 - 2√3)R/13 になる。

res[2][1] |> sympy.sqrtdenest |> simplify |> println
res[2][2] |> sympy.sqrtdenest |> simplify |> println

   R/3
   -2*sqrt(3)*R/13 + 5*R/13

外円の直径が 3 のとき，等円の直径は 1，小円の直径は 3(5 - 2√3)/13 = 0.3544380888143644 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 3/2
   (r1, r2) = R .* [1/3, (5 - 2√3)/13]
   @printf("R = %g;  r1 = %g;  r2 = %g\n", R, r1, r2)
   plot()
   circle(0, 0, R, :green)
   circle(0, 0, r1)
   rotate(√3r1, r1, r1, angle=60)
   rotate(R - r2, 0, r2, :blue, angle=60)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(√3r1, r1, "等円:r1,(√3r1,r1)", :red, :center, delta=-delta/2)
       point(R - r2, 0, "小円:r2\n(R-r2,0)", :blue, :right, delta=-delta/2)
   end
end;

算額（その652）

算額（その652）

群馬県吾妻郡東吾妻町 一宮神社 明治5年(1872)3月
http://www.wasan.jp/gunma/itinomiya.html

この算額の問・答・術はほとんど判読できない。図はかろうじて見ることができるので，以下のようなものでもあろうかと推測して問を建ててみた。

正三角形と甲円，乙円，丙円がある。甲円，乙円はそれぞれ正三角形の斜辺に 2 点で接しており，甲円と乙円は外接している。丙円は甲円と正三角形の底辺が作る弧の中にあり，底辺と甲円に接している。
正三角形の一辺の長さと甲円の直径がわかっているときに，乙円の直径を求めよ。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, r1::positive, y1::positive, r2::positive, r3::positive
eq1 = (y1 - r1) + 2r3
eq2 = dist(0, sqrt(Sym(3))a, a, 0, 0, y1) - r1^2
eq3 = r2/(sqrt(Sym(3))a - y1 - r1 - r2) - 1//2
res = solve([eq1, eq2, eq3], (y1, r2, r3))

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

乙円は甲円の 1/3 の大きさである。
甲円の直径が 150 のとき，乙円の直径は 50 である。

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

a = 100;  r1 = 75;  y1 = 23.2051;  r2 = 25;  r3 = 25.8975

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 100
   r1 = 75
   (y1, r2, r3) = (√3a - 2r1, r1/3, (3r1 - √3a)/2)
   @printf("a = %g;  r1 = %g;  y1 = %g;  r2 = %g;  r3 = %g\n", a, r1, y1, r2, r3)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   circle(0, y1, r1)
   circle(0, -r3, r3, :blue)
   circle(0, r2 + r1 + y1, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, y1, " 甲円:r1,(0, y1)", :red, :left, :vcenter)
       point(0, y1 + r1 + r2, " 乙円:r2,(0,y1+r1+r2)", :green, :left, :vcenter)
       point(0, -r3, " 丙円:r3,(0,-r3)", :blue, :left, :vcenter)
       point(0, √3a, " √3a", :black, :left, :vcenter)
       point(a, 0, "a", :black, :left, :bottom, delta=delta/2)
   end
end;

算額（その651）

算額（その651）

長野県飯山市静間大久保 静間観音堂（静観庵）弘化5年(1848)

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

日堂薫，宮崎雄也，鷲森勇希，伊藤栄一：「飯山市静間観音堂の算額」の復元
https://sbc21.co.jp/gakkokagaku/2015/27.pdf

直径と2本の弦を隔てて 4 個の等円が外円に内接し，直径と弦に接している。等円の直径を外円の直径で表わせ。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r, x, y1, x0, y2
y0 = sqrt(R^2 - x0^2)
eq1 = x^2 + y2^2 - (R - r)^2
eq2 = r^2 + y1^2 - (R - r)^2
eq3 = dist(0, R, x0, -y0, r, y1) - r^2
eq4 = dist(0, R, x0, -y0, x, y2) - r^2
eq5 = y2/x * (y0 + R)/x0 - 1;

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)
   (r, x, y1, x0, y2) = u
   t = -R - sqrt(R^2 - x0^2)
   t2 = x0^2 + t^2
   return [
       x^2 + y2^2 - (R - r)^2,  # eq1
       r^2 + y1^2 - (R - r)^2,  # eq2
       -r^2 + (r - x0*(r*x0 + (-R + y1)*t)/t2)^2 + (-R + y1 - t*(r*x0 + (-R + y1)*t)/t2)^2,  # eq3
       -r^2 + (x - x0*(x*x0 + (-R + y2)*t)/t2)^2 + (-R + y2 - t*(x*x0 + (-R + y2)*t)/t2)^2,  # eq4
       -1 + y2*(R + sqrt(R^2 - x0^2))/(x*x0),  # eq5
   ]
end;
R = 1
iniv = BigFloat[8, 17, -18, 18, 8] ./26
iniv = BigFloat[0.313, 0.637, -0.612, 0.694, 0.257] .* R
res = nls(H, ini=iniv)

   (BigFloat[0.3129063194259740325268103642250497724124515879919851568449522877018154669364085, 0.6371784719628400853411527943465112615089968502655608047796456430352258943786332, -0.6117085589952554644133700611927907136554589335804161767166315804026984221703721, 0.6940076375173001867788131238285698449294652216735132602273056924441179784109038, 0.2571017711954972908126786596219933186205074882118679191481568858007998140638422], true)

等円の直径は外円の直径の 0.312906 倍である。

   r = 0.3129063194259740325268103642250497724124515879919851568449522877018154669364085
   R = 1;  r = 0.312906;  x = 0.637178;  y1 = -0.611709;  y2 = 0.257102;  x0 = 0.694008

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x, y1, x0, y2) = res[1]
   println("r = $r")
   @printf("R = %g;  r = %g;  x = %g;  y1 = %g;  y2 = %g;  x0 = %g\n", R, r, x, y1, y2, x0)
   y0 = -sqrt(R^2 - x0^2) 
   plot([x0, 0, -x0], [y0, R, y0], color=:black, lw=0.5)
   circle(0, 0, R)
   circle(x, y2, r, :blue)
   circle(-x, y2, r, :blue)
   circle(r, y1, r, :blue)
   circle(-r, y1, r, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(x, y2, "(x,y2)")
       point(r, y1, "(r,y1)")
       point(x0, y0, " (x0,y0)", :red)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
   end
end;

『「飯山市静間観音堂の算額」の復元』には，記号を上に合わせると，r について以下の 3 次式が示されている。
R = 1 にして方程式を解くと，実数解 1 つが得られ，
-(1133/8 + 15*sqrt(114))^(1/3)/3 + 71/(12*(1133/8 + 15*sqrt(114))^(1/3)) + 5/3
評価すると上と同じく r = 0.312906319425974 が得られる。

@syms R, r
eq = 4r^3 - 20R*r^2 + 57R^2*r - 16R^3;

res2 = solve(eq(R=>1));
res2[3] |> println

   -(1133/8 + 15*sqrt(114))^(1/3)/3 + 71/(12*(1133/8 + 15*sqrt(114))^(1/3)) + 5/3

res2[3].evalf() |> println

   0.312906319425974