算額（その2117）

算額（その2117）

百四十七 群馬県甘楽郡妙義町下高田 高太神社 大正12年(1923)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に大円 2 個，小円 2 個を容れる。台形の高さが 20 寸，大円の直径が 14.4 寸のとき，小円の直径はいかほどか。

等脚台形の上底，下底の長さをそれぞれ 2b, 2a 
等脚台形の斜辺を延長してできる二等辺三角形の高さを h0
大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (r2, h - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive, h0::positive,
     r1::positive, r2::positive, y2::positive;
eq1 = r1/(h0 - r1) - r2/(h0 - (h - r2))
eq2 = b/(h0 - h) - a/h0
eq3 = dist2(0, h0, a, 0, r1, r1, r1)
eq4 = (r1 - r2)^2 + (h - r2 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (r2, a, b, h0))[2]

   ((-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2), 4*(h^(17/2)*r1^(3/2) - 26*h^(15/2)*r1^(5/2) + 249*h^(13/2)*r1^(7/2) - 1192*h^(11/2)*r1^(9/2) + 3160*h^(9/2)*r1^(11/2) - 4768*h^(7/2)*r1^(13/2) + 3984*h^(5/2)*r1^(15/2) - 1664*h^(3/2)*r1^(17/2) + 256*sqrt(h)*r1^(19/2) - h^8*r1^2 + 26*h^7*r1^3 - 249*h^6*r1^4 + 1192*h^5*r1^5 - 3160*h^4*r1^6 + 4768*h^3*r1^7 - 3984*h^2*r1^8 + 1664*h*r1^9 - 256*r1^10)/(h^9 - 28*h^8*r1 + 301*h^7*r1^2 - 1690*h^6*r1^3 + 5544*h^5*r1^4 - 11088*h^4*r1^5 + 13520*h^3*r1^6 - 9632*h^2*r1^7 + 3584*h*r1^8 - 512*r1^9), 4*(-h^(9/2)*sqrt(r1) + 13*h^(7/2)*r1^(3/2) - 52*h^(5/2)*r1^(5/2) + 68*h^(3/2)*r1^(7/2) - 16*sqrt(h)*r1^(9/2) - h^4*r1 + h^3*r1^2 + 24*h^2*r1^3 - 52*h*r1^4 + 16*r1^5)/(h^4 - 18*h^3*r1 + 84*h^2*r1^2 - 120*h*r1^3 + 32*r1^4), (sqrt(h)*r1*(h - 6*r1) + 2*r1^(3/2)*(-h + 2*r1))/(sqrt(h)*(h - 4*r1)))

res[1] |> println

   (-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2)

小円の半径 r2 は，h, r1 を用いてい以下のように計算できる。

f(h, r1) = (-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2)
f(20, 14.4/2) * 2

   6.399999999999852

台形の高さが 20 寸，大円の直径が 14.4 寸のとき，小円の直径は 6.4 寸である。

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

   h = 20;  r1 = 7.2;  r2 = 3.2;  a = 24.6857;  b = 4.51765;  h0 = 24.48

function draw(h, r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, a, b, h0) = ((-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2), 4*(h^(17/2)*r1^(3/2) - 26*h^(15/2)*r1^(5/2) + 249*h^(13/2)*r1^(7/2) - 1192*h^(11/2)*r1^(9/2) + 3160*h^(9/2)*r1^(11/2) - 4768*h^(7/2)*r1^(13/2) + 3984*h^(5/2)*r1^(15/2) - 1664*h^(3/2)*r1^(17/2) + 256*sqrt(h)*r1^(19/2) - h^8*r1^2 + 26*h^7*r1^3 - 249*h^6*r1^4 + 1192*h^5*r1^5 - 3160*h^4*r1^6 + 4768*h^3*r1^7 - 3984*h^2*r1^8 + 1664*h*r1^9 - 256*r1^10)/(h^9 - 28*h^8*r1 + 301*h^7*r1^2 - 1690*h^6*r1^3 + 5544*h^5*r1^4 - 11088*h^4*r1^5 + 13520*h^3*r1^6 - 9632*h^2*r1^7 + 3584*h*r1^8 - 512*r1^9), 4*(-h^(9/2)*sqrt(r1) + 13*h^(7/2)*r1^(3/2) - 52*h^(5/2)*r1^(5/2) + 68*h^(3/2)*r1^(7/2) - 16*sqrt(h)*r1^(9/2) - h^4*r1 + h^3*r1^2 + 24*h^2*r1^3 - 52*h*r1^4 + 16*r1^5)/(h^4 - 18*h^3*r1 + 84*h^2*r1^2 - 120*h*r1^3 + 32*r1^4), (sqrt(h)*r1*(h - 6*r1) + 2*r1^(3/2)*(-h + 2*r1))/(sqrt(h)*(h - 4*r1)))
   @printf("台形の高さが %g，大円の直径が %g のとき，小円の直径は %g である。\n", h, 2r1, 2r2)
   @printf("h = %g;  r1 = %g;  r2 = %g;  a = %g;  b = %g;  h0 = %g\n", h, r1, r2, a, b, h0)
   plot([a, 0, -a, a], [0, h0, 0, 0], color=:green, lw=0.5)
   segment(-b, h, b, h, :green)
   circle2(r1, r1, r1)
   circle2(r2, h - r2, r2, :blue)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(r2, h - r2, " 小円:r2,(r2,h-r2)", :blue, :left, :vcenter)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(0, h, " h", :green, :left, :bottom, delta=delta/2)
       point(0, h0, " h0", :green, :left, :bottom, delta=delta/2)
       point(b, h, " (b,h)", :green, :left, :vcenter)
   end
end;

draw(20, 14.4/2, true)

 

算額（その2116）

算額（その2116）

百四十七 群馬県甘楽郡妙義町下高田 高太神社 大正12年(1923)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，外円，正三角形
#Julia, #SymPy, #算額, #和算

外円の中に大円 2 個，小円 1 個，直角三角形 1 個を容れる。正三角形の一辺の長さが 7.5 寸，小円の直径が 4.755 寸のとき，大円の直径はいかほどか。

正三角形の一辺の長さを 2a，頂点と円の接点座標を (a, x0, y0);x0 = a, y0 = -sqrt(R^2 - a^2)
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, y1); x1 = r1, y1 = √Sym(3)a + y0
小円の半径と中心座標を r2, (0, R - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, R::positive, r1::positive,
     x1::positive, y1::positive, r2::positive;
x0 = a
y0 = -sqrt(R^2 - a^2)
x1 = r1
y1 = √Sym(3)a + y0
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = R - 2r2 - √Sym(3)a - y0
res = solve([eq1, eq2], (r1, R))[1]

   (a*(-3*a^3 - 4*sqrt(3)*a^2*r2 - 4*a*r2^2 + 3*a*sqrt(a^4 + 4*sqrt(3)*a^3*r2 + 16*a^2*r2^2 + 8*sqrt(3)*a*r2^3 + 4*r2^4) + 2*r2*sqrt(3*a^4 + 12*sqrt(3)*a^3*r2 + 48*a^2*r2^2 + 24*sqrt(3)*a*r2^3 + 12*r2^4))/(2*(sqrt(3)*a^3 + 5*a^2*r2 + 3*sqrt(3)*a*r2^2 + 2*r2^3)), 2*(a^2 + sqrt(3)*a*r2 + r2^2)/(sqrt(3)*a + 2*r2))

正三角形の一辺の長さが 7.5 寸，小円の直径が 4.755 寸のとき，大円の直径は 5.89244501321687 寸である。

「答」には「大円径六寸二分五厘」と書いているが，それは多分，外円の半径である。

res[1](a => 7.5/2, r2 => 4.755/2).evalf() * 2 |> println

   5.89244501321687

f(a, r2) = a*(-3*a^3 - 4*sqrt(3)*a^2*r2 - 4*a*r2^2 + 3*a*sqrt(a^4 + 4*sqrt(3)*a^3*r2 + 16*a^2*r2^2 + 8*sqrt(3)*a*r2^3 + 4*r2^4) + 2*r2*sqrt(3*a^4 + 12*sqrt(3)*a^3*r2 + 48*a^2*r2^2 + 24*sqrt(3)*a*r2^3 + 12*r2^4))/(2*(sqrt(3)*a^3 + 5*a^2*r2 + 3*sqrt(3)*a*r2^2 + 2*r2^3))
f(7.5/2, 4.755/2) * 2

   5.892445013216866

function draw(a, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, R) = (a*(-3*a^3 - 4*sqrt(3)*a^2*r2 - 4*a*r2^2 + 3*a*sqrt(a^4 + 4*sqrt(3)*a^3*r2 + 16*a^2*r2^2 + 8*sqrt(3)*a*r2^3 + 4*r2^4) + 2*r2*sqrt(3*a^4 + 12*sqrt(3)*a^3*r2 + 48*a^2*r2^2 + 24*sqrt(3)*a*r2^3 + 12*r2^4))/(2*(sqrt(3)*a^3 + 5*a^2*r2 + 3*sqrt(3)*a*r2^2 + 2*r2^3)), 2*(a^2 + sqrt(3)*a*r2 + r2^2)/(sqrt(3)*a + 2*r2))
   x1 = r1
   y1 = √3a - sqrt(R^2 - a^2)
   x0 = a
   y0 = -sqrt(R^2 - x0^2)
   @printf("正三角形の一辺の長さが %g，小円の直径が %g のとき，大円の直径は %g である。\n", 2a, 2r2, 2r1)
   @printf("a = %g;  r2 = %g;  R = %g;  r1 = %g;  x1 = %g;  y1 = %g;  x0 = %g;  y0 = %g\n", a, r2, R, r1, x1, y1, x0, y0)
   plot([a, 0, -a, a], [y0, √3a + y0, y0, y0], color=:green, lw=0.5)
   circle(0, 0, R, :magenta)
   circle2(x1, y1, r1)
   circle(0, R - r2, r2, :blue)
   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(0, √3a + y0, "√3a+y0  ", :green, :right, delta=-delta/2)
       point(x1, y1, "大円:r1,(x1,y1)", :red, :center, delta=-delta)
       point(0, R - r2, "小円:r2,(0,R-r2)", :blue, :center, :bottom, delta=2delta)
       point(x0, y0, " (x0,y0)", :green, :left, :vcenter)
   end
end;

draw(7.5/2, 4.755/2, true)

算額（その2115）

算額（その2115）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円15個，外円，菱形3個
#Julia, #SymPy, #算額, #和算

外円の中に菱形 3 個，等円 4 個を容れる。等円の直径が 1 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
菱形の頂点と円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     y::positive, x0::positive, y0::positive;
y0 = sqrt(R^2 - x0^2)
eq1 = x^2 + y^2 - (R - r)^2
eq2 = dist2(-R, 0, x0, y0, x, y, r)
eq3 = y /R - (R - y0)/y0
eq4 = 2y0 - y - R;

eq1, eq3, eq4 を連立させて解き，x, y, x0 を求める。

res = solve([eq1, eq3, eq4], (x, y, x0))[1]

   (sqrt(-(sqrt(2)*R - r)*(-2*R + sqrt(2)*R + r)), R*(-1 + sqrt(2)), sqrt(2)*R/2)

eq2 の x, y, x0 に代入し，simplify する。

eq12 = eq2(x => res[1], y => res[2], x0 => res[3]) |> simplify
eq12 |> println

   -2*R^2 + 3*sqrt(2)*R^2/2 - R*r + sqrt(2)*R*r/2 - r^2/2 - sqrt(2)*r^2/4

R について解く。

ans_R = solve(eq12, R)[1] |> factor |> println

   r*(4 + 3*sqrt(2))/2

外円の半径 R は，等円の半径 r の (4 + 3√2)/2 倍である。
等円の直径が 1 寸のとき，外円の直径は (4 + 3√2)/2 = 4.121320343559643 寸である。

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

draw(1/2, true)

算額（その2114）

算額（その2114）

百四十 群馬県甘楽郡妙義町下高田 高太神社 明治22年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円16個，菱形
#Julia, #SymPy, #算額, #和算

菱形の中に等円を 16 個容れる。菱形の対角線の長いほうが 80 寸，菱形の一辺の長さが 50 寸のとき，等円の直径はいかほどか。

本問は，算額（その928）の類題である。

菱形の対角線を 2a, 2b;  a > b
等円の半径と中心座標を r, (x, r), (r, y), ((2r + x)/3, (r + 2y)/3), ((r + 2x)/3, (2r + y)/3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     r::positive, x::positive, y::positive
b = sqrt(c^2 - a^2) 
eq1 = dist2(a, 0, 0, b, r, y, r)
eq2 = dist2(a, 0, 0, b, x, r, r)
eq3 = (x - r)^2 + (y - r)^2 - 36r^2;

eq1, eq2 から y, x を求める。

ans_y = solve(eq1, y)[1]
ans_y |> println
ans_x = solve(eq2, x)[1]
ans_x |> println

   (-c*r + (a - r)*sqrt(-a^2 + c^2))/a
   (a^2 - a*c + r*sqrt(-a^2 + c^2))/(a - c)

y, x を eq3 に代入する。

eq13 = eq3(x => ans_x, y => ans_y) |> expand |> simplify |> numerator
eq13 |> println

   a^3*c^2 - 36*a^3*r^2 - a^2*c^3 - 2*a^2*c^2*r + 36*a^2*c*r^2 + 2*a*c^3*r + 2*a*c^2*r*sqrt(-a^2 + c^2) - 2*c^3*r^2 - 2*c^2*r^2*sqrt(-a^2 + c^2)

r について解く。

ans_r = solve(eq13, r)[2]  # 2 of 2
ans_r |> println

   a*c*(6*a*(a - c) + c*(-a + c + sqrt(-a^2 + c^2)))/(2*(18*a^3 - 18*a^2*c + c^3 + c^2*sqrt(-a^2 + c^2)))

等円の直径は 2*4.54545454545454 = 9.09090909090908 = 100/11 寸である。

ans_r(a => 80/2, c=> 50, b => 30) |> println

   4.54545454545454

ans_x, ans_y は a, c の他に r を含む式であるが，上で得られた ans_r を代入すればよい。

ans_x(a => 80/2, c=> 50, r => 50/11) |> println  # 290/11
ans_y(a => 80/2, c=> 50, r => 50/11) |> println  # 230/11

   26.3636363636364
   20.9090909090909

function draw(a, c, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = sqrt(c^2 - a^2)
   r = a*c*(6*a*(a - c) + c*(-a + c + sqrt(-a^2 + c^2)))/(2*(18*a^3 - 18*a^2*c + c^3 + c^2*sqrt(-a^2 + c^2)))
   y = (-c*r + (a - r)*sqrt(-a^2 + c^2))/a
   x = (a^2 - a*c + r*sqrt(-a^2 + c^2))/(a - c)  
   @printf("a = %g;  b = %g;  c = %g\n", a, b, c)
   @printf("r = %g;  x = %g;  y = %g\n", r, x, y)
   @printf("等円の直径は %g 寸である。\n", 2r)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
   circle4(r, y, r)
   circle4((2r + x)/3, (r + 2y)/3, r)
   circle4((r + 2x)/3, (2r + y)/3, r)
   circle4(x, r, r)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(r, y, "(r,y)", :red, :right, delta=-delta, deltax=4delta)
       point(x, r, "(x,r)", :red, :right, delta=-delta, deltax=4delta)
       point(0, b, " b", :green, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
   end
end;

draw(80/2, 50, true)

算額（その2113）

算額（その2113）

百三十八 群馬県利根郡月夜野町上津 八幡神社 明治22年(1889)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円39個，円環
#Julia, #SymPy, #算額, #和算

大円の内外にそれぞれ 39 個の等円（外等円，内等円と呼ぶ）が互いに接している。
外等円と内等円の直径の差が 8.454 寸，大円と内等円の直径の和が 171.305 寸，「外余積が 109.688535160125 歩」のとき，大円，外等円，内等円の直径はいかほどか。

算額（その672）と同じであるが，問題の提示法に難がある。

大円の半径と中心座標を R, (0, 0)
外等円の半径と中心座標を r1, ((R + r1)*cos(pi/19), (R + r1)*sin(pi/19)); (R + r1)*sin(pi/19) = r1
内等円の半径と中心座標を r2, ((R - r2)*cos(pi/19), (R - r2)*sin(pi/19)); (R - r2)*sin(pi/19) = r2
とする。
未知数は R, r1, r2 の 3 個なので，方程式も 3 本用意すれば解ける。
外等円と内等円は互いに接し合っているという 2 条件は eq1, eq2 で表されている。
残り1つの条件が必要であるが，「問」には 3 個の条件が提示されている。3 番目の条件に出てくる「外余積」は有効数字が多いので，この条件を使えば正確な解が得られると思われるが，この「外余積」が何であるかわからない。
残る 2 条件のいずれを使うかで解が異なる。それは「外等円と内等円の直径の差が 8.454 寸」も「大円と内等円の直径の和が 171.305 寸」も有理数ではなく誤差を含むのが原因である。
更に，後で分かるが，eq1, eq2 に使われている「円周率」は 3.16 である（円積率=0.79）。

1.　「外等円と内等円の直径の差が 8.454 寸」を使う場合

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 =3.16
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
eq3 = 2(r1 - r2) - 8.454;

res = solve([eq1, eq2, eq3], (r1, r2, R))
res |> println

   Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(r2 => 10.6530275137650, R => 75.0022912276876, r1 => 14.8800275137650)

大円の直径 = 2*75.0022912276876 = 150.0045824553752
外等円の直径 = 2*14.8800275137650 = 29.76005502753
内等円の直径 = 2*10.6530275137650 = 21.30605502753

2.　「大円と内等円の直径の和が 171.305 寸」を使う場合

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 =3.16
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
eq3 = 2(R + r2) - 171.305;

res = solve([eq1, eq2, eq3], (r1, r2, R))
res |> println

   Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(r2 => 10.6526769444053, R => 74.9998230555947, r1 => 14.8795378424756)

大円の直径 = 2*74.9998230555947 = 149.9996461111894
外等円の直径 = 2*14.8795378424756 = 29.7590756849512
内等円の直径 = 2*10.6526769444053 = 21.3053538888106

3.　大円の直径を与える場合

r1, r2 は 大円 R の関数なので，条件は eq1, eq2 だけで十分である。求められた r1, r2 は R を含む式なので，R として任意の値を設定すれば，その場合の外等円，内等円の半径（直径）はすぐに得られる。

大円の直径 2R は，前 2 項から 150 寸と推定される（実際に「答」で 150 寸とされている）。

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 = Sym(316)/100
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
ans_r1 = solve(eq1, r1)[1]
ans_r2 = solve(eq2, r2)[1]
(ans_r1, ans_r2) |> println
(2ans_r1(R => 150/2).evalf(), 2ans_r2(R => 150/2).evalf())

   (R*sin(79/475)/(1 - sin(79/475)), R*sin(79/475)/(sin(79/475) + 1))

   (29.7591458944761, 21.3054041537713)

外等円，内等円の半径はそれぞれ
R*sin(79/475)/(1 - sin(79/475))
R*sin(79/475)/(sin(79/475) + 1)
であり，R = 150/2 を代入し 2 倍すると，それぞれの直径は 29.7591458944761, 21.3054041537713 となり，小数点以下 3 桁までとると，29.759, 21.305 となり，「答」と一致する。

4.　eq1, eq2 を，円周率 = π で計算する場合

近似計算ではなく正確な解を求める。

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 = PI
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
ans_r1 = solve(eq1, r1)[1]
ans_r2 = solve(eq2, r2)[1]
(ans_r1, ans_r2) |> println
(2ans_r1(R => 150/2).evalf(), 2ans_r2(R => 150/2).evalf())

   (R*sin(pi/19)/(1 - sin(pi/19)), R*sin(pi/19)/(sin(pi/19) + 1))

   (29.5535416156890, 21.1998138649765)

大円の直径が 150 寸のとき，外等円の直径は 29.5535 寸, 内等円の直径は 21.1998 寸である。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (R*sin(pi/19)/(1 - sin(pi/19)), R*sin(pi/19)/(sin(pi/19) + 1))
   @printf("大円の直径が %g のとき，外等円の直径は %g, 内等円の直径は %g である。\n", 2R, 2r1, 2r2)
   @printf("R = %g;  r1 = %g;  r2 = %g\n", R, r1, r2)
   plot()
   circle(0, 0, R, :magenta)
   circle(0, 0, R + r1, :gray90)
   circle(0, 0, R - r2, :gray90)
   rotate(0, -r1 - R, r1, :blue, angle=360/19) 
   rotate(0, r2 - R, r2, angle=360/19) 
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;

draw(150/2, true)

算額（その2112）

算額（その2112）

百三十七 群馬県藤岡市藤岡 金光寺 明治21年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，円弧，最大化
#Julia, #SymPy, #算額, #和算

団扇の中に円弧，甲円，乙円を容れる。団扇の直径，円弧の直径がそれぞれ 4 寸，3 寸のとき，乙円の直径が最大になるのは甲円の直径がいかほどのときか。

団扇の直径と中心座標を R, (0, 0)
円弧の直径と中心座標を r0, (0, y0); y0 < 0
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r0::positive, y0::positive,
     r1::positive, r2::positive, x2::positive, y2::positive;
y0 = 2r1 - R - r0
eq1 = x2^2 + (y2 - y0)^2 - (r0 - r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x2^2 + (y2 - r1 + R)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r2, x2, y2))[2]  # 2 of 2

   (-r1*(-R + r1)*(r0 - r1)/(R*r0 - r1^2), 2*sqrt(R)*sqrt(r0)*r1*sqrt(-(-R + r1)*(r0 - r1))/(R*r0 - r1^2), (-R^2*r0 + R*r0*r1 + r0*r1^2 - r1^3)/(R*r0 - r1^2))

ans_r2 = res[1];

乙円の半径 r2 は -r1*(-R + r1)*(r0 - r1)/(R*r0 - r1^2) である。
R, r0 が与えられるので，変化しうるのは甲円の半径 r1 である。
r2 と r1 の関係を下図に示す。r1 が 0.5〜1 の間で r2 が最大値をとることがわかる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(ans_r2(R => 4/2, r0 => 3/2), xlims=(0,1.5), xlabel="r1", ylabel="r2")
vline!([1.37879462408182/2])
hline!([0])

r2 が最大になるときの r1 は，r2 の式を r1 で偏微分し，その式が 0 になるときの r1 を求めればよい。

r10 = solve(diff(ans_r2, r1), r1)[2];  # 2 of 4
r10 |> println
2r10(R => 4/2, r0 => 3/2) |> println

   -sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2 + sqrt(-2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 16*R*r0/3 + (4*R^2*r0 + 4*R*r0^2)/sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) - 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2

   1.37879462408182

乙円の直径が最大になるときの甲円の直径を表す式はとても長く，SymPy では簡約化できないが，数値は求めることができる。

団扇の直径が 4 寸, 円弧の直径が 3 寸で，甲円の直径が 1.37879 寸のときに乙円の直径は最大値 0.580181 寸になる。

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

    R = 2;  r0 = 1.5;  y0 = -2.12121;  r1 = 0.689397;  r2 = 0.290091;  x2 = 0.974955;  y2 = -1.40473

function draw(R, r0, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = -sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2 + sqrt(-2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 16*R*r0/3 + (4*R^2*r0 + 4*R*r0^2)/sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) - 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2
   (r2, x2, y2) = (-r1*(-R + r1)*(r0 - r1)/(R*r0 - r1^2), 2*sqrt(R)*sqrt(r0)*r1*sqrt(-(-R + r1)*(r0 - r1))/(R*r0 - r1^2), (-R^2*r0 + R*r0*r1 + r0*r1^2 - r1^3)/(R*r0 - r1^2))
   y0 = 2r1 - R - r0
   @printf("団扇の直径が %g, 円弧の直径が %g で，甲円の直径が %g のときに乙円の直径は最大値 %g になる。\n", 2R, 2r0, 2r1, 2r2)
   @printf("R = %g;  r0 = %g;  y0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", R, r0, y0, r1, r2, x2, y2)
   plot()
   circle(0, 0, R)
   circle(0, y0, r0, :blue)
   circle(0, r1 - R, r1, :green)
   circle2(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(0, 0, "団扇:R,(0,0)", :red, :center, delta=-delta/2)
       point(0, y0, "円弧:r0,(0,y0)", :blue, :center, delta=-delta/2)
       point(0, r1 - R, "甲円:r1\n(0,r1-R)", :green, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :magenta, :left, :bottom, delta=4delta, deltax=-4delta)
   end
end;

draw(4/2, 3/2, true)

算額（その2111）

算額（その2111）

百三十七 群馬県藤岡市藤岡 金光寺 明治21年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，楕円，最大化
#Julia, #SymPy, #算額, #和算

大円と楕円が交わり，其の隙間に等円 3 個が入っている。楕円の長径が与えられたとき，等円の直径が最大になるのは楕円の短径がいかほどのときか。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
大円の半径と中心座標を r1, (0, r1 - b)
等円の半径と中心座標を r2, (x2, 0), (0, b + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive,
     r2::positive, x2::positive;
r1 = b + r2
eq1 = x2^2 + (r1 -b)^2 - (r1 + r2)^2
eq2 = x2^2 - (a^2 - b^2)*(b^2 - r2^2)/b^2
res = solve([eq1, eq2], (r2, x2))[1]

   ((8*a^2*b^2*(a - b)*(a + b)/(a^2 + 2*b^2) - 4*a^2*b^2 + 4*b^4)/(4*b*(a - b)*(a + b)), 2*sqrt(2)*a*b*sqrt(a - b)*sqrt(a + b)/(a^2 + 2*b^2))

x2 は簡約化できる。

ans_x2 = res[2] |> simplify
ans_x2 |> println

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

等円の半径も簡約化され，a, b を含む関数である。a は与えられるので，変化できるのは b である。 0 < b < √2a/2

ans_r2 = res[1] |> simplify
ans_r2 |> println

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

a = 1 のとき，楕円の短半径 b と小円の半径 r2 は以下の図のようになる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(ans_r2(a => 1), xlims=(0,√2/2), xlabel="b", ylabel="r2")
vline!([0.34356074972251255])
hline!([0])

r2 が最大になるときの b は，r2 の式を b で偏微分し，その式が 0 になるときの b を求めればよい。

b0 = solve(diff(ans_r2, b), b)[1]
b0 |> println
b0(a => 1).evalf() |> println

   a*sqrt(-1 + sqrt(5)/2)
   0.343560749722512

楕円の長半径 a のとき，等円の半径が最大になるのは，楕円の短半径が b = a*sqrt(√5/2 - 1) のときである。

a = 1 のとき，b = sqrt(√5/2 - 1) = 0.343560749722512 である。

そのときの等円の半径は最大値 sqrt(√5 - 2)*(√10 - √2)/4 = 0.212332220528909 になる。

r2_max = ans_r2(b => b0) |> simplify |> factor
r2_max |> println
r2_max(a => 1).evalf() |> println

   a*sqrt(-2 + sqrt(5))*(-sqrt(2) + sqrt(10))/4
   0.212332220528909

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = a*sqrt(-1 + sqrt(5)/2)
   r2 = a*sqrt(√5 - 2)*(√10 - √2)/4
   r1 = b + r2
   x2 = 2a*b*sqrt(2a^2 - 2b^2)/(a^2 + 2b^2)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle(0, r1 - b, r1, :green)
   circle2(x2, 0, r2, :blue)
   circle(0, b + r2, r2, :blue)
   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", :red, :left, :bottom, delta=delta/2, deltax=delta)
       point(0, b, "b", :red, :center, :bottom, delta=2delta)
       point(0, r1 - b, "大円:r1,(0,r1-b)", :green, :center, delta=-2delta)
       point(0, b + r2, "等円:r2\n(0,b+r2)", :blue, :center, :bottom, delta=2delta)
       point(x2, 0, "等円:r2\n(x2,0)", :blue, :center, :bottom, delta=2delta)
   end
end;

draw(30, true)

算額（その2110）

算額（その2110）

百三十七 群馬県藤岡市藤岡 金光寺 明治21年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，扇
#Julia, #SymPy, #算額, #和算

扇面に大円 1 個，等円 3 個を容れる。扇長（要から先端までの長さ）が与えられたとき，大円の直径を求める術を述べよ。

扇長を R，扇の端（図参照）の座標を (x0, sqrt(R^2 - x0^2))
大円の半径と中心座標を r1, (0, R - r1)
等円の半径と中心座標を r2, (0, R - 2r1 - r2), (x2, y2)

残念ながら，答えを得るための「術」は SymPy の能力では得られないようだ（何らかの手立てはあるはず）。
以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms R::positive, x0::positive, r1::positive,
     r2::positive, x2::positive, y2::positive;
sinθ = x0/R
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x2^2 + (y2 - R + r1)^2 - (r1 + r2)^2
eq3 = dist2(0, 0, x0, sqrt(R^2 - x0^2), 0, R - r1, r1)
eq4 = dist2(0, 0, x0, sqrt(R^2 - x0^2), 0, R - 2r1 - r2, r2)
eq5 = dist2(0, 0, x0, sqrt(R^2 - x0^2), x2, y2, r2);

using NLsolve

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

function H(u)
   (r1, r2, x2, y2, x0) = u
   return [
       x2^2 + y2^2 - (R - r2)^2,  # eq1
       x2^2 - (r1 + r2)^2 + (-R + r1 + y2)^2,  # eq2
       -R^2*r1^2 + R^2*x0^2 - 2*R*r1*x0^2 + r1^2*x0^2,  # eq3
       -R^2*r2^2 + R^2*x0^2 - 4*R*r1*x0^2 - 2*R*r2*x0^2 + 4*r1^2*x0^2 + 4*r1*r2*x0^2 + r2^2*x0^2,  # eq4
       (R^2*(-r2^2 + x2^2) - x0^2*x2^2 + x0^2*y2^2 - 2*x0*x2*y2*sqrt(R^2 - x0^2))/R^2,  # eq5
   ]
end;

R = 30
iniv = BigFloat[11, 3, 13, 24, 17]
res = nls(H, ini=iniv)

   ([10.98076211353316, 2.942286340599478, 13.126775968941422, 23.660254037844386, 17.320508075688775], true)

たとえば，扇長が 30 寸のとき，大円の直径は 2*10.98076211353316 = 21.96152422706632 寸である。

「術」は「3 の平方根から 1 を引き，扇長を掛ける」と，実に簡明な答えが示されている。
上の数値解は正しく，正確なものであることがわかる。

2*10.98076211353316 |> println
(√3 - 1)*30 |> println

   21.96152422706632
   21.961524227066317

function draw(R, more)
    pyplot(size=(700, 700), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x2, y2, x0) = res[1]
   y0 = sqrt(R^2 - x0^2)
   θ = atand(y0, x0)
   plot([-x0, 0, x0], [y0, 0, y0], color=:red, lw=0.5)
   circle(0, R - r1, r1, :green)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ)
   circle2(x2, y2, r2, :blue)
   circle(0, R - 2r1 - r2, r2, :blue)
   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 - r1, "大円:r1,(0,R-r1)", :green, :center, delta=-delta/2)
       point(0, R - 2r1 - r2, "等円:r2\n(0,R-2r1-r2)", :black, :center, :bottom, delta=delta/2)
       point(x2, y2, "等円:r2\n(x2,y2)", :black, :center, :bottom, delta=delta/2)
       point(0, R, "扇長:R", :red, :center, :bottom, delta=delta/2)
       point(x0, sqrt(R^2 - x0^2), " (x0,√(R^2-x0^2))", :red, :left, :vcenter)
       xlims!(-x0 - 2delta, x0 + 30delta)
   end
end;

draw(30, true)

算額（その2109）

算額（その2109）

百三十四 群馬県富岡市一ノ宮 貫前神社 明治20年(1887)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円8個
#Julia, #SymPy, #算額, #和算

大円 3 個と小円 1 個が交わり，その隙間に等円 3 個，甲円 1 個を容れる。小円の直径が 5 寸，等円の直径が 2 寸のとき，甲円の直径はいかほどか。

大円の半径と中心座標を r1, (0, 0), (r3 + r1)
小円の半径と中心座標を r2, (0, r3 + r2)
甲円の半径と中心座標を r3, (0, 0)
等円の半径と中心座標を r4, (r3 + r4, 0), (0, r3 + r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive,
     r3::positive, r4::positive;
eq1 = r3 + 2r4 - r1
eq2 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, r3))[1]

   (sqrt(r4)*sqrt(4*r2 + r4)/2 + 3*r4/2, sqrt(r4)*sqrt(4*r2 + r4)/2 - r4/2)

大円の半径は (√r4*sqrt(4r2 + r4) + 3r4)/2
甲円の半径は (√r4*sqrt(4r2 + r4) - r4)/2
である。

小円の直径が 5 寸，等円の直径が 2 寸のとき，大円の直径は 6.31662479035540 寸，甲円の直径は 2.31662479035540 寸である。

2res[1](r2 => 5/2, r4 => 2/2) |> println
2res[2](r2 => 5/2, r4 => 2/2) |> println

   6.31662479035540
   2.31662479035540

function draw(r2, r4, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (sqrt(r4)*sqrt(4*r2 + r4)/2 + 3*r4/2, sqrt(r4)*sqrt(4*r2 + r4)/2 - r4/2)
   r1 = (√r4*sqrt(4r2 + r4) + 3r4)/2
   r3 = (√r4*sqrt(4r2 + r4) - r4)/2
   @printf("小円の直径が %g，等円の直径が %g のとき，大円の直径は %g，甲円の直径は %g である。", 2r2, 2r4, 2r1, 2r3)
   plot()
   circle(0, 0, r1, :blue)
   circle2(r3 + r1, 0, r1, :blue)
   circle(0, r3 + r2, r2, :magenta)
   circle(0, 0, r3, :green)
   circle2(r3 + r4, 0, r4)
   circle(0, r3 + r4, r4)
   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(r3 + r1, 0, "大円:r1\n(r3+r1,0)", :blue, :center, delta=-delta/2)
       point(0, r3 + r2, "小円:r2\n(0,r3+r2)", :magenta, :center, :bottom, delta=delta/2)
       point(r3, 0, "r3 ", :green, :right, :vcenter)
       point(0, 0, "甲円", :green, :center, :vcenter, mark=false)
       point(r3 + r4, 0, "r3+r4", :red, :center, delta=-delta/2)
       point(r3 + r4, 0, "等円", :red, :center, :bottom, delta=delta/2, mark=false)
   end
end;

draw(5/2, 2/2, true)

 

算額（その2108）

算額（その2108）

百三十四 群馬県富岡市一ノ宮 貫前神社 明治20年(1887)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，直角三角形
#Julia, #SymPy, #算額, #和算

直角三角形に内接する大円は等円の直径の計算には無関係なお飾りである。大円を除けば，算額（その928）の第一象限の図と同じになる。

直角三角形の 3 辺を，「鈎」，「股」，「弦」
大円の半径と中心座標を r0, (r0, r0)
等円の半径と中心座標を r1, (r1, y1), (x1, r1), (x, y); x = (r1 + x1)/2, y = (y1 + r1)/2
とおき，以下の連立方程式を解く。

1. 大円の半径 r0 は以下の式で求められる。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive,
     x1::positive, y1::positive,
     x::positive, y::positive,
     鈎::positive, 股::positive, 弦::positive;
eq0 = 鈎 + 股 - sqrt(鈎^2 + 股^2) ⩵ 2r0
ans_r0 = solve(eq0, r0)[1]
ans_r0 |> println
ans_r0(鈎 => 5.4, 股 => 7.2) |> println

   股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
   1.80000000000000

2. r1,x1,y1 は SymPy の能力的に，同時には解けないので，逐次解いていく。

eq1 = dist2(0, 鈎, 股, 0, r1, y1, r1)
eq2 = dist2(0, 鈎, 股, 0, x1, r1, r1)
eq3 = (x1 - r1)^2 + (y1 - r1)^2 - (4r1)^2;

まず，eq1, eq2 から x1, y1 を求める。

(ans_x1, ans_y1) = solve([eq1, eq2], (x1, y1))[1]  # 1 of 4

   (-r1*sqrt(股^2 + 鈎^2)/鈎 - 股*(r1 - 鈎)/鈎, -r1*sqrt(股^2 + 鈎^2)/股 - 鈎*(r1 - 股)/股)

eq3 の x1, y1 に代入し，r1 を求める。

eq13 = eq3(x1 => ans_x1, y1 => ans_y1)
ans_r1 = solve(eq13, r1)[2]  # 2 of 2
ans_r1 |> println

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

鈎，股が 5.4 寸, 7.2 寸のとき，等円の直径は 2 寸である。

ans_r1(鈎 => 5.4, 股 => 7.2) * 2 |> println

   2.00000000000000

sqrt(股^2 + 鈎^2) が何回も出てくるので，「弦」で置き換える。

ans_r1 = ans_r1(sqrt(股^2 + 鈎^2) => 弦) |> simplify
ans_r1 |> println

   股*鈎*(弦*股^2 - 4*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/(2*(弦*股^3 + 弦*股^2*鈎 + 弦*股*鈎^2 + 弦*鈎^3 + 股^4 + 股^3*鈎 - 6*股^2*鈎^2 + 股*鈎^3 + 鈎^4))

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

   鈎 = 5.4;  股 = 7.2;  弦 = 9;  r0 = 1.8;  r1 = 1;  x1 = 4.2;  y1 = 3.4;  x = 2.6;  y = 2.2

function draw(鈎, 股, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   弦 = sqrt(鈎^2 + 股^2)
   r0 = 股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
   r1 = 股*鈎*(弦*股^2 - 4*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/(2*(弦*股^3 + 弦*股^2*鈎 + 弦*股*鈎^2 + 弦*鈎^3 + 股^4 + 股^3*鈎 - 6*股^2*鈎^2 + 股*鈎^3 + 鈎^4))
   (x1, y1) = (-r1*sqrt(股^2 + 鈎^2)/鈎 - 股*(r1 - 鈎)/鈎, -r1*sqrt(股^2 + 鈎^2)/股 - 鈎*(r1 - 股)/股)
   (x, y) = ((r1 + x1)/2, (y1 + r1)/2)
   @printf("鈎，股が %g, %g のとき，等円の直径は %g である。\n", 鈎, 股, 2r1)
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  r0 = %g;  r1 = %g;  x1 = %g;  y1 = %g;  x = %g;  y = %g\n", 鈎, 股, 弦, r0, r1, x1, y1, x, y)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   circle(r0, r0, r0, :blue)
   circle(r1, y1, r1)
   circle(x, y, r1)
   circle(x1, r1, r1)
   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, 鈎, "鈎", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
       point(股, 0, "股", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
       point(r0, r0, "大円:r0,(r0,r0)", :blue, :center, delta=-delta)
       point(r1, y1, "等円:r1,(r1,y1)", :red, :center, delta=-delta)
       point(x, y, "等円:r1,(x,r)", :red, :center, delta=-delta)
       point(x1, r1, "等円:r1,(x1,r1)", :red, :center, delta=-delta)
   end
end;

draw(5.4, 7.2, true)

#SymPy, #算額, #和算

算額（その2107）

算額（その2107）

百三十四 群馬県富岡市一ノ宮 貫前神社 明治20年(1887)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個，正三角形

正三角形の内外に甲円，乙円を容れ，両円を包括する円（無名）の上に正三角形の頂点を通る天円を置く。また，無名円と正三角形の間に丙円を置く。天円の直径が 6 寸，乙円の直径が 1 寸のとき，丙円の直径はいかほどか。

正三角形の一辺の長さを 2a
無名円の半径と中心座標を r0, (0, √3a - 2r1 - r0)
天円の半径と中心座標を r1, (0, √3a - r1)
甲円の半径と中心座標を r2, (0, r2)
乙円の半径と中心座標を r3, (0, -r3)
丙円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
SymPy の能力上，一度に解くことができないので，逐次解いていく。

include("julia-source.txt");

using SymPy

@syms a::positive, r0::positive,
     r1::positive, r2::positive,
     r3::positive, r4::positive, x4::positive;
r0 = r2 + r3
eq1 = 2r1 + 2r2 - √Sym(3)a
eq2 = x4^2 + (√Sym(3)a - 2r1 - r0 - r4)^2 - (r0 + r4)^2
eq3 = dist2(a, 0, 0, √Sym(3)a, 0, √Sym(3)a - 2r1 - r0, r0)
eq4 = dist2(a, 0, 0, √Sym(3)a, x4, r4, r4);

ans_r2 = solve(eq3, r2)[1]
ans_r2 |> println

   2*r1 - r3

eq11 = eq1(r2 => ans_r2)
ans_a = solve(eq11, a)[1]
ans_a |> println

   2*sqrt(3)*(3*r1 - r3)/3

eq12 = eq2(r2 => ans_r2, a => ans_a) |> simplify
ans_x4 = solve(eq12, x4)[1]
ans_x4 |> println

   2*sqrt(2*r1 - r3)*sqrt(r3 + r4)

eq14 = eq4(r2 => ans_r2, a => ans_a, x4 => ans_x4) |> simplify
ans_r4 = solve(eq14, r4)[2]  # 2 of 3
ans_r4 |> println

   -4*sqrt(2)*sqrt(r1)*sqrt(2*r1 - r3)/3 + 10*r1/3 - 4*r3/3

eq1, eq2, eq3 を連立させて a, r2, x4 を同時に解くことはできる。前述の最終段階と同じく，得られた解を eq4 に代入して解けば r4 が得られる。

(ans_a, ans_r2, ans_x4) = solve([eq1, eq2, eq3], (a, r2, x4))[1]

   (2*sqrt(3)*(3*r1 - r3)/3, 2*r1 - r3, sqrt(8*r1*r3 + 8*r1*r4 - 4*r3^2 - 4*r3*r4))

ans_x4 は以下のように簡約化できる。

ans_x4 |> factor |> println

   2*sqrt(2*r1 - r3)*sqrt(r3 + r4)

最後に r4 を求める。

eq14 = eq4(r2 => ans_r2, a => ans_a, x4 => ans_x4) |> simplify
ans_r4 = solve(eq14, r4)[2]  # 2 of 3
ans_r4 |> println

   -4*sqrt(2)*sqrt(r1)*sqrt(2*r1 - r3)/3 + 10*r1/3 - 4*r3/3

いずれにせよ，天円の直径が 6 寸，乙円の直径が 1 寸のとき，丙円の直径は 3.34783294256526 寸である。

ans_r4(r1 => 6/2, r3 => 1/2).evalf() * 2 |> println

   3.34783294256526

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

   r1 = 3;  r3 = 0.5;  r2 = 5.5;  a = 9.81495;  r4 = 1.67392;  x4 = 6.91565;  r0 = 6

「答」では「丙円径三寸」となっている？
天円の直径が 6，乙円の直径が 1.55 のとき，丙円の直径は 3.00238 になる。

function draw(r1, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2r1 - r3
   a = 2√3(3r1 - r3)/3
   r4 = 10r1/3 - 4r3/3 - 4√2*sqrt(r1)*sqrt(2r1 - r3)/3
   x4 = 2sqrt(2r1 - r3)*sqrt(r3 + r4)
   r0 = r2 + r3
   @printf("天円の直径が %g，乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r3, 2r4)
   @printf("r1 = %g;  r3 = %g;  r2 = %g;  a = %g;  r4 = %g;  x4 = %g;  r0 = %g\n", r1, r3, r2, a, r4, x4, r0)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:green, lw=0.5)
   circle(0, √3a - r1, r1)
   circle(0, r2, r2, :blue)
   circle(0, √3a - 2r1 - r0, r0, :magenta)
   circle(0, - r3, r3, :orange)
   circle2(x4, r4, r4, :brown)
   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", :green, :left, :bottom, delta=delta/2, deltax=delta/2) 
       point(0, √3a, "√3a", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
       point(0, √3a - r1, "天円:r1\n(0,√3a-r1)", :red, :center, delta=-delta)
       point(0, r2, "甲円:r2,(0,r2)", :blue, :center,:bottom, delta=delta)
       point(0, -r3, "乙円:r3,(0,-r3)", :black, :left, :vcenter, deltax=2delta)
       point(x4, r4, "丙円:r4\n(x4,r4)", :brown, :center, deltax=-delta)
       point(0, √3a - 2r1 - r0, "無名円:r0,(0,√3a-2r1-r0)", :magenta, :center, delta=-delta)
   end
end;

draw(6/2, 1/2, true)

#SymPy #算額 #和算

算額（その2106）

算額（その2106）

百三十四 群馬県富岡市一ノ宮 貫前神社 明治20年(1887)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，直角三角形，中鈎

直角三角形の中に中鈎（直角の頂点から斜辺への垂線）を引き，区画された領域に大円 1 個，等円 3 個を容れる。鈎が 5.4 寸，股が 7.2 寸，大円の直径が 2.88 寸のとき，等円の直径はいかほどか。

鈎，股を「鈎」，「股」
中鈎と斜辺（弦）の交点座標を (x, y)
大円の半径と中心座標を r1, (x1, r1)
等円の半径と中心座標を r2, (r2, y2), (r2, y2 + 2r2), (r2, y2 + 4r2)
とおき，以下の連立方程式を解く。
SymPy の能力的に，一度に解けないので，逐次解いていく。

まず，中鈎と斜辺の交点座標 (x, y) を求める。
eq5, eq6 を解いて x = 股*鈎^2/(股^2 + 鈎^2), y = 股^2*鈎/(股^2 + 鈎^2) を得る。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive,
     x::positive, y::positive
eq5 = y*股 - 鈎*(股 - x);
eq6 = x/y - 鈎/股;
res = solve([eq5, eq6], (x, y))
res[x] |> println
res[y] |> println

   股*鈎^2/(股^2 + 鈎^2)
   股^2*鈎/(股^2 + 鈎^2)

@syms 鈎::positive, 股::positive,
     x::positive, y::positive,
     r1::positive, x1::positive,
     r2::positive, y2::positive
x = 股*鈎^2/(股^2 + 鈎^2)
y = 股^2*鈎/(股^2 + 鈎^2)

eq1 = dist2(0, 0, x, y, x1, r1, r1)
eq2 = dist2(0, 0, x, y, r2, y2, r2)
eq4 = dist2(股, 0, 0, 鈎, r2, y2 + 4r2, r2);

この問題では，「鈎が 5.4 寸，股が 7.2 寸，大円の直径が 2.88 寸」というように，条件が 3 つ与えられている。しかし，鈎と股が決まれば中鈎が決まり，大円が入る直角三角形が決まれば大円の直径は自ずと決まる。つまり，大円の直径の情報は余分である上にもし，個の数値を誤って与えてしまうと図形が成立しない。

鈎，股の 2 条件から，内接する大円の直径をは以下のようにして求めることができる。
中鈎の長さを「中鈎」，中鈎の脚で分割される斜辺（弦）の長い方を「長弦」とする。

@syms 中鈎::positive, 長弦::positive
中鈎 = sqrt(x^2 + y^2);

長弦を求める。

eq01 = sqrt(中鈎^2 + 長弦^2) - 股
ans_長弦 = solve(eq01, 長弦)[1]

直角三角形に内接する円の直径は「鈎 + 股 - 弦 = 直径」の関係があるので，以下の方程式を解く。鈎，股を与えれば半径が求まる。

eq0 = 中鈎 + ans_長弦 - 股 - 2r1
ans_r1 = solve(eq0, r1)[1](鈎 => 5.4, 股 => 7.2)
ans_r1 |> println

   1.44000000000000

さらに，この算額問題を解くにあたり，大円は必要ない。
しかし，図を描くために大円の中心座標（x 座標）を求める。

ans_x1 = solve(eq1, x1)[2]  # 2 of 2
ans_x1 |> println

   r1*(鈎 + sqrt(股^2 + 鈎^2))/股

さて，本来の目的である，等円の直径（と中心座標）を求める。
まず，eq2 に，最初に求めた鈎一番下の等円の中心座標（y 座標）を求める。

ans_y2 = solve(eq2, y2)[2]
ans_y2 |> println

   r2*(股 + sqrt(股^2 + 鈎^2))/鈎

最後に，eq4 の y2 に上で求めた ans_y2 を代入して，r2 を求める。

eq14 = eq4(y2 => ans_y2)
ans_r2 = solve(eq14, r2)[1]  # 1 of 2
ans_r2 |> println

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

式は長いが，以下のように鈎が 5.4 寸, 股が 7.2 寸を代入すれば，等円の直径は 1.2 寸と求まる。

ans_r2(鈎 => 5.4, 股 => 7.2) * 2 |> println

   1.20000000000000

以上をまとめると，鈎，股が与えられれば，以下の順に計算を進める。
   x = 股*鈎^2/(股^2 + 鈎^2)
   y = 股^2*鈎/(股^2 + 鈎^2)
   x1 = r1*(鈎 + sqrt(股^2 + 鈎^2))/股
   r2 = 鈎^2*(股^2 + 4*股*鈎 + 股*sqrt(股^2 + 鈎^2) + 鈎^2 - 鈎*sqrt(股^2 + 鈎^2))/(2*(股^3 + 4*股^2*鈎 + 股^2*sqrt(股^2 + 鈎^2) + 9*股*鈎^2 + 4*股*鈎*sqrt(股^2 + 鈎^2) + 4*鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2)))
   y2 = r2*(股 + sqrt(股^2 + 鈎^2))/鈎

function draw(鈎, 股, r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #(x1, r2, y2, x, y) = res[1]
   x = 股*鈎^2/(股^2 + 鈎^2)
   y = 股^2*鈎/(股^2 + 鈎^2)
   x1 = r1*(鈎 + sqrt(股^2 + 鈎^2))/股
   r2 = 鈎^2*(股^2 + 4*股*鈎 + 股*sqrt(股^2 + 鈎^2) + 鈎^2 - 鈎*sqrt(股^2 + 鈎^2))/(2*(股^3 + 4*股^2*鈎 + 股^2*sqrt(股^2 + 鈎^2) + 9*股*鈎^2 + 4*股*鈎*sqrt(股^2 + 鈎^2) + 4*鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2)))
   y2 = r2*(股 + sqrt(股^2 + 鈎^2))/鈎
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   circle(x1, r1, r1)
   circle(r2, y2, r2, :blue)
   circle(r2, y2 + 2r2, r2, :blue)
   circle(r2, y2 + 4r2, r2, :blue)
   segment(0, 0, x, y, :magenta)
   segment(x, y, 股, 0, :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(0, 鈎, " 鈎", :green, :left, :bottom, delta=delta/2)
       point(股, 0, " 股", :green, :left, :bottom, delta=delta/2)
       point(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(r2, y2, " 等円:r2,(r2,y2)", :blue, :left, :vcenter)
       point(r2, y2 + 2r2, " 等円:r2,(r2,y2+2r2)", :blue, :left, :vcenter)
       point(r2, y2 + 4r2, " 等円:r2,(r2,y2+4r2)", :blue, :left, :vcenter)
       point(x, y, "(x,y)", :magenta, :left, :bottom, delta=delta/2)
       point(x/4, y/5, "中鈎", :magenta, :left, :vcenter, mark=false)
       point((x + 股)/2, y/2, "  長弦", :black, :left, :vcenter, mark=false)
       xlims!(-2delta, 股 + 5delta)
   end
end;

draw(5.4, 7.2, 2.88/2, true)

算額（その2105）

算額（その2105）

百十六 群馬県吾妻郡吾妻町金井 一宮神社 明治5年(1872)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円1個，四分円，斜線

四分円と，四分円の端を通る傾きが -1 の斜線と，直線に囲まれた領域に円を容れる。四分円の半径が 0.5 寸のとき，円の直径はいかほどか。

注：斜線と四分円の交点の y 座標を「高」と呼んでいる。四分円の半径と同じ長さである。また，斜線を単に「斜」と呼んでいるが，算額の図からもわかるように，傾き -1 の斜線である。

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

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive
eq1 = dist2(0, r1, r1, 0, x2, r2, r2)
eq2 = (x2 + r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r2, x2))[3]  # 3 of 3

   ((r1*(3 - 2*sqrt(2)) + r1)^2/(4*r1), r1*(3 - 2*sqrt(2)))

res[1] |> simplify |> x -> x/r1 |> expand |> x -> r1*x |> println

   r1*(6 - 4*sqrt(2))

円の半径は四分円の半径の (6 - 4√2) 倍である。
四分円の半径が 0.5 寸のとき，円の直径は 2*0.5*(6 - 4√2) = 0.3431457505076194 寸である。

2 * 0.5*(6 - 4√2)

   0.3431457505076194

function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2) = ((r1*(3 - 2*sqrt(2)) + r1)^2/(4*r1), r1*(3 - 2*sqrt(2)))
   plot()
   circle(-r1, r1, r1, beginangle=270, endangle=360)
   segment(0, r1, r1, 0, :green)
   circle(x2, r2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r1, "r1", :red, :center, :bottom, delta=delta)
       point(-r1, 0, "-r1", :red, :center, delta=-delta)
       point(r1, 0, "r1", :green, :center, delta=-delta)
       point(-r1, r1, "四分円の中心", :red, :left, delta=-delta)
       point(x2, r2, "r2,(x2,r2)", :blue, :center, delta=-delta)
       ylims!(-6delta, r1 + 4delta)
   end
end;

draw(5, true)

算額（その2104）

算額（その2104）

百十五 群馬県富岡市一宮 貫前神社 明治4年(1871)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円1個，正三角形，正五角形

正三角形の中に円，長方形，五角形を容れる。円の直径が 10 寸のとき，正五角形の一辺の長さはいかほどか。

1. 円の半径を r とおく。IE = r
2. 円を内包する正三角形の一辺の長さを 2a とおく。EF = a = √3r
3. 正五角形を内包する円の半径を R とおく。AE = AB = AD = R
4. EF = AB*cosd(18°) = R*cosd(18°) = a = √3r
5. R*cosd(18°) = √3r より，R = √3r/cosd(18°)
6. 正五角形の一辺の長さ = 2CD = 2*R*sind(36°)
7. 2*R*sind(36°) = 2 * (√3r/cosd(18°)) * sind(36°) 
8 r = 10/2 のとき s = 2 * (√3(10/2)/cosd(18)) * sind(36) = 10.704662693192699
9. 図を描くために必要な正三角形の一辺の長さを 2b とすると，相似関係から CH = b = a + (yo + R)/√3
以上により，円の直径が 10 寸のとき，正五角形の一辺の長さは 10.704662693192699 寸である。

2 * (√3(10/2)/cosd(18)) * sind(36) 

   10.704662693192699

s は三角関数を含むが，cosd(18), sind(36) はルートを含む式で表せるので，以下のように簡略化される。
s = 5(√15 - √3) = 10.704662693192699

include("julia-source.txt");

using SymPy
@syms d, r
s = 2r * (√Sym(3)/cosd(Sym(18))) * sind(Sym(36)) |> simplify
s/r |> x -> apart(x, d) |> simplify |> x -> x*r |> println

   r*(-sqrt(3) + sqrt(15))

円の直径が 10 寸のとき，正五角形の一辺の長さは 5(√15 - √3) = 10.704662693192699。

function draw(r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = √3r
   R = a/cosd(18)
   xo = R*sind(36)
   yo = R*cosd(36)
   plot()
   circle(0, yo + R + r, r, :blue)
   plot!([a, 0, -a, a], (yo + R) .+ [0, √3a, 0, 0], color=:green, lw=0.5)
   θ = 18:72:378
   x = R.*cosd.(θ)
   y = R.*sind.(θ)
   circle(0, yo, R, :pink)
   plot!(x, yo .+ y, color=:red, lw=0.5)
   plot!([a, a, -a, -a, a], [0, yo + R, yo + R, 0, 0], color=:orange, lw=0.5)
   b = a + (yo + R)/√3
   plot!([b, 0, -b, b], [0, yo + R + √3a, 0, 0], color=:black, 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(0, yo, " A", :green, :left, :vcenter)
       point(x[1], yo + y[1], " B", :green, :left, :vcenter)
       point(0, 0, "C", :green, :center, delta=-delta)
       point(R*sind(36), 0, "D", :green, :center, delta=-delta)
       point(0, yo + R, "E", :green, :center, delta=-delta)
       point(R*cosd(18), yo + R, " F", :green, :left, :vcenter)
       point(a, 0, "G", :green, :left, delta=-delta)
       point(b, 0, "H", :green, :left, delta=-delta)
       point(0, yo + R + r, " I", :green, :left, :vcenter)
       point(0, yo + R + √3a, " J", :green, :left, :vcenter)
       point(0, yo + R + r, "")
   end
end;

draw(10/2, true)

算額（その2103）

算額（その2103）

百十五 群馬県富岡市一宮 貫前神社 明治4年(1871)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，外円，斜線，弦2本

外円の中に等斜を 2 本引き，全円を容れる。外円の直径が 10 寸，斜の長さが 8 寸のとき，全円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
全円の半径と中心座標を r, (0, R - r)
斜と外円の交点座標を (0, -R), (x, sqrt(R^2 - x^2)
斜の長さを 「斜」
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     斜::positive;
eq1 = (x^2 + (R - sqrt(R^2 - x^2))^2) - 斜^2
eq2 = dist2(0, -R, x, sqrt(R^2 - x^2), 0, R - r, r)
res = solve([eq1, eq2], (r, x))[1]  # 1 of 2

   (2*R*斜*(2*R*sqrt(128*R^8 - 128*R^6*斜^2 + 64*R^6*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) + 48*R^4*斜^4 - 32*R^4*斜^2*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) - 6*R^2*斜^6 + 8*R^2*斜^4*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) - 斜^6*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)) - 斜*(2*R - 斜)*(2*R + 斜)*(2*R^2 + sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)))/(32*R^6 - 24*R^4*斜^2 + 16*R^4*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) + 6*R^2*斜^4 - 4*R^2*斜^2*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) + 斜^4*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)), 斜*sqrt(4*R^2 - 斜^2)/(2*R))

全円の半径を表す式は長くなり，SymPy では簡約化できないが，式は正しい。
外円の直径が 10 寸，斜が 8 寸のとき，全円の直径は 15/2 = 7.5 寸である。

2res[1](R => 10//2, 斜 => 8) |> println

   15/2

function draw(R, 斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x) = (2*R*斜*(2*R*sqrt(128*R^8 - 128*R^6*斜^2 + 64*R^6*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) + 48*R^4*斜^4 - 32*R^4*斜^2*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) - 6*R^2*斜^6 + 8*R^2*斜^4*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) - 斜^6*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)) - 斜*(2*R - 斜)*(2*R + 斜)*(2*R^2 + sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)))/(32*R^6 - 24*R^4*斜^2 + 16*R^4*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) + 6*R^2*斜^4 - 4*R^2*斜^2*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) + 斜^4*sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)), 斜*sqrt(4*R^2 - 斜^2)/(2*R))
   @printf("外円の直径が %g，等斜の長さが %g のとき，全円の直径は %g である。\n", 2R, 斜, 2r)
   y = sqrt(R^2 - x^2)
   plot([-x, 0, x], [y, -R, y], color=:green, lw=0.5)
   circle(0, 0, R, :blue)
   circle(0, R - r, 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(0, R - r, "全円:r,(0,R-r)", :red, :center, delta=-delta/2)
       point(x, sqrt(R^2 - x^2), "(x,sqrt(R^2-x^2)) ", :green, :right, :vcenter)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
   end
end;

draw(10/2, 8, true)