算額（その680）

算額（その680）

長野県信州新町日名 日置神社 明治44年(1911)
中村信弥「改訂増補 長野県の算額」県内の算額(270)
http://www.wasan.jp/zoho/zoho.html

第7問　正方形内に大円と小円がある正方形の面積が 235平方寸，外積（正方形と二円の間の面積）が 132.45平方寸，大小の円の直径の差が 2 寸のとき，小円の直径を求めよ。ただし円積率を 0.7 とする（円周率を 4*0.7 とするということ）。

正方形の一辺の長さを a, 二円の半径を r1, r2 とし，以下の連立方程式を解く。

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

using SymPy
@syms r1::positive, r2::positive, 円周率::positive,
     a::positive
円周率 = 4 * 0.7
a = √235
eq2 = a^2 - 円周率*(r1^2 + r2^2) - 132.45
eq3 = 2r1 - 2r2 - 2;
res = solve([eq2, eq3], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (4.75000000000000, 3.75000000000000)

小円の直径は 2*3.75 = 7.5 寸である。
この設定で図を描いても，2 つの円は外接しない。

---

円周率として π を使い，正方形内の2円が外接するとして，外積がいくつになるか計算してみる。

using SymPy
@syms r1::positive, r2::positive, 円周率::positive,
     a::positive
円周率 = PI
a = √Sym(235)
eq2 = √Sym(2)*(a - r2 - r1) - (r1 + r2) 
eq3 = 2r1 - 2r2 - 2;
res = solve([eq2, eq3], (r1, r2))

   Dict{Any, Any} with 2 entries:
     r2 => -sqrt(470)/2 - 1/2 + sqrt(235)
     r1 => -sqrt(470)/2 + 1/2 + sqrt(235)

r1 = -sqrt(470)/2 + 1/2 + sqrt(235)
r2 = -sqrt(470)/2 - 1/2 + sqrt(235)
外積 = (a^2 - 円周率*(r1^2 + r2^2))
(r1, r2, 外積.evalf())

   (4.989968022416491, 3.989968022416491, 106.761363826833)

外積はほぼ 106.76 である。この結果だと，ちゃんと二円は外接する。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1, r2) = (√235, -sqrt(470)/2 + 1/2 + sqrt(235), -sqrt(470)/2 - 1/2 + sqrt(235))
   @printf("a = %g;  r1 = %g;  r2 = %g\n", a, r1, r2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circlef(r1, r1, r1, :orange)
   circlef(a - r2, a - r2, r2, :tomato)
   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)
   end
end;

算額（その679）

算額（その679）

長野県飯山市下木島 鳥出神社 天保14年(1843)
中村信弥「改訂増補 長野県の算額」県内の算額（166）
http://www.wasan.jp/zoho/zoho.html

三角形内に三辺に接する全円と，中鈎（頂点から底辺への垂線）と二辺に接する中円と小円が入っている。全円，中円，小円の直径がそれぞれ 7寸，6寸，4寸のとき，中鈎の長さ（三角形の高さ）はいかほどか。

中鈎と底辺の交点を原点と定め，右の頂点と左の頂点の座標を (a, 0), (-b, 0)，高さを h とする。
全円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (r2, r2)
小円の半径と中心座標を r3, (-r3, r3)
とおき，以下の連立方程式を解く。
4 元連立方程式として解くと不適切な解になるので，eq1, eq2, eq3 から a, b, h を求め，そのあと eq4 から x1 を求める。

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

using SymPy
@syms a::positive, b::negative, h::positive,
     r1::positive, x1::positive,
     r2::positive, r3::positive, l1::positive, l2::positive
(r1, r2, r3) = (7, 6, 4) .// 2
l1 = sqrt(a^2 + h^2)
l2 = sqrt(b^2 + h^2)
eq1 = (a - b + l1 + l2)*r1 - (a - b)*h
eq2 = a + h - l1 - 2r2
eq3 = -b + h - l2 - 2r3;
res1 = solve([eq1, eq2, eq3], (a, b, h))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (15, -6, 8)

中鈎の長さ（三角形の高さ）は 8 寸である。

@syms a::positive, b::negative, h::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, r3::positive, l1::positive, l2::positive
(r1, r2, r3) = (7, 6, 4) .// 2
(a, b, h) = (15, -6, 8)
eq4 = dist(a, 0, 0, h, x1, r1) - r1^2
res2 = solve(eq4, x1)
res2 |> println

   Sym[0.999999999999994, 15.8750000000000]

x1 = 1 が適解である。

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

   r1 = 3.5;  r2 = 3;  r3 = 2;  a = 15;  b = -6;  h = 8;  x1 = 1

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (7, 6, 4) .// 2
   (a, b, h) = res1[1]
   x1 = res2[1]
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  a = %g;  b = %g;  h = %g;  x1 = %g\n",
       r1, r2, r3, a, b, h, x1)
   plot([a, 0, b, a], [0, h, 0, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(-r3, r3, r3, :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(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(b, 0, "b", :black, :right, :bottom, delta=delta)
       point(0, h, " h", :black, :left, :bottom, delta=delta)
       point(x1, r1, "    全円:r1\n    (x1,r1)", :red, :center, :bottom, delta=delta)
       point(r2, r2, "中円:r2\n(r2,r2)", :blue, :center, :top, delta=-delta)
       point(-r3, r3, "小円:r3\n(-r3,r3)", :green, :center, :top, delta=-delta)
   end
end;

 

算額（その678）

算額（その678）

長野県小諸市諸 金毘羅社 文政7年(1824)
中村信弥「改訂増補 長野県の算額」県内の算額（106）
http://www.wasan.jp/zoho/zoho.html

円内に等楕円が 3 個入っている。円の直径が 99 寸のとき，楕円の面積が最大になるときの楕円の長径を求めよ。

円の半径と中心座標を R, (0, 0)
楕円の長径，短径と中心座標を a, b, (0, y)
円と楕円の共通接点を (x1, y1)
上と右の楕円の共通接点を (x2, y2)
とおき，以下の連立方程式を解く。

eq1, eq2, eq3 は円と楕円が共通接点 (x1, y1) を持つこと，
eq4, eq5 は 2 つの楕円が共通接点 (x2, y2) を持つことを意味している。
これらの方程式群から，楕円の中心の y 座標を求める。

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

using SymPy
@syms x1::positive, y1::positive, x2::positive, y2::positive, y::positive, a::positive, b::positive, R::positive
y2 = x2/sqrt(Sym(3))
eq1 = x1^2/a^2 + (y1 - y)^2/b^2 - 1
eq2 = b^2*x1/(a^2*(y - y1)) + x1/y1
eq3 = x1^2 + y1^2 - R^2
eq4 = x2^2/a^2 + (y2 - y)^2/b^2 - 1
eq5 = b^2*x2/(a^2*(y - y2))- 1/sqrt(Sym(3));

res1 = solve([eq1, eq2, eq3], (x1, y1, y))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt((-R^2*b^2 + a^4)/(a - b))/sqrt(a + b), a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2), sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a)
    (-sqrt(-(R^2*b^2 - a^4)/(a - b))/sqrt(a + b), a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2), sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a)

res2 = solve([eq4, eq5], (x2, y))

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

両方の解 y が等しいとすると以下のようになる。

eq = res1[1][3] - res2[1][2]
eq |> println

   -sqrt(3*a^2 + 9*b^2)/3 + sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a

方程式を b について解く。

res3 = solve(eq,  b)
res3[1] |> println

   a*sqrt(9*R^2 - 12*a^2)/(3*R)

楕円の面積 S は長軸の長さ a の関数になり，以下のような外形である。

S = PI * a * a*sqrt(9*R^2 - 12*a^2)/(3*R);

円の半径を R = 99/2 として，図を描いてみる。
a が 30〜40 の間で面積が最大になることがわかる。

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(S(R => 99/2), xlims=(0, 43), xlabel="長径 a", ylabel="楕円の面積 S")

S の導関数を求め，接線の傾きが 0 になるときの a を求める。

max_at_a = solve(diff(S, a), a)
max_at_a[1] |> println

   sqrt(2)*R/2

a の値が 35.0017856687341 のとき，面積が最大値 8888.52378442978 になる。

max_at_a[1](R => 99/2).evalf() |> println

   35.0017856687341

S(R => 99/2, a => 35.0017856687341).evalf() |> println

   2222.13094610745

そのときの短径 b は 20.2082903779612 である。

(a*sqrt(9*R^2 - 12*a^2)/(3*R))(R => 99/2, a => 35.0017856687341).evalf() |> println

   20.2082903779612

算額において，長径，短径は差渡し径なので，現代の用語で表すものの 2 倍である。
よって，算額の答えとしては「長径は 70.0035713374682，短径は 40.4165807559224」である。

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

   x1 = 24.75;  y1 = 42.8683;  x2 = 24.75; y2 = 14.2894;  y = 28.5788;  a = 35.0018;  b = 20.2083

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 99/2
   a = √2R/2
   b = a*sqrt(9*R^2 - 12*a^2)/(3*R)
   (x1, y1, y) = (sqrt((-R^2*b^2 + a^4)/(a - b))/sqrt(a + b),
       a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2),
       sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a)
   x2 = a^2/sqrt(a^2 + 3*b^2)
   y2 = x2/sqrt(3)
   @printf("x1 = %g;  y1 = %g;  x2 = %g; y2 = %g;  y = %g;  a = %g;  b = %g\n", x1, y1, x2, y2, y, a, b)
   plot()
   circle(0, 0, R)
   ellipse(0, y, a, b, color=:blue)
   ellipse(√3y/2, -y/2, a, b, φ=240, color=:blue)
   ellipse(-√3y/2, -y/2, a, b, φ=300, color=: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(x1, y1, " (x1,y1)", :blue, :left, :bottom)
       point(x2, y2, " (x2,y2)", :blue, :left)
       point(0, y, " y", :blue, :left, :vcenter)
   end
end;

算額（その677）

算額（その677）

群馬の算額 148-7 赤城神社 明治20年
http://takasakiwasan.web.fc2.com/gunnsann/g148−7.html

菱形の中に等円が内接する小さな菱形を 12 個入れる。（大きな）菱形の一辺の長さが 20 寸，等円の直径が 4.8 寸のとき，（大きな）菱形の長い対角線の長さを求めよ。

小さな菱形の長い対角線を 2a, 短い対角線を 2b とする。小さい菱形の一辺の長さを c とする。
等円の半径を r とする。

以下の連立方程式を解き，a, b を求めれば，（大きな）菱形の長い対角線の長さは 2 * 4a である。

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

using SymPy
@syms a::positive, b::positive, r::positive, c::positive

eq1 = r/a - b/c
eq2 = a^2 + b^2 - c^2
res = solve([eq1, eq2], (a, b))

   2-element Vector{Tuple{Sym, Sym}}:
    (sqrt(2)*sqrt(c)*sqrt(c - sqrt(c - 2*r)*sqrt(c + 2*r))*(c + sqrt(c^2 - 4*r^2))/(4*r), sqrt(2)*sqrt(c)*sqrt(c - sqrt(c - 2*r)*sqrt(c + 2*r))/2)
    (sqrt(2)*sqrt(c)*sqrt(c + sqrt(c - 2*r)*sqrt(c + 2*r))*(c - sqrt(c^2 - 4*r^2))/(4*r), sqrt(2)*sqrt(c)*sqrt(c + sqrt(c - 2*r)*sqrt(c + 2*r))/2)

2 組の解が得られるが，a > b なので最初のものが適解である。

b = sqrt(2c)*sqrt(c - sqrt(c^2 - 4r^2))/2
a = b*(c + sqrt(c^2 - 4r^2))/2r;

r = 4.8, c = 5 のとき，a = 4, b = 3 である。

a(r => 48/20, c => 5).evalf() |> println
b(r => 48/20, c => 5).evalf() |> println

   4.00000000000000
   3.00000000000000

（大きな）菱形の長い対角線の長さは 2 * 4a なので 32 寸である。

function unit(x, y, a, b, r)
   circle(x, y, r)
   plot!([x, x + a, x, x - a, x], [y - b, y, y + b, y, y - b], color=:blue, lw=0.5)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 4 .* [0, 1, 2, 3, 2, 1, 0, -1, -2, -3, -2, -1]
   y = 3 .* [6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5]
   (a, b, r) = (4, 3, 4.8/2)
   plot()
   for i = 1:12
       unit(x[i], y[i], a, b, r)
   end
   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(a, 0, " a", :blue, :left, delta=-delta/2)
       point(0, b, " b", :blue, :left, :vcenter)
   end
end;

算額（その676）

算額（その676）

長野県長野市松代町西条 西楽寺 弘化5年(1848)
中村信弥「改訂増補 長野県の算額」県内の算額（176）
http://www.wasan.jp/zoho/zoho.html

甲円は外円に内接し，乙円，丙円...は隣の円と甲円に外接し，外円に内接している。外円，甲円，乙円の直径がそれぞれ 168 寸，88 寸，3 寸のとき，乙円から始まる累円のうち最も大きい円の直径を求めよ。

甲円，乙円のパラメータを初期値として，丙円，丁円...のパラメータは順次計算できる。
デカルトの円定理で直径だけを求めることはできるが，図を描くためには累円の中心座標も必要なので，逐次的に数値解を求める。

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

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, x2::positive, y2::negativez

(R, r1, r2) = (168, 88, 3) .// 2
eq1 = x2^2 + (r1 - R - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
res1 = solve([eq1, eq2], (x2, y2))

   1-element Vector{Tuple{Sym, Sym}}:
    (231/10, -396/5)

現在の円のパラメータ（初期状態では r2, x2, y2）から，次円のパラメータを求める関数は以下のようになる。
解は 2 通り求まる。一方は累円が大きくなる方向，もう一方は累円が小さくなる方向である。
累円が一番大きくなると解は振動する。
次に，最大の累円のパラメータを再設定して今度は累円が小さくなる方向で累円のパラメータを求める。

@syms R::positive, r1::positive, r2::positive, x2::positive, y2
@syms r3, x3, y3
function nextcircle(R, r1, r2, x2, y2)
   eq3 = x3^2 + (r1 - R - y3)^2 - (r1 + r3)^2
   eq4 = x3^2 + y3^2 - (R - r3)^2
   eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
   res2 = solve([eq3, eq4, eq5], (r3, x3, y3))
end;

実際に計算すると，累円が最大になるのは甲円から数えて 10 番目の円である。
円の半径は 77/2，中心座標は (-231/10, 196/5) である。

n = 2  (3//2, 231//10, -396//5) ## 乙円
n = 3  (231/118, 15477/590, -22932/295)
n = 4  (77/29, 4389/145, -10948/145)
n = 5  (231/61, 10857/305, -21924/305)
n = 6  (231/40, 8547/200, -1638/25)
n = 7  (77/8, 2079/40, -266/5)
n = 8  (231/13, 3927/65, -1764/65)
n = 9  (33, 231/5, 108/5)
n = 10  (77/2, -231/10, 196/5)
n = 11  (231/10, -3003/50, -252/25)
n = 12  (231/19, -5313/95, -4284/95)
n = 13  (7, -231/5, -308/5)
n = 14  (231/52, -9933/260, -4536/65)
n = 15  (231/76, -12243/380, -7056/95)
n = 16  (11/5, -693/25, -1924/25)
n = 17  (231/139, -16863/695, -54684/695)
n = 18  (231/178, -19173/890, -35532/445)
n = 19  (77/74, -7161/370, -14924/185)
n = 20  (231/271, -23793/1355, -110124/1355)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (168, 88) .// 2
   plot()
   circle(0, 0, R)
   circlef(0, r1 - R, r1, :green)
   point(0, r1 - R, " 1", :white, :left, :vcenter, mark=false)
   (r2, x2, y2) = (3//2, 231//10, -396//5)
   for n = 2:10
       println("n = $n  ", (r2, x2, y2))
       circlef(x2, y2, r2, n)
       point(x2, y2, string(n), :white, :center, :vcenter, mark=false)
       (r2, x2, y2) = nextcircle(168//2, 88//2, r2, x2, y2)[2]
   end
   (r2, x2, y2) = (77//2, -231//10, 196//5)
   for n = 11:20
       (r2, x2, y2) = nextcircle(168//2, 88//2, r2, x2, y2)[1]
       point(x2, y2, string(n), :white, :center, :vcenter, mark=false)
       println("n = $n  ", (r2, x2, y2))
       circlef(x2, y2, r2, n)
   end
end;

算額（その675）

算額（その675）

新潟県長岡市与板町本与板 与板八幡宮 文化元年(1804)
涌田和芳，外川一仁: 与板八幡宮の紛失算額（2），長岡工業高等専門学校研究紀要，第48巻，2012.
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_46-50/vol_48/48_07wakuta.pdf

三角形内に全円が内接している。三角形の辺の長い方から，大斜，中斜，小斜と名付ける。
2 本の斜線を入れ，区画された領域に甲円，乙円を入れる。甲円，乙円は斜線と三角形の辺に接している。
三角形の高さ（中鈎と呼ぶ）が 3 寸，全円の直径が 2 寸，乙円の直径が 1 寸のとき，甲円の直径はいかほどか。

大斜の長さを a，中鉤と大斜の交点を (b, 0)，中鈎の長さを h，斜線と中斜，小斜の交点座標を (c, cy), (d, dy) とする。
全円の半径と中心座標を r1, (x1, r1)
甲円の半径と中心座標を r2, (x2, r2);  x2 = x1 である。
乙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。
なお，x2 = x1, cy = h*(a - c)/(a - b), dy = h*d/b である。

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

using SymPy
@syms a, b, c, d, h, cy, dy, r1, x1, y1, r2, x2, y2, r3, x3, y3
x2 = x1
cy = h*(a - c)/(a - b)
dy = h*d/b
eq1 = dist(b, h, a, 0, x1, r1) - r1^2
eq2 = dist(0, 0, c, cy, x2, r2) - r2^2
eq3 = dist(0, 0, b, h, x3, y3) - r3^2
eq4 = dist(0, 0, b, h, x1, r1) - r1^2
eq5 = dist(a, 0, d, dy, x3, y3) - r3^2
eq6 = dist(a, 0, d, dy, x2, r2) - r2^2
eq7 = (a + sqrt((a - b)^2 + h^2) + sqrt(b^2 + h^2))*r1 - a*h
eq8 = (h - y3)/(x3 - b) - (h - r1)/(x1 - b);

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, r2, x3, y3, a, b, c, d) = u
   return [
       -r1^2 + (-b + x1 - (a - b)*(-h*(-h + r1) + (a - b)*(-b + x1))/(h^2 + (a - b)^2))^2 + (-h + h*(-h*(-h + r1) + (a - b)*(-b + x1))/(h^2 + (a - b)^2) + r1)^2,  # eq1
       -r2^2 + (-c*(c*x1 + h*r2*(a - c)/(a - b))/(c^2 + h^2*(a - c)^2/(a - b)^2) + x1)^2 + (-h*(a - c)*(c*x1 + h*r2*(a - c)/(a - b))/((a - b)*(c^2 + h^2*(a - c)^2/(a - b)^2)) + r2)^2,  # eq2
       -r3^2 + (-b*(b*x3 + h*y3)/(b^2 + h^2) + x3)^2 + (-h*(b*x3 + h*y3)/(b^2 + h^2) + y3)^2,  # eq3
       -r1^2 + (-b*(b*x1 + h*r1)/(b^2 + h^2) + x1)^2 + (-h*(b*x1 + h*r1)/(b^2 + h^2) + r1)^2,  # eq4
       -r3^2 + (y3 - d*h*((-a + d)*(-a + x3) + d*h*y3/b)/(b*((-a + d)^2 + d^2*h^2/b^2)))^2 + (-a + x3 - (-a + d)*((-a + d)*(-a + x3) + d*h*y3/b)/((-a + d)^2 + d^2*h^2/b^2))^2,  # eq5
       -r2^2 + (r2 - d*h*((-a + d)*(-a + x1) + d*h*r2/b)/(b*((-a + d)^2 + d^2*h^2/b^2)))^2 + (-a + x1 - (-a + d)*((-a + d)*(-a + x1) + d*h*r2/b)/((-a + d)^2 + d^2*h^2/b^2))^2,  # eq6
       -a*h + r1*(a + sqrt(b^2 + h^2) + sqrt(h^2 + (a - b)^2)),  # eq7
       (h - y3)/(-b + x3) - (h - r1)/(-b + x1),  # eq8
   ]
end;

(h, r1, r3) = (3, 2//2, 1//2)
iniv = BigFloat[1.6, 0.6, 1.53, 2, 3.47, 1.47, 2.2, 0.9]
res = nls(H, ini=iniv)

   (BigFloat[1.608783345088622385337178167670691516795402078864423523680572389575119930785486, 0.6000000000000000000000000000000000000000000000000000000000000000000814904650529, 1.544788378049916344806170511615455296433424123818581654504997356980354393989754, 2.000000000000000000000000000000000000000000000000000000000000000000199561380369, 3.473546558332068932798386959562327915038715977912214524063444909529663746416152, 1.480793411011210304275162855560219076071446168772739785329422324385407988362897, 2.204699202026184681037300584878993043245344709257443697861753272879706442814819, 0.9087534413959433109983030676937383181821297022187909525620532161139324795037814], true)

   甲円の直径 = 1.2
   x1 = 1.60878;  r2 = 0.6;  x2 = 1.60878;  x3 = 1.54479;  y3 = 2;  a = 3.47355;  b = 1.48079;  c = 2.2047;  cy = 1.91019;  d = 0.908753;  dy = 1.84108

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, r2, x3, y3, a, b, c, d) = res[1]
   x2 = x1
   cy = h*(a - c)/(a - b)
   dy = h*d/b
   @printf("甲円の直径 = %g\n", 2r2)
   @printf("x1 = %g;  r2 = %g;  x2 = %g;  x3 = %g;  y3 = %g;  a = %g;  b = %g;  c = %g;  cy = %g;  d = %g;  dy = %g\n", x1, r2, x2, x3, y3, a, b, c, cy, d, dy)
   x2 = x1
   cy = h*(a - c)/(a - b)
   dy = h*d/b
   plot([0, a, b, 0], [0, 0, h, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :blue)
   circle(x2, r2, r2, :green)
   circle(x3, y3, r3)
   segment(0, 0, c, cy)
   segment(a, 0, d, dy)
   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)
       plot!(ylims=(-0.2, 3.1))
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(c, cy, " (c,cy)", :black, :left, :bottom)
       point(d, dy, "(d,dy) ", :black, :right, :bottom)
       point(b, h, "(b,h)", :black, :center, :bottom, delta=delta)
       point(x1, r1, "全円:r1,(x1,r1)", :blue, :center, delta=-delta)
       point(x2, r2, "甲円:r2,(x1,r2)", :green, :center, delta=-delta)
       point(x3, y3, "乙円:r3,(x3,y3)", :red, :center, :bottom, delta=delta)
   end
end;

算額（その674）

算額（その674）

『不朽算法』上巻第六問 [日下貞八郎誠嗣辺・安島万蔵直円遺稿]
米山忠興: 側円術(IV) ---『安子遺稿側円解二条第一』の現代語訳 ---，東洋大学紀要. 自然科学篇，Vol. 52, pp. 97-116, 2008-03
http://id.nii.ac.jp/1060/00002535/

直角三角形（鈎股弦）内に円と楕円が入っている。円は直角三角形の三辺に接しており，楕円と外接している。楕円の長径は底辺（股）に平行で，直角三角形の斜辺（弦）に接している。
直角三角形の辺の長さと楕円の短径を与えて，長径を求めよ。

直角三角形の直角を挟む二辺のうち，短い方を「鈎」，長い方を「股」とする（斜辺「弦」は sqrt(鈎^2 + 股^2) である）。
円の半径 r は「(鈎 + 股 - 弦)/2」，中心座標は (r, r)
楕円の長径と短径を 2a, 2b，中心座標を (x0, y0) とする。
楕円と斜辺，楕円と円の接点座標を (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms 鈎, 股, 弦, a, b, r, x, x1, y1, x2, y2, x0, y0, slope0

(鈎, 股, b) = (2240, 2808, 455//2)
弦 = 3592  # sqrt(鈎^2 + 股^2)
r = (鈎 + 股 - 弦)//2
y0 = b
長径，短径が a b，中心座標が (x0, y0) の楕円
(x1, y1) が楕円上にある
eq1 = (x1 - x0)^2/a^2 + (y1 - y0)^2/b^2 - 1
(x1, y1) における楕円の接線の傾き
slope0 = -鈎//股  # 目的とする接線。符号に注意
eq2 = b^2*(x1 - x0)/(a^2*(y0 - y1)) - slope0  # slope
もう一個の条件
eq3 = y1/(股 - x1) + slope0
楕円と円が外接する
eq4 = (x2 - x0)^2/a^2 + (y2 - y0)^2/b^2 - 1
slope2 = (x2 - r)/(r - y2)
eq5 = b^2*(x2 - x0)/(a^2*(y0 - y2)) - slope2
eq6 = (x2 - r)^2 + (y2 - r)^2 - r^2;

using NLsolve

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

function H(u)
   (x1, y1, x0, x2, y2, a) = u
   return [
       4*(y1 - 455/2)^2/207025 - 1 + (-x0 + x1)^2/a^2,  # eq1
       280/351 + (-207025*x0/4 + 207025*x1/4)/(a^2*(455/2 - y1)),  # eq2
       y1/(2808 - x1) - 280/351,  # eq3
       4*(y2 - 455/2)^2/207025 - 1 + (-x0 + x2)^2/a^2,  # eq4
       -(x2 - 728)/(728 - y2) + (-207025*x0/4 + 207025*x2/4)/(a^2*(455/2 - y2)),  # eq5
       (x2 - 728)^2 + (y2 - 728)^2 - 529984,  # eq6
   ]
end;

iniv = BigFloat[2398, 327, 1870, 1300, 230, 580]
res = nls(H, ini=iniv)

   (BigFloat[2397.842696629213493816709019056387960138168817704006327516310149222101255756478, 327.1910112359550540245449460804349431777480338181361396304024428466875180291161, 1872.000000000000004115470053719782306992126301163385447666330986804201857308351, 1310.399999999999999733344750238321544467653247568361438815996098107649959134677, 291.1999999999999996444596669844287261829629732073189591407411915504426299859508, 585.0000000000000042869479726247732354323812232955671058370412692641114894978273], true)

楕円の長径(a) = 585（差渡し径(2a) = 1170）である。

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

   鈎 = 2240;  股 = 2808;  b = 227.5
   x1 = 2397.84;  y1 = 327.191;  x0 = 1872;  y0 = 227.5
   x2 = 1310.4;  y2 = 291.2;  r = 728

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, b) = (2240, 2808, 455//2)
   弦 = sqrt(鈎^2 + 股^2)
   r = (鈎 + 股 - 弦)/2
   y0 = b
   (x1, y1, x0, x2, y2, a) = res[1]
   @printf("鈎 = %g;  股 = %g;  b = %g\n", 鈎, 股, b)
   @printf("楕円の長径(a) = %g　（差渡し径(2a) = %g\n", a, 2a)
   @printf("x1 = %g;  y1 = %g;  x0 = %g;  y0 = %g\n", x1, y1, x0, y0)
   @printf("x2 = %g;  y2 = %g;  r = %g\n", x2, y2, r)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   ellipse(x0, b, a, b, color=:red)
   circle(r, r, 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(x0, y0, "(x0,y0)")
       point(r, r, "(r,r)")
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
       point(x1, y1, "(x1,y1) ", :red, :right)
       point(x2, y2, "(x2,y2) ", :red, :right)
   end
end;

 

算額（その673）

算額（その673）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問題. 100　19✕19魔方陣

関孝和が考案した方式があるらしいが，多分追跡できないので，ヒンズー方式，バシェー方式の二通りで書いてみた。

function generate_odd_magic_square(n)
   # ヒンズーの連続方式
   m = zeros(Int, n, n)
   (i, j) = (div(n + 1, 2), 1)
   m[j, i] = 1
   for num in 2:n^2
       (next_i, next_j) = (i + 1, j - 1)
       next_i > n && (next_i = 1)
       next_j < 1 && (next_j = n)
       m[next_j, next_i] != 0 && ((next_i, next_j) = (i, j + 1))
       m[next_j, next_i] = num
       (i, j) = (next_i, next_j)
   end
   return m
end;

function generate_odd_magic_square2(n)
   m = zeros(Int, n, n)
   (i, j) = (n - 1, div(n - 1, 2))
   for num in 1:n^2
       m[i + 1, j + 1] = num
       (next_i, next_j) = ((i + 1) % n + 1, (j + 1) % n + 1)
       if m[next_i, next_j] == 0
           (i, j) = (next_i - 1, next_j - 1)
       else
           i -= 1
       end
   end
   return m
end;

function generate_odd_magic_square3(n)
   # バシェー方式
   a = zeros(Int, n, n)
   m = n ÷ 2
   for i = 1:n
       (sx, sy) = (i - 1 - m, n - i - m)
       for j = 1:n
           (x, y) = (sx + j, sy + j)
           x <= 0 && (x += n)
           x >  n && (x -= n)
           y <= 0 && (y += n)
           y >  n && (y -= n)
           a[x, y] = (i - 1) * n + j
       end
   end
   return a
end;

function print_magic_square(m)
   n = size(m, 1)
   for i in 1:n
       for j in 1:n
           print(m[i, j], "\t")
       end
       println()
   end
end;

function check_magic_square(magic_square)
   n = size(magic_square, 1)
   sums = n*(n^2 + 1)/2
   allequal(sum(magic_square, dims=1) .== sums) &&
   allequal(sum(magic_square, dims=2) .== sums) &&
   sum(magic_square[i, i] for i in 1:n) == sums &&
   sum(magic_square[i, n + 1 -i] for i in 1:n) == sums
end;

n = 19
magic_square = generate_odd_magic_square(n);
magic_square2 = generate_odd_magic_square2(n);
magic_square3 = generate_odd_magic_square3(n);

print_magic_square(magic_square);
println("--")
print_magic_square(magic_square2);
println("--")
print_magic_square(magic_square3);
println("--")

check_magic_square(magic_square) |> println
check_magic_square(magic_square2) |> println
check_magic_square(magic_square3) |> println

   192 213 234 255 276 297 318 339 360 1 22 43 64 85 106 127 148 169 190 
   212 233 254 275 296 317 338 359 19 21 42 63 84 105 126 147 168 189 191 
   232 253 274 295 316 337 358 18 20 41 62 83 104 125 146 167 188 209 211 
   252 273 294 315 336 357 17 38 40 61 82 103 124 145 166 187 208 210 231 
   272 293 314 335 356 16 37 39 60 81 102 123 144 165 186 207 228 230 251 
   292 313 334 355 15 36 57 59 80 101 122 143 164 185 206 227 229 250 271 
   312 333 354 14 35 56 58 79 100 121 142 163 184 205 226 247 249 270 291 
   332 353 13 34 55 76 78 99 120 141 162 183 204 225 246 248 269 290 311 
   352 12 33 54 75 77 98 119 140 161 182 203 224 245 266 268 289 310 331 
   11 32 53 74 95 97 118 139 160 181 202 223 244 265 267 288 309 330 351 
   31 52 73 94 96 117 138 159 180 201 222 243 264 285 287 308 329 350 10 
   51 72 93 114 116 137 158 179 200 221 242 263 284 286 307 328 349 9 30 
   71 92 113 115 136 157 178 199 220 241 262 283 304 306 327 348 8 29 50 
   91 112 133 135 156 177 198 219 240 261 282 303 305 326 347 7 28 49 70 
   111 132 134 155 176 197 218 239 260 281 302 323 325 346 6 27 48 69 90 
   131 152 154 175 196 217 238 259 280 301 322 324 345 5 26 47 68 89 110 
   151 153 174 195 216 237 258 279 300 321 342 344 4 25 46 67 88 109 130 
   171 173 194 215 236 257 278 299 320 341 343 3 24 45 66 87 108 129 150 
   172 193 214 235 256 277 298 319 340 361 2 23 44 65 86 107 128 149 170 
   --
   172 193 214 235 256 277 298 319 340 361 2 23 44 65 86 107 128 149 170 
   171 173 194 215 236 257 278 299 320 341 343 3 24 45 66 87 108 129 150 
   151 153 174 195 216 237 258 279 300 321 342 344 4 25 46 67 88 109 130 
   131 152 154 175 196 217 238 259 280 301 322 324 345 5 26 47 68 89 110 
   111 132 134 155 176 197 218 239 260 281 302 323 325 346 6 27 48 69 90 
   91 112 133 135 156 177 198 219 240 261 282 303 305 326 347 7 28 49 70 
   71 92 113 115 136 157 178 199 220 241 262 283 304 306 327 348 8 29 50 
   51 72 93 114 116 137 158 179 200 221 242 263 284 286 307 328 349 9 30 
   31 52 73 94 96 117 138 159 180 201 222 243 264 285 287 308 329 350 10 
   11 32 53 74 95 97 118 139 160 181 202 223 244 265 267 288 309 330 351 
   352 12 33 54 75 77 98 119 140 161 182 203 224 245 266 268 289 310 331 
   332 353 13 34 55 76 78 99 120 141 162 183 204 225 246 248 269 290 311 
   312 333 354 14 35 56 58 79 100 121 142 163 184 205 226 247 249 270 291 
   292 313 334 355 15 36 57 59 80 101 122 143 164 185 206 227 229 250 271 
   272 293 314 335 356 16 37 39 60 81 102 123 144 165 186 207 228 230 251 
   252 273 294 315 336 357 17 38 40 61 82 103 124 145 166 187 208 210 231 
   232 253 274 295 316 337 358 18 20 41 62 83 104 125 146 167 188 209 211 
   212 233 254 275 296 317 338 359 19 21 42 63 84 105 126 147 168 189 191 
   192 213 234 255 276 297 318 339 360 1 22 43 64 85 106 127 148 169 190 
   --
   172 353 154 335 136 317 118 299 100 281 82 263 64 245 46 227 28 209 10 
   11 173 354 155 336 137 318 119 300 101 282 83 264 65 246 47 228 29 191 
   192 12 174 355 156 337 138 319 120 301 102 283 84 265 66 247 48 210 30 
   31 193 13 175 356 157 338 139 320 121 302 103 284 85 266 67 229 49 211 
   212 32 194 14 176 357 158 339 140 321 122 303 104 285 86 248 68 230 50 
   51 213 33 195 15 177 358 159 340 141 322 123 304 105 267 87 249 69 231 
   232 52 214 34 196 16 178 359 160 341 142 323 124 286 106 268 88 250 70 
   71 233 53 215 35 197 17 179 360 161 342 143 305 125 287 107 269 89 251 
   252 72 234 54 216 36 198 18 180 361 162 324 144 306 126 288 108 270 90 
   91 253 73 235 55 217 37 199 19 181 343 163 325 145 307 127 289 109 271 
   272 92 254 74 236 56 218 38 200 1 182 344 164 326 146 308 128 290 110 
   111 273 93 255 75 237 57 219 20 201 2 183 345 165 327 147 309 129 291 
   292 112 274 94 256 76 238 39 220 21 202 3 184 346 166 328 148 310 130 
   131 293 113 275 95 257 58 239 40 221 22 203 4 185 347 167 329 149 311 
   312 132 294 114 276 77 258 59 240 41 222 23 204 5 186 348 168 330 150 
   151 313 133 295 96 277 78 259 60 241 42 223 24 205 6 187 349 169 331 
   332 152 314 115 296 97 278 79 260 61 242 43 224 25 206 7 188 350 170 
   171 333 134 315 116 297 98 279 80 261 62 243 44 225 26 207 8 189 351 
   352 153 334 135 316 117 298 99 280 81 262 63 244 45 226 27 208 9 190 
   --
   true
   true
   true

 

算額（その672）

算額（その672）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問. 61　大円の内外に中円，小円を 19 個配置する。大円の直径が 1 尺のとき，中円，小円の直径はいかほどか。

プログラム的には外円に内接する小円を求めるのが簡単である。
外円（半径 R，中心座標 (0, 0)）に外接する中円の半径，内接する小円の半径をそれぞれ r1, r2 とする。以下の方程式を解く。

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

using SymPy
@syms θ::positive, R::positive, r1::positive, r2::positive, r3::positive, n::integer

θ = 360/n
eq1 = (R + r1)*sind(θ/2) - r1
eq2 = (R - r2)*sind(θ/2) - r2
solve([eq1, eq2], (r1, r2)) |> println

   Dict{Any, Any}(r2 => R*sin(pi/n)/(sin(pi/n) + 1), r1 => -R*sin(pi/n)/(sin(pi/n) - 1))

r1 = R*sin(pi/n)/(1 - sin(pi/n))
r2 = R*sin(pi/n)/(1 + sin(pi/n))
である。
外円の直径が = 1 尺 = 10 寸のとき，中円の直径は 1.9702361077126025 寸，小円の直径は 1.413320924331764 寸である。

R = 10/2
n = 19
r1 = R*sin(pi/n)/(1 - sin(pi/n))
r2 = R*sin(pi/n)/(1 + sin(pi/n))
(r1, r2) .* 2 |> println

   (1.9702361077126025, 1.413320924331764)

同じプログラムで，外円の直径が R, 円の数が n の図形を順次縮小して図を描くことができる(maxn=2 を指定する）。

function draw(R =10, n=19, maxn=2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   θ = 360/n
   x = cosd.(range(0, 360, n + 1))
   y = sind.(range(0, 360, n + 1))
   r1 = -R*sin(pi/n)/(sin(pi/n) - 1)
   R += 2r1
   plot()
   for j = 1:maxn
       if maxn == 2 || j != 1
           println("外円の直径 = $R,  外円に内接する円の直径 = $r1")
           for i in 1:n
               circle((R - r1)*x[i], (R - r1)*y[i], r1, :blue)
           end
       end
       maxn == 2 && j == 2 && break
       R -= 2r1
       r1 = R*sin(pi/n)/(sin(pi/n) + 1)
       circle(0, 0, R)
   end
   if more && maxn == 2
       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(R, 0, " R", :red, :left, :vcenter)
       point(R + 1.97, 0, " R+r1", :red, :left, :vcenter)
       point(R - 1.41, 0, "R-r2 ", :red, :right, :vcenter)
   end
end;

算額（その671）

算額（その671）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問. 25　正三角形に内接する大小 2 個の円がある。小円の直径が 1 寸のとき，正三角形の一辺の長さはいかほどか。

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

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

using SymPy
@syms r1::positive, r2::positive, a::positive
s3 = sqrt(Sym(3))
eq1 = r2/(s3*a - (2r1 + r2)) - 1//2
eq2 = r1/(s3*a - r1) - 1//2
res = solve([eq1, eq2], (r1, a))

   Dict{Any, Any} with 2 entries:
     r1 => 3*r2
     a  => 3*sqrt(3)*r2

正三角形の一辺の長さは小円の半径の 6√3 倍である。
小円の直径が 1 寸のとき，三角形の一辺の長さは 1/2 * 6√3 = 5.196152422706632 寸である。
ちなみに，そのときの大円の直径は小円の直径の 3 倍なので，3 寸である。

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

 

算額（その670）

算額（その670）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf

米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

岡山県瀬戸内市長船町土師宮森 片山日子神社 明治6年(1873)
http://www.wasan.jp/okayama/katayamahiko.html
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 145，森北出版株式会社，2010年４月22日.

キーワード：円4個，菱形
#Julia, #SymPy, #算額, #和算, #数学

問. 19　菱形に大小各 2 個の円が内接している，菱形の長径が 24 寸，短径が 18 寸のとき，大小の円の直径はいくらか。

算額（その669）では大円の中心が y 軸上にあるが，本問では大円の中心が x 軸上にある。

長径，短径を 2a, 2b とおく。
大円の半径と中心座標を r1, (r1, 0)
小円の半径と中心座標を r2, (0, y2)
とおき，以下の連立方程式を解く。
菱形の一辺の長さは c = sqrt(a^2 + b^2) を使わないで，未定義のままにしておく。そうしないと　r2, y2 の式に何度も出てきて式がとんでもなく長くなる。

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

using SymPy
@syms r1::positive, r2::positive, y2::positive,
     a::positive, b::positive, c::positive
c = sqrt(a^2 + b^2)
eq1 = r1^2 + y2^2 - (r1 + r2)^2
eq2 = r1/(a - r1) - b/c
eq3 = r2/(b - y2) - a/c
res = solve([eq1, eq2, eq3], (r1, r2, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (a*b/(b + c), a*b*(-a^2 + a*sqrt(a^2 + b^2 + 4*b*c + 3*c^2) - b*c - c^2)/(a^2*b + a^2*c - b*c^2 - c^3), a*b*(a*b + 2*a*c - c*sqrt(a^2 + b^2 + 4*b*c + 3*c^2))/(a^2*b + a^2*c - b*c^2 - c^3))
    (a*b/(b + c), -a*b*(a^2 + a*sqrt(a^2 + b^2 + 4*b*c + 3*c^2) + b*c + c^2)/(a^2*b + a^2*c - b*c^2 - c^3), a*b*(a*b + 2*a*c + c*sqrt(a^2 + b^2 + 4*b*c + 3*c^2))/(a^2*b + a^2*c - b*c^2 - c^3))

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

   長径，短径が 24 寸，18 寸のとき，大円の直径は 9 寸，小円の直径は 5.40356 寸である。
   r1 = 4.5;  r2 = 2.70178;  y2 = 5.62278

function draw(n, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (24, 18) .// 2
   c = sqrt(a^2 + b^2)
   d = sqrt(a^2 + b^2 + 4*b*c + 3*c^2)
   e = a^2*b + a^2*c - b*c^2 - c^3
   (r1, r2, y2) = (
       a*b/(b + c),
       a*b*(-a^2 + a*d - b*c - c^2)/e,
       a*b*(a*b + 2*a*c - c*d)/e)
   @printf("大円の直径 = %g;  小円の直径 = %g\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  y2 = %g\n", r1, r2, y2)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:black, lw=0.5)
   circle(r1, 0, r1)
   circle(-r1, 0, r1)
   circle(0, y2, r2, :blue)
   circle(0, -y2, r2, :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(r1, 0, "大円:r1,(r1,0)", :blue, :center, delta=-delta)
       point(0, y2, " 小円:r2\n (0,y2)", :red, :left, :vcenter)
       point(a, 0, "a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
   end
end;

算額（その669）

算額（その669）

関孝和: 闕擬抄一百問答術
http://hyonemitsu.web.fc2.com/Ketsugishyotojutsu.pdf
米光丁: 和算への旅
http://hyonemitsu.web.fc2.com/

問. 18　菱形に大小各 2 個の円が内接している，菱形の長径が 24 寸，短径が 18 寸のとき，大小の円の直径はいくらか。

算額（その670）では大円の中心が x 軸上にあるが，本問では大円の中心が y 軸上にある。

長径，短径を 2a, 2b とおく。
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (x2, 0)
とおき，以下の連立方程式を解く。
菱形の一辺の長さは c = sqrt(a^2 + b^2) を使わないで，未定義のままにしておく。そうしないと　r2, x2 の式に何度も出てきて式がとんでもなく長くなる。

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

using SymPy
@syms r1::positive, r2::positive, x2::positive,
     a::positive, b::positive, c::positive
c = sqrt(a^2 + b^2)
eq1 = x2^2 + r1^2 - (r1 + r2)^2
eq2 = r1/(b - r1) - a/c
eq3 = r2/(a - x2) - b/c
res = solve([eq1, eq2, eq3], (r1, r2, x2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (a*b/(a + c), a*b*(-a*c - b^2 + b*sqrt(a^2 + 4*a*c + b^2 + 3*c^2) - c^2)/(a*b^2 - a*c^2 + b^2*c - c^3), a*b*(a*b + 2*b*c - c*sqrt(a^2 + 4*a*c + b^2 + 3*c^2))/(a*b^2 - a*c^2 + b^2*c - c^3))
    (a*b/(a + c), -a*b*(a*c + b^2 + b*sqrt(a^2 + 4*a*c + b^2 + 3*c^2) + c^2)/(a*b^2 - a*c^2 + b^2*c - c^3), a*b*(a*b + 2*b*c + c*sqrt(a^2 + 4*a*c + b^2 + 3*c^2))/(a*b^2 - a*c^2 + b^2*c - c^3))

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

   長径，短径が 24 寸，18 寸のとき，大円の直径は 8 寸，小円の直径は 6.87539 寸である。
   r1 = 4;  r2 = 3.43769;  x2 = 6.27051

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (24, 18) .// 2
   c = sqrt(a^2 + b^2)
   d = sqrt(a^2 + 4*a*c + b^2 + 3*c^2)
   e = a*b^2 - a*c^2 + b^2*c - c^3
   (r1, r2, x2) = (
       a*b/(a + c),
       a*b*(-a*c - b^2 + b*d - c^2)/e,
       a*b*(a*b + 2*b*c - c*d)/e)
   @printf("大円の直径 = %g;  小円の直径 = %g\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g\n", r1, r2, x2)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, -r1, r1)
   circle(x2, 0, r2, :blue)
   circle(-x2, 0, r2, :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(0, r1, " 大円:r1\n (0,r1)", :red, :left, :vcenter)
       point(x2, 0, "小円:r2,(x2,0)", :blue, :center, delta=-delta)
       point(a, 0, "a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
   end
end;

算額（その668）

算額（その668）

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

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

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

正方形の中に，四分円，半円（大円），天円，累円（甲円，乙円，丙円 ...）があり，天円，累円は隣り合う円と外接しさらに正方形の一辺と半円にも外接している。さらに，全円，中円，小円と斜線があり互いに接している（斜線と接しているのは全円と右側の小円のみ）。
累円の個数 n を，全円，中円，小円の直径と n 番目の累円の直径であらわせ。
注：半円と小円は外接していない。

正方形の一辺の長さを a，斜線と正方形の辺の交点座標を (0, b), (c, 0)
全円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (r2, r2), (3r2 - 2r3, r2), (5r2 - 4r3, r2)
小円の半径と中心座標を r3, (2r2 - r3, r2), (4r2 - 3r3, r2)
とおく。

順次以下の手順にしたがって，解を求める。

まず，以下の連立方程式を解き a, b, c を求める（r1, r2, r3 を含む式で表される）。

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

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, r2::positive, r3::positive,
     x4::positive, r4::positive
eq1 = r2/(a - (5r2 - 4r3)) - r1/(a - r1)
eq2 = r2/(a - c) - b/a
eq3 = a + b - sqrt(a^2 + b^2) - 2r1
res = solve([eq1, eq2, eq3], (a, b, c))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*r1*(r2 - r3)/(r1 - r2), r1*(r1 - 5*r2 + 4*r3)/(r1 - 3*r2 + 2*r3), 4*(r1^2*r2 - r1^2*r3 - 6*r1*r2^2 + 10*r1*r2*r3 - 4*r1*r3^2 + 3*r2^3 - 5*r2^2*r3 + 2*r2*r3^2)/(r1^2 - 6*r1*r2 + 4*r1*r3 + 5*r2^2 - 4*r2*r3))

正方形の一辺の長さは 4r1*(r2 - r3)/(r1 - r2) である。

次に，大円の半径と中心座標を r4, (x4, a - r4) とおき，以下の連立方程式を解く。

eq4 = (a - x4)^2 + r4^2 - (a - r4)^2
eq5 = (a - x4)^2 + (a/2 - r4)^2 - (a/2 + r4)^2
res2 = solve([eq4, eq5], (r4, x4))

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

2 組の解が得られるが，最初のものが適解である。
天円の半径は，正方形の一辺の長さの 1/4 である。

次に累円間の関係式を求める。
天円を初円(x4, a - r4) とし，次の甲円 (x5, a - r5) を決める。

@syms r5, x5
eq6 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
eq7 = (a - x5)^2 + (a/2 - r5)^2 - (a/2 + r5)^2
solve([eq6, eq7], (r5, x5))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-a + x4)*(-2*sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) + a + 2*a*(-2*r4 + x4)/(-a + 2*r4) + x4)/(2*(-a + 2*r4)), sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) - a*(-2*r4 + x4)/(-a + 2*r4))
    ((-a + x4)*(2*sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) + a + 2*a*(-2*r4 + x4)/(-a + 2*r4) + x4)/(2*(-a + 2*r4)), -sqrt(2)*sqrt(a)*sqrt(r4)*(-a + x4)/(-a + 2*r4) - a*(-2*r4 + x4)/(-a + 2*r4))

2 組の解が得られるが，最初のものが適解である。
この解は，左側にある円に接する円のパラメータなので，これにより求められた円とその右側にある円とも同じ関係が成り立つ。

以下の関数は現在の累円から次の累円を求める関数である。

nextcircle(a, r4, x4) = ((x4 - a)*(-2sqrt(2a*r4)*(x4 - a)/(2r4 - a) + a + 2a*(x4 - 2r4)/(2r4 - a) + x4)/(2(2r4 - a)),
       sqrt(2a*r4)*(x4 - a)/(2r4 - a) - a*(x4 - 2*r4)/(2r4 - a))

ただし，このまま長々しい数式のままでの取り扱いは難しい。

直線の上に互いに接する3個の円の半径についてはデカルトの円定理より，以下のような関係式が成り立つ。

大円・天円・甲円の半径においては，1/√甲円の半径 = 1/√大円の半径 + 1/√天円の半径
大円・甲円・乙円の半径においては，1/√乙円の半径 = 1/√大円の半径 + 1/√甲円の半径
大円・乙円・丙円の半径においては，1/√丙円の半径 = 1/√大円の半径 + 1/√乙円の半径
大円・丙円・丁円の半径においては，1/√丁円の半径 = 1/√大円の半径 + 1/√丙円の半径
　：
下から順次代入していくと，甲円から数えて4番目の丁円の半径について
1/√丁円の半径 = 1/√大円の半径 + (1/√大円の半径 + 1/√乙円の半径)
   = 2/√大円の半径 + 1/√乙円の半径 = 2/√大円の半径 + (1/√大円の半径 + 1/√甲円の半径)
   = 3/√大円の半径 + 1/√甲円の半径 = 3/√大円の半径 + (1/√大円の半径 + 1/√天円の半径)
   = 4/√大円の半径 + 1/√天円の半径
甲円から数えて n 番目の累円の半径 rn については
1/√(rn) = n/√大円の半径 + 1/√天円の半径
が成り立つ。

大円，天の半径は上述の通り，a/2, a/4 なので
1/√(rn) = n/√(a/2) + 1/√(a/4)

@syms rn, n, r1, r2, r3
a = 4*r1*(r2 - r3)/(r1 - r2)
eq = 1/√rn - (n/√(a/2) + 1/√(a/4))

方程式を解いて n を求める。

solve(eq, n)[1] |> println

   sqrt(2)*(-sqrt(rn) + sqrt(r1*(r2 - r3)/(r1 - r2)))/sqrt(rn)

SymPyではこれ以上簡約化できないが，以下のように簡約化できる。
sqrt(2)*(sqrt(r1*(r2 - r3)/(r1 - r2)/(rn)) - 1)

実行例を見れば，19 番目の累円の半径は 0.015874124835590503，でこのとき n は 19 となる。

(r1, r2, r3) = (4.3, 3, 1) ./ 2
rn = 0.015874124835590503
sqrt(2)*(sqrt(r1*(r2 - r3)/(r1 - r2)/(rn)) - 1)

   19.000000000000043

   r1 = 2.15;  r2 = 1.5;  r3 = 0.5;  a = 13.2308;  b = 5.33519;  c = 9.51091
   大円の半径 = 6.615384615384616
   天円の半径 = 3.307692307692308
   -19.619822485207095
   n = 0;  n = 0.0;  r = 3.307692307692308;  x = 3.875202587377986
   n = 1;  n = 1.0000000000000009;  r = 1.1350205593713572;  x = 7.750405174755974
   n = 2;  n = 1.9999999999999998;  r = 0.5675102796856791;  x = 9.355566643391244
   n = 3;  n = 2.9999999999999996;  r = 0.33950675324251667;  x = 10.233458601408486
   n = 4;  n = 3.999999999999998;  r = 0.22567545882547532;  x = 10.787058971033915
   n = 5;  n = 4.999999999999996;  r = 0.16079341811242404;  x = 11.168042584375158
          :
   n = 18;  n = 18.00000000000003;  r = 0.01755155072579467;  x = 12.549270122486426
   n = 19;  n = 19.000000000000043;  r = 0.015874124835590503;  x = 12.5826536817502

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (4.3, 3, 1) ./ 2
   (a, b, c) = (4*r1*(r2 - r3)/(r1 - r2), r1*(r1 - 5*r2 + 4*r3)/(r1 - 3*r2 + 2*r3), 4*(r1^2*r2 - r1^2*r3 - 6*r1*r2^2 + 10*r1*r2*r3 - 4*r1*r3^2 + 3*r2^3 - 5*r2^2*r3 + 2*r2*r3^2)/(r1^2 - 6*r1*r2 + 4*r1*r3 + 5*r2^2 - 4*r2*r3))
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  a = %g;  b = %g;  c = %g\n",
       r1, r2, r3, a, b, c)
   println("大円の半径 = $(a/2)")
   println("天円の半径 = $(a/4)")
   println((a - (5r2 - 4r3))^2 + (a - r2)^2 - (a + r2)^2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(3r2 - 2r3, r2, r2, :blue)
   circle(5r2 - 4r3, r2, r2, :blue)
   circle(2r2 - r3, r2, r3, :green)
   circle(4r2 - 3r3, r2, r3, :green)
   segment(0, b, a, 0, :orange)
   circle(a, a, a, :magenta, beginangle=180, endangle=270)
   circle(a, a/2, a/2, :blue, beginangle=90, endangle=270)
   (r4, x4) = (a/4, a*(1 - sqrt(2)/2))
   for i = 1:20
       n = sqrt(2)*(-sqrt(r4) + sqrt(r1*(r2 - r3)/(r1 - r2)))/sqrt(r4)
       println("n = $(round(Int, n));  n = $n;  r = $r4;  x = $(x4)")
       circle(x4, a - r4, r4)
       (r4, x4) = nextcircle(a, r4, x4)
   end
   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(a, a/2, "(a,a/2)")
       point(r1, r1, "全円:r1,(r1,r1)", :black, :center, :bottom, delta=delta)
       point(5r2 - 4r3, r2, "中円:r2,(5r2-4r3,r2)", :black, delta=-delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom)
       point(c, 0, " c", :black, :left, :bottom, delta=delta/2)
       segment(0, r2, c, r2, :gray80)
       segment(c, 0, c, r2, :gray80)
   end
end;