算額（その1175）

算額（その1175）

三八 諏訪大神社 熊谷市新堀 弘化4年(1847)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，外円，斜線4本

外円の中に原点の中心を原点として，x 軸と y 軸を線対称とする 4 本の斜線を引き，区画された領域に等円 4 個を容れる。外円の直径を 1.6 寸としたとき，等円が最も大きくなるときの等円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
とおき，以下の連立方程式を立てる。

include("julia-source.txt");

using SymPy

@syms R::positive, r::positive, x::positive, y::positive
eq1 = x^2 + y^2 - (R - r)^2
eq2 = r^2*(R^2 + y^2) - x^2*y^2;

eq1 は x^2 = (R - r)^2 - y^2 なので，eq2 の x^2 に代入して整理する。

eq2 = eq2(x^2 => (R - r)^2 - y^2) |> expand
eq2 |> println

   R^2*r^2 - R^2*y^2 + 2*R*r*y^2 + y^4

eq2 は，R, y で r が決まるという式なので，r について解けば R, y の関数になる。

f = solve(eq2, r)[1]
f |> println

   y*(R - y)/R

これは単純に R が定数のときは放物線の式であり，y がどのような値を取ると f(y) が最大になるかは計算しなくてもわかる。
たとえば R = 1.6/2 のとき y = 0.4 のとき r = 0.2 が最大値になる。

pyplot(size=(300, 200), grid=true, aspectratio=:none, label="")
plot(f(R => 1.6/2), xlims=(0, 0.8), xlabel="y", ylabel="r") 

r が最大になるときの y，そして y がその値を取るときの r は以下のようにして求める。

g = diff(f, y)
y0 = solve(g, y)[1]
y0 |> println

   R/2

r0 = f(y => y0)
r0 |> println

   R/4

R = 1.6/2 寸 のとき，y0 = 0.4 寸, r0 = 0.2 寸ゆえ，等円の直径は 0.4 寸である。

@syms x0, y0, x, y, r, R
y0 = sqrt(R^2 - x0^2)
x = sqrt((R - r)^2 - y^2)
eq = dist2(-R, 0, x0, y0, x, y, r)
res = solve(eq, x0)[1] |> println

   (R*(8*R^4 - 16*R^3*r + 8*R^3*sqrt(R^2 - 2*R*r + r^2 - y^2) + 8*R^2*r^2 - 8*R^2*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 8*R^2*y^2 + 4*R*r*y^2 - 4*R*y^2*sqrt(R^2 - 2*R*r + r^2 - y^2) - r^4 + 2*r^2*y^2) - 4*sqrt(2)*r*y*sqrt(R^3*(4*R^3 - 8*R^2*r + 4*R^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + 5*R*r^2 - 4*R*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 3*R*y^2 - r^3 + r^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + r*y^2 - y^2*sqrt(R^2 - 2*R*r + r^2 - y^2))))/(8*R^4 - 16*R^3*r + 8*R^3*sqrt(R^2 - 2*R*r + r^2 - y^2) + 12*R^2*r^2 - 8*R^2*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 4*R^2*y^2 - 4*R*r^3 + 4*R*r^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + r^4)

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y, r) = (R/2, R/4)
   x = sqrt((R - r)^2 - y^2)
   t = sqrt(R^2 - 2*R*r + r^2 - y^2)
   u1 = R*(8R^4 - 16R^3*r + 8R^3*t + 8R^2*r^2 - 8R^2*r*t - 8R^2*y^2 + 4R*r*y^2 - 4R*y^2*t - r^4 + 2r^2*y^2)
   u2 = 4√2r*y*sqrt(R^3*(4R^3 - 8R^2*r + 4R^2*t + 5R*r^2 - 4R*r*t - 3R*y^2 - r^3 + r^2*t + r*y^2 - y^2*t))
   v = 8R^4 - 16R^3*r + 8R^3*t + 12R^2*r^2 - 8R^2*r*t - 4R^2*y^2 - 4R*r^3 + 4R*r^2*t + r^4
   x0 = (u1 - u2)/v
   y0 = sqrt(R^2 - x0^2)
   plot()
   circle(0, 0, R, :red)
   circle4(x, y, r, :blue)
   segment(-R, 0, x0, y0, :green)
   segment(-R, 0, x0, -y0, :green)
   segment(R, 0, -x0, y0, :green)
   segment(R, 0, -x0, -y0, :green)
   #=
   abline(0, y, -y/R, -R, R, :green)
   abline(0, y, y/R, -R, R, :green)
   abline(0, -y, -y/R, -R, R, :green)
   abline(0, -y, y/R, -R, R, :green)
   =#
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
       point(0, y, "y", :blue, :center, delta=-delta/2)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

draw(1.6/2, true)

算額（その1174）

算額（その1174）

一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円1個，二等辺三角形，正方形

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

図形を時計回りに倒し，二等辺三角形の右側の等辺を水平にして考える。
二等辺三角形の頂角をθ，その半分を θ/2
正方形の一辺の長さを s
円の半径を r
とおいて，以下の連立方程式を解く。

メモ：方程式で tan() を使って書くと計算時間がかかり，有限の時間内に答えが求まらない。

include("julia-source.txt");

using SymPy

@syms s::positive, r::positive,
     a::positive, b::positive,
     θ::positive;

eq1 = s*tand(θ/2) + s + s/tand(θ) - a
eq2 = s + (a - s - s*tand(θ/2)) - sqrt(s^2 + (a - s - s*tand(θ/2))^2) - 2r
res = solve([eq1, eq2], (a, θ))[1]

   ((2*r*(r - s) + s*(2*r - s)*(tan(atan(s*(1 - s/(2*r))/(r - s))/2) + 1))/(2*r - s), 180*atan(s*(2*r - s)/(2*r*(r - s)))/pi)

a を求める計算式は複雑であるが，実際の計算には支障がない。

正方形の一辺の長さが 3 寸，円の直径が 2 寸のとき，二等辺三角形の等辺の長さは 8 寸である。

res[1](s => 3, r => 1).evalf() |> println

   8.00000000000000

function draw(s, r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, θ) = ((2*r*(r - s) + s*(2*r - s)*(tan(atan(s*(1 - s/(2*r))/(r - s))/2) + 1))/(2*r - s), 180*atan(s*(2*r - s)/(2*r*(r - s)))/pi)
   @printf("a = %g;  θ = %g\n", a, θ)
   b = s*tand(θ/2)
   (x, y) = float.(intersection(0, 0, b, s, a, 0, b + s, s))
   @printf("b = %g;  s = %g;  x = %g;  y = %g\n", b, s, x, y)
   plot()
   rect(b, 0, b + s, s, :blue)
   circle(b + s + r, r, r, :red)
   plot!([0, a, x, 0], [0, 0, y, 0], color=:green, lw=0.5)
   segment(x/2, y/2, a, 0, :gray)
   circle(a, 0, 0.1a, beginangle=180 - θ, endangle=180, lw=2)
   circle(a, 0, 0.12a, :magenta, beginangle=180 - θ/2, endangle=180, lw=2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       #point(x, y, " (x,y)", :green, :left, :vcenter)
       point(a, 0, "a", :green, :center, delta=-delta)
       point(b, 0, "b", :blue, :center, delta=-delta)
       point(b + s, 0, "b+s", :blue, :center, delta=-delta)
       point(b + s + r, 0, "b+s+r", :blue, :center, delta=-delta)
       point(b + s + r, r, "", :blue, :center, delta=-delta)
       point(0.88a, delta/2, "θ/a ", :magenta, :right, :bottom, delta=delta, mark=false)
       point(0.96a, -0.05a*tan(θ), "θ ", :red, :right, :bottom, delta=2delta, mark=false)
       ylims!(-8delta, y + 3delta)
   end
end;

 

算額（その1173）

算額（その1173）

『算法用学精 巻之二』 文久2年(1862)
http://www.wasan.jp/terasima/terasima3.pdf
キーワード：円4個，外円，弦

外円の中に，水平な弦，大円 1 個，小円 2 個を容れる。外円の直径が 1 寸，大円の直径が 0.81 寸のとき，小円の直径を求めよ。

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

図形的には，算額（その750）と同等である。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative;

eq1 = x1^2 + (2r2 - R + r1)^2 - (R - r1)^2
eq2 = x2^2 + (3r2 - R)^2 - (R - r2)^2
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r2, x2, x1))[1];

SymPy では十分な簡約化ができないので，手動で簡約化する。

@syms tR, tr1
res[1](√R => tR, √r1 => tr1) |> simplify |> (x -> x(tR => √R, tr1 => √r1)) |> println

   2*sqrt(r1)*(sqrt(R) - sqrt(r1))

res[2](√R => tR, √r1 => tr1) |> simplify |> (x -> x(tR => √R, tr1 => √r1)) |> println

   2*sqrt(2)*r1^(1/4)*(R^(3/2) + 5*sqrt(R)*r1 - 4*R*sqrt(r1) - 2*r1^(3/2))/sqrt(R^(3/2) + 3*sqrt(R)*r1 - 3*R*sqrt(r1) - r1^(3/2))

res[3](√R => tR, √r1 => tr1) |> simplify |> (x -> x(tR => √R, tr1 => √r1)) |> println

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

r2 = 2*sqrt(r1)*(sqrt(R) - sqrt(r1))
t = sqrt(R^(3/2) + 3*sqrt(R)*r1 - 3*R*sqrt(r1) - r1^(3/2))
x2 = 2*sqrt(2)*r1^(1/4)*(R^(3/2) + 5*sqrt(R)*r1 - 4*R*sqrt(r1) - 2*r1^(3/2))/t
x1 = 2*sqrt(2)*r1^(1/4)*t

小円の半径 r2 は 2√r1*(√R - √r1) で求めることができる。
外円，大円の直径をそれぞれ 1 寸，0.81 寸のとき，小円の直径は 0.18 寸である。

R = 1/2
r1 = 0.81/2
r2 = 2√r1*(√R - √r1)

   0.0900000000000001

function draw(R, r1, more, plot_flag)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2*sqrt(r1)*(sqrt(R) - sqrt(r1))
   t = sqrt(R^(3/2) + 3*sqrt(R)*r1 - 3*R*sqrt(r1) - r1^(3/2))
   u = 2*sqrt(2)*r1^(1/4)
   x1 = u*t
   x2 = u*(R^(3/2) + 5*sqrt(R)*r1 - 4*R*sqrt(r1) - 2*r1^(3/2))/t
   @printf("外円の直径が %g，大円の直径が %g のとき，小円の直径は %g である。\n", 2R, 2r1, 2r2)
   y = 2r2 - R
   x = sqrt(R^2 - y^2)
   plot()
   circle(0, 0, R, :magenta)
   circle(x1, y + r1, r1)
   circle(x2, y + r2, r2, :blue)
   circle(0, r2 - R, r2, :blue)
   segment(-x, y, x, y, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
       point(x1, y + r1, "大円:r1,(x1,2r2-R+r1)", :red, :center, delta=-delta/2)
       point(x2, y + r2, " 小円:r2,(x2,3r2-R)", :blue, :left, :vcenter)
       point(0, r2 - R, " 小円:0,(0,r2-R)", :blue, :left, :vcenter)
   end
end;

算額（その1172）

算額（その1172）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，外円

外円の中に，甲円 1 個，乙円 2 個，丙円 2 個，丁円 2 個を容れる。丁円の直径が 1 寸のとき，甲円の直径はいかほどか。

注：「問」には明示的に書かれておらず，『埼玉の算額』の図では下にある乙円の中心が甲円の円周上にないように見えるが，実際は「円周上にある」としないと（条件不足で）連立方程式が解けない。この制約がないと甲円の大きさに制限がつかない。甲円の大きさに制限をつけるのはおかしい気もする。甲円の大きさに制約をつけたとき，実は丙円と丁円が同じ大きさになる。「問」でも『埼玉の算額』の図でも，丙円と乙円の大きさは違うと示唆している。しかし，どの程度違うかを明示しないと甲円の大きさは決まらない。

以下では，条件不足を解消するために，丙円の大きさも指定する。

外円の半径と中心座標を R, (0, 0); R = 2r2
甲円の半径と中心座標を r1, (0, R - r1); 2r1 = 3r2
乙円の半径と中心座標を r2, (0, r2), (0, -r2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4); y4 = 0
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::negative,
     r4::positive, x4::positive;

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

   (3*r3/2 + 3*r4/2 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)/2, (3*r3 + 3*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))*(3*r3 + 5*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))/(4*(3*r3 + r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))), 3*r3/4 + 3*r4/4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)/4, sqrt(2)*sqrt(r3)*sqrt(-r3 + 3*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)), 3*r3/2 - 3*r4/2 - sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)/2, sqrt(2)*sqrt(r4)*sqrt(3*r3 + 5*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))/2)

@syms t
println("t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)")
var = ("R", "r1", "r2", "x3", "y3", "x4")
for i = 1:6
   @printf("%s = %s\n", var[i], res[i](sqrt(9r3^2 - 2r3*r4 + 9r4^2) => t))
end

   t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)
   R = 3*r3/2 + 3*r4/2 + t/2
   r1 = (3*r3 + 3*r4 + t)*(3*r3 + 5*r4 + t)/(4*(3*r3 + r4 + t))
   r2 = 3*r3/4 + 3*r4/4 + t/4
   x3 = sqrt(2)*sqrt(r3)*sqrt(-r3 + 3*r4 + t)
   y3 = 3*r3/2 - 3*r4/2 - t/2
   x4 = sqrt(2)*sqrt(r4)*sqrt(3*r3 + 5*r4 + t)/2

丙円と丁円の直径が等しい 2r3 = 2r4 = 1 寸のときにのみ，算額の「答」通り，甲円の直径 2r1 = 3.75 寸になる。

r3 = r4 = 1/2
t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)
r1 = (3*r3 + 3*r4 + t)*(3*r3 + 5*r4 + t)/(4*(3*r3 + r4 + t))
2r1

   3.75

function draw(r4, r3, more, plot_flag)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)
   R = 3*r3/2 + 3*r4/2 + t/2
   r1 = (3*r3 + 3*r4 + t)*(3*r3 + 5*r4 + t)/(4*(3*r3 + r4 + t))
   r2 = 3*r3/4 + 3*r4/4 + t/4
   x3 = sqrt(2)*sqrt(r3)*sqrt(-r3 + 3*r4 + t)
   y3 = 3*r3/2 - 3*r4/2 - t/2
   x4 = sqrt(2)*sqrt(r4)*sqrt(3*r3 + 5*r4 + t)/2
   if !plot_flag
       @printf("丁円の直径が %g のとき，甲円の直径は %g である。\n", 2r4, 2r1)
       @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g\n", R, r1, r2, r3, x3, y3, r4, x4)
   end
   p1 = plot(showaxis=!plot_flag)
   circle(0, 0, R, :magenta)
   circle(0, R - r1, r1)
   circle22(0, r2, r2, :blue)
   circle2(x3, y3, r3, :orange)
   circle2(x4, 0, r4, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
       !plot_flag && point(0, R - r1, @sprintf("甲円の直径 = %g", 2r1), :red, :center, :bottom, delta=delta/2)
       point(0, R - r1, "r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(0, r2, "乙円:r2,(0,r2)", :blue, :center, delta=-delta/2)
       point(0, -r2, "", :blue, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3\n(x3,y3)", :orange, :center, delta=-delta/2)
       point(x4, 0, "丁円:r4\n(x4,0)", :green, :center, delta=-delta/2)
   end
   return p1
end;

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

算額（その1171）

算額（その1171）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円1個，正三角形，菱形，斜線

正三角形の中に斜線と菱形と円を容れる。正三角形の一辺の長さを 8 寸としたとき，斜線の長さはいかほどか。

正三角形の一辺の長さを a
円の半径と中心座標を r, (x, r)
菱形の頂点の座標を (c, 2r), (a - 2r/√3, 0)
斜線と正三角形の底辺との交点座標を (b, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r::positive, x::positive, a::positive, b::positive, c::positive
c = a - 2r/√Sym(3)
eq1 = dist2(0, 0, a/2, √Sym(3)a/2, x, r, r)
eq2 = dist2(b, 0, a/2, √Sym(3)a/2, x, r, r)
eq3 = dist2(2c - a, 0, c, 2r, x, r, r)
res = solve([eq1, eq2, eq3], (b, x, r))[1]

   (5*a/8, 3*a/8, sqrt(3)*a/8)

斜線と正三角形の底辺との交点の座標は，(5a/8, 0) である。
正三角形の頂点 (a/2, √3a/2) からの線分の長さは 
sqrt((5a/8 - a/2)^2 + (√3a/2)^2) = 7a/8 である。
正三角形の一辺の長さが 8 寸のときは，斜線の長さは 7 寸である。

図を描いて，正三角形の辺に平行で円の中心を通る平行線や接線を描くと一辺が2の小さな正三角形ができる。三角形 ABC の斜辺の長さ AB を求めることになる。AC = 1, BC = 4√3 なのでピタゴラスの定理から AB = sqrt(AC^2 + BC^2) = sqrt(1 + 48) = sqrt(49) = 7 となる。

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

   a = 8;  b = 5;  c = 6;  x = 3;  r = 1.73205

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, x, r) = (5*a/8, 3*a/8, sqrt(3)*a/8)
   c = a - 2r/√3
   @printf("正三角形の一辺の長さが %g のとき，斜線の長さは %g である。\n", a, 7a/8)
   @printf("a = %g;  b = %g;  c = %g;  x = %g;  r = %g\n", a, b, c, x, r)
   plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:blue, lw=0.5)
   circle(x, r, r)
   segment(a/2, √3a/2, b, 0)
   segment(2c - a, 0, c, 2r, :green)
   segment(c, 2r, a - c, 2r, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a/2, √3a/2, " (a/2,√3a/2)", :blue, :left, :vcenter)
       point(c, 2r, " (c,2r)", :green, :left, :vcenter)
       point(a - c, 2r, "(a-c,2r) ", :green, :right, :vcenter)
       point(x, r, "(x,r)", :red, :center, delta=-delta/2)
       point(2c - a, 0, "2c-a", :green, :center, delta=-delta/2)
       point(b, 0, "b", :black, :center, delta=-delta/2)
       point(a, 0, "a", :blue, :center, delta=-delta/2)
       ylims!(-5delta, √3a/2 + 2delta)
   end
end;

算額（その1170）

算額（その1170）

　七九 春日部市 香取神社 明治9年(1876)
　九九 春日部市小渕 観音院 明治30年(1897)
一〇三 春日部市 東福寺 明治40年(1907)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：面積

面積が 120408472.1265051159 坪の円形の土地がある。円の直径を求めよ。

注：10間四方=1坪

面積を「円法 = 0.79」で割り，平方根を求めればよい。

数値演算的には問題はないが，答えが 123456.78899999999間 となり，出題者が期待した 12345.6789 にならないのは元の面積の数値が精度不足なのである。

sqrt(120408472.1265051159/0.79)*10  # 10間四方=1坪

   123456.78899999999

sqrt(big"120408472.1265051159"/0.79)*10

   123456.7889999999972240087150564101168588130753116810662338684335743893401965257

あまり意味はないが，以下のように長い数値だと 123456.789 に極めて近い数値が得られる。

sqrt(big"120408472.12650512131489653123169070170206396142020821571350097656251"/0.79)*10

   123456.7890000000000000000000000000000000000000000000000000000000000051265823247

ところで，「円法 = 0.79」はなにかといえば，一辺の長さが a の正方形の面積 s1 = a^2 と，正方形の中に内接する円，すなわち直径が a の円の面積 s2 = (a/2)^2 * 円周率 = a^2 * (円周率/4) の比である。
s1/s2 = (a^2) / (a^2 * (円周率/4)) = 円周率/4 = 円法
円法 = 0.79 ということは，円周率 = 4*0.79 = 3.16 である。当時は円周率は「円率」と呼ばれていた。
今はよく使われている 3.14 ではないかといえば，もともとは 3.14 に近い近似値は知られていたのであろうが pi/4 ≒ 0.7853981633974483 の小数点以下 3 桁目で四捨五入して 0.79 としたのであるが，0.7853... を 0.79 はちょっと大雑把すぎたのである。

現在用いられている円周率 π を用いたときに問題文を書き直すと以下のようになる。

pi/4

   0.7853981633974483

big"123456.789"^2*pi/4

   1.197070795767721080355012226425152943759227052416650126777784107916604885390605e+10

sqrt(big"119707079.576772108"/(pi/4))*10

   123456.7890000000023879722548741154634802177298661877133836534897753633868719298

面積が 1億1970万7079歩5分7厘6毛7糸7忽2微1繊0沙8塵の円形の土地の直径はいかほどか。
答　123456.789 間

算額（その1169）

算額（その1169）

　九九 春日部市小渕 観音院 明治30年(1897)
一〇三 春日部市 東福寺 明治40年(1907)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：面積

1,5241,5787.50190521 坪の土地がある。
これを正方形にしたとき，その一辺はいかほどか。
注：10 間四方を 1 坪とする。

平方根を求めればよい。
そのまま計算しても精度上の問題はない。
12万3456.789間

sqrt(1_5241_5787.50190521)*10  # 10間四方=1坪

   123456.78899999999

123456.789^2

   1.5241578750190521e10

sqrt(big"1_5241_5787.50190521")*10

   123456.7890000000000000000000000000000000000000000000000000000000000000000000004

big"123456.789"^2

   1.524157875019052100000000000000000000000000000000000000000000000000000000000004e+10

 

算額（その1168）

算額（その1168）

十七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

(10) 京都市中京区三条大宮西二筋目下ル 武信稲荷神社 嘉永6年(1853)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

キーワード：円6個，二等辺三角形

二等辺三角形と甲円，乙円，丙円が 2 個ずつある。乙円の直径が 65 寸のとき，丙円の直径はいくつか。

算額（その360）に丙円を加えたものである。

甲円の直径，中心座標を r1, (0, r1) および (0, 3r1)
乙円の直径，中心座標を r2, (x2，y2)
丙円の半径，中心座標を r3, (x3, y3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive, a::positive, r3::positive, x3::positive;

eq1 = x2^2 + (y2 - 3r1)^2 - (r1 - r2)^2
eq2 = a^2 + (4r1)^2 - (3a)^2
eq3 = 3(r1 - 2r2) - r1
eq4 = 4r1*(y2 - 3r1) - a*x2
eq5 = x3^2 + r1^2 - (r1 + r3)^2
eq6 = dist2(0, 4r1, a, 0, x3, 2r1, r3)
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r3, x3, r1, a, x2, y2))[2]

   (6*r2*(7 - 4*sqrt(3)), -6*sqrt(2)*r2*(35 - 21*sqrt(3))/7, 3*r2, 3*sqrt(2)*r2, 4*sqrt(2)*r2/3, 29*r2/3)

丙円の半径 r3 は乙円の半径 r2 の 6(7 - 4√3) である。
よって，乙円の直径が 65 寸なら，丙円の直径は 65*6(7 - 4√3) = 28.00074019255158 寸である。

65*6(7 - 4√3)

   28.00074019255158

ちなみに，甲円の直径は乙円の直径の 3 倍である。

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

  r2 = 32.5;  r3 = 14.0004;  x3 = 54.0933;  r1 = 97.5;  a = 137.886;  x2 = 61.2826;  y2 = 314.167

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, x3, r1, a, x2, y2) = (6*r2*(7 - 4*sqrt(3)), -6*sqrt(2)*r2*(35 - 21*sqrt(3))/7, 3*r2, 3*sqrt(2)*r2, 4*sqrt(2)*r2/3, 29*r2/3)
   @printf("乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r2, 2r3)
   @printf("r2 = %g;  r3 = %g;  x3 = %g;  r1 = %g;  a = %g;  x2 = %g;  y2 = %g\n", r2, r3, x3, r1, a, x2, y2)
   plot([a, 0, -a, a], [0, 4r1, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, 3r1, r1)
   circle2(x2, y2, r2, :blue)
   circle2(x3, 2r1, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, 4r1, " 4r1", :red, :left, :bottom, delta=delta/2)
       point(0, 3r1, " 3r1", :red, :left, :vcenter)
       point(0, 2r1, " 2r1", :red, :left, :bottom, delta=delta/2)
       point(0, r1, "甲円:r1,(0,r1)", :red, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :center, delta=-delta/2)
       point(x3, 2r1, "丙円:r3 \n(x3,2r1)", :green, :right, delta=-2delta)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
   end
end;

和算の心（その006）

和算の心（その006）

直線上の 3 個の円

直線上に，甲円，乙円，丙円が互いに接しており，この順序で載っている。甲円，乙円と直線の接点間の距離（以下では，水平距離と呼ぶ）を求めよ。

甲円，乙円，丙円の半径を r1，r2，r3 とする。

1.　一般的な場合

次節で述べる特別な場合，すなわち 3 円が互いに接している（甲円と乙円も接している）場合よりも甲円と丙円が離れている場合である。

甲円と乙円，乙円と丙円の水平距離は ， である。（和算の心（その002））

結果として，甲円と丙円の水平距離は になる。

甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解き x3 を求める。

include("julia-source.txt");

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

SymPy では，これ以上は自動的に簡約化されないが，2√r2(√r1 + √r3) である。

res[2][2] |> simplify

   

res[2][2](√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

using LaTeXStrings
function draw(r1=3, r2=2, r3=1, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   x2 = 2√(r1*r2)
   x3 = 2√r2*(√r1 + √r3)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   p1 = plot(ylims=(-0.05r1, 1.05*2r1), showaxis=false)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   segment(-r1, 0, x3+r3, 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, r1, L"r_1,(0,r_1)", :red, :center, delta=-delta)
       point(x2, r2, L"r_2,(x_2,r_2)", :blue, :center, delta=-delta)
       point(x3, r3, L"r_3,(x_3,r_3)", :green, :center, delta=-delta)
       point(x3 - r3, 1.3*2r3, L"x_3 = 2\sqrt{r_2}(\sqrt{r_1} + \sqrt{r_3})", mark=false)
       point(x3 - r3, 1.15*2r3, L"x_2 = 2\sqrt{r_1 r_2})", mark=false)
   end
   display(p1)
end;

draw(3, 1, 2, true)

   r1 = 3;  r2 = 1;  x2 = 3.4641;  r3 = 2;  x3 = 6.29253

2.　 3 個の円が互いに接する場合

この状況では，乙円の半径には制約が付く（どのような値でも取れるわけではない）。そのため，甲円と乙円，乙円と丙円の水平距離は簡単な数式では表せない。しかし，甲円と丙円の水平距離は と簡単なままである。

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

include("julia-source.txt");

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

SymPy でこの連立方程式を解くと，r2 と x2 は複雑な式になる。

res[1][1]  # r2

   

√r1 => r1, √r3 => r3 と変換したうえで因数分解し，r1 => √r1, r3 => √r3 でもとに戻すと簡約化された式になる。

res[1][1](√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

res[1][2]  # x2

   

res[1][2](√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

x3 は簡潔な式のままである。

res[1][3]  #x3

   

別の関係式を使って連立方程式を解くと，特別なことをしないでも簡潔な式が得られる。

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive
eq1 = x2 - 2√(r1*r2)
eq2 = x3 - 2√(r1*r3)
eq3 = (x3 - x2) - 2√(r2*r3)
res = solve([eq1, eq2, eq3], (r2, x2, x3));

res[1][1] |> display
res[1][2] |> display
res[1][3] |> display

   

   

   

デカルトの円定理による場合

デカルトの円定理は，半径の関係のみである。半径を求めた後で，位置を求めればよい。

@syms r1::positive, r2::positive, r3::positive
r2 = 1/(1/r1 + 1/r3 + 2√(1/(r1*r3)))

   

これも，√r1 => r1, √r3 => r3 と変換したうえで因数分解し，r1 => √r1, r3 => √r3 でもとに戻すと簡約化された式になる。

r2(√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

using LaTeXStrings
function draw(r1=3, r3=1, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   (r2, x2, x3) = (r1*r3/(sqrt(r1) + sqrt(r3))^2, 2*r1*sqrt(r3)/(sqrt(r1) + sqrt(r3)), 2*sqrt(r1)*sqrt(r3))
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   p1 = plot(ylims=(-0.05r1, 1.05*2r1), showaxis=false)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   segment(-r1, 0, x3+r3, 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, r1, L"r_1,(0,r_1)", :red, :center, delta=-delta)
       point(x2, r2, L"r_2,(x_2,r_2)", :blue, :center, delta=-delta)
       point(x3, r3, L"r_3,(x_3,r_3)", :green, :center, delta=-delta)
       point(x3 - r3, 1.45*2r3, L"r_2 = \frac{r_1 r_3}{(\sqrt{r1} + \sqrt{r3})^2}", mark=false)
       point(x3 - r3, 1.15*2r3, L"x_3 = 2\sqrt{r_1 r_3}", mark=false)
   end
   display(p1)
end;

draw(3, 2, true)

   r1 = 3;  r2 = 0.606123;  x2 = 2.69694;  r3 = 2;  x3 = 4.89898

 

3.　3 個の円が互いに接し，両端の円の半径が等しい場合

乙円，丙円の中心座標を (x2, r2), (x3, r3) とするが，2√(r1*r2) + 2√(r2*r3) = 2√(r1*r3) となる。
x2 = r1, x3 = 2r1 は自明である。
さらに，r1 = r3 である。
以下の方程式を解くと r2 = r1/4 であることがわかる。

include("julia-source.txt");

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive
x2 = r1
eq1 = r1^2 + (x2 - r2)^2 - (r1 + r2)^2
solve(eq1, r2)[1]  # r2

   

または，前節で得られた r2 の式を r3 = r1 として簡約化すればよい。

r3 = r1
r1*r3/(sqrt(r1) + sqrt(r3))^2  # r2

   

または，もう少し前の方程式をr3 = r1 として解いてもよい。

r3 = r1
eq = 2√(r1*r2) + 2√(r2*r3) - 2√(r1*r3) 
solve(eq, r2)[1]  # r2

   

using LaTeXStrings
function draw(r1, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   r3 = r1
   x3 = 2r1
   r2 = r1/4
   x2 = r1
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   p1 = plot(ylims=(-0.05r1, 1.05*2r1), showaxis=false)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   segment(-r1, 0, x3 + r3, 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, r1, L"r_1,(0,r_1)", :red, :center, delta=-delta)
       point(x2, r2, L"r_2,(x_2,r_2)", :blue, :center, delta=-delta)
       point(x3, r3, L"r_3,(x_3,r_3)", :green, :center, delta=-delta)
       point(x3 - r3, 1.15*2r3, L"$r_3 = r_1,\ x_3 = 2r_1$", mark=false)
       point(x3 - r3, 1.3*2r3, L"$r_2 = \frac{r_1}{4},\ x_2 = r_1$", mark=false)
   end
   display(p1)
end;

draw(3, true)

   r1 = 3;  r2 = 0.75;  x2 = 3;  r3 = 3;  x3 = 6

 

算額（その1167）

算額（その1167）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

長方形内に交差する 2 個の楕円を容れる。楕円の長径が 4 寸，短径が 2.4 寸，長方形の短辺が 3 寸のとき，甲斜の長さはいかほどか。
注：甲斜とは，原点対称の 2 個の楕円の交点間の距離の長い方である。

図に示す x1, y1, x2, y2, x01, y01, x02, y02, x03, y03 を未知数として以下の連立方程式を立て，それを解く。

甲斜の長さは 2sqrt(x03^2 + y03^2) で求める。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, side_a::positive, side_b::positive, K::positive,
     x1::positive, y1::positive, x2::positive, y2::negative,
     x01::positive, y01::positive, 
     x02::positive, y02::positive,
     x03::positive, y03::positive
side_b = sqrt((x2 - x1)^2 + (y1 - y2)^2)
side_a = sqrt((x2 + x1)^2 + (y2 + y1)^2)
eq1 = x01^2/a^2 + y01^2/b^2 - 1
eq2 = -b^2*x01/(a^2*y01) - (y2 - y1)/(x2 - x1)
eq7 = (y2 - y01)/(x2 - x01)  - (y2 - y1)/(x2 - x1)
eq3 = x02^2/a^2 + y02^2/b^2 - 1
eq4 = -b^2*x02/(a^2*y02) - (y2 + y1)/(x2 + x1)
eq8 = (y2 - y02)/(x2 - x02) - (y2 + y1)/(x2 + x1)
eq5 = (4a^2 + 4b^2) - (side_a^2 + side_b^2)  # formula-89
eq6 = (y2 - y1)/(x2 - x1) * (y2 + y1)/(x2 + x1) + 1
eq9 = side_b - K
eq10 = y03/x03 - (y2 + y1)/(x2 + x1)
eq11 = x03^2/a^2 + y03^2/b^2 - 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 Float64.(v), r.f_converged
end;

function H(u)
   (x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = u
   return [
       -1 + y01^2/b^2 + x01^2/a^2,  # eq1
       -(-y1 + y2)/(-x1 + x2) - b^2*x01/(a^2*y01),  # eq2
       -1 + y02^2/b^2 + x02^2/a^2,  # eq3
       -(y1 + y2)/(x1 + x2) - b^2*x02/(a^2*y02),  # eq4
       4*a^2 + 4*b^2 - (-x1 + x2)^2 - (x1 + x2)^2 - (y1 - y2)^2 - (y1 + y2)^2,  # eq5
       1 + (-y1 + y2)*(y1 + y2)/((-x1 + x2)*(x1 + x2)),  # eq6
       -(-y1 + y2)/(-x1 + x2) + (-y01 + y2)/(-x01 + x2),  # eq7
       -K + sqrt((-x1 + x2)^2 + (y1 - y2)^2),  # eq9
       -(y1 + y2)/(x1 + x2) + y03/x03,  # eq10
       -1 + y03^2/b^2 + x03^2/a^2,  # eq11
   ]
end;

(a, b, K) = (4/2, 2.4/2, 3)
iniv = BigFloat[0.614, 2.254, 2.315, -0.228, 1.789, 0.539, 1.505, -0.788, 1.315, 0.91]
res = nls(H, ini=iniv)

   ([0.6329571688388319, 2.244853051407938, 2.3204571688388325, -0.2355388027151154, 1.8516704311458714, 0.45351293387424213, 1.5000000000000002, -0.793725393319377, 1.322875655532295, 0.9000000000000001], true)

(x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = res[1]
甲斜 = 2sqrt(x03^2 + y03^2)

   3.1999999999999993

甲斜は 3.2 寸である。

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

  長辺 = 3.57211;  短辺 = 3;  x1 = 0.632957;  y1 = 2.24485;  x2 = 2.32046;  y2 = -0.235539;  x01 = 1.85167;  y01 = 0.453513;  x02 = 1.5;  y02 = -0.793725;  x03 = 1.32288;  y03 = 0.9

function draw(p, q, a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (4/2, 2.4/2)
   (x1, y1) = (0.614, 2.254)
   (x2, y2) = (2.315, -0.228)
   (x01, y01) = (1.789, 0.539)
   (x02, y02) = (1.505, -0.788)
   (x1, y1, x2, y2, x01, y01, x02, y02, x03, y03) = res[1]
   side_b = sqrt((x2 - x1)^2 + (y1 - y2)^2)
   side_a = sqrt((x2 + x1)^2 + (y2 + y1)^2)
   @printf("甲斜 = %g;  長辺 = %g;  短辺 = %g\n", 2sqrt(x03^2 + y03^2), side_a, side_b)
   @printf("長辺 = %g;  短辺 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g;  x01 = %g;  y01 = %g;  x02 = %g;  y02 = %g;  x03 = %g;  y03 = %g\n", side_a, side_b, x1, y1, x2, y2, x01, y01, x02, y02, x03, y03)
   plot([x1, -x2, -x1, x2, x1], [y1, -y2, -y1, y2, y1], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   ellipse(0, 0, a, b, φ= 2atand(y03/x03), color=:red)
   segment(x03, y03, -x03, -y03, :green)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x1, y1, " (x1,y1)", :blue, :left, :vcenter)
       point(x2, y2, " (x2,y2)", :blue, :left, :vcenter)
       point(-x2, -y2, " (-x2,-y2)", :blue, :left, :vcenter)
       point(-x1, -y1, " (-x1,-y1)", :blue, :left, :vcenter)
       point(x01, y01, " (x01,y01)", :red, :left, :vcenter)
       point(x02, y02, " (x02,y02)", :red, :left, :vcenter)
       point(x03, y03, " (x03,y03)", :green, :left, :vcenter)
       point(-x03, -y03, " (-x03,-y03)", :green, :left, :vcenter)
       xlims!(-x2 - 5delta, x2+10delta)
   end
end;

算額（その1166）

算額（その1166）

六四 加須市不動岡 総願寺 慶応2年(1866)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：楕円4個，長方形，対角線

長方形の中に対角線を引き，分割された区画に，合同な楕円を 4 個容れる。楕円の短径が 1 寸のとき，長径はいかほどか。

注：対角線を引いてできる三角形のうち，長方形の短辺を一辺とする三角形は正三角形である。

二等辺三角形に内接する楕円の長径，短径については算法助術の公式97による（「和算の心（005）」参照）。

長方形の中心を原点とし，長方形の長辺，短辺を 2a, 2b とする。
正三角形の底辺と高さは，2b, √3b
二等辺正三角形の底辺と高さは，2a, b
以上に基づき，
(1) 2 つの楕円の長径を求め，それが等しいこと
(2) a, b の関係式
の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, side_a::positive, height::positive, q::positive
p = side_a*sqrt((height - 2*q)/height)/2
p1 = p(side_a => 2b, height => √Sym(3)b);
p1 |> println

   3^(3/4)*sqrt(b)*sqrt(sqrt(3)*b - 2*q)/3

p2 = p(side_a => 2a, height => b)
p2 |> println

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

eq1 = p1 - p2
eq2 = 4b^2 - (a^2 + b^2);
res = solve([eq1, eq2], (a, b))[1]

   (q*(-1 + 3*sqrt(3)), q*(9 - sqrt(3))/3)

@syms q
side_a = 2q*(9 - √Sym(3))/3
height = √Sym(3)side_a/2
p = side_a*sqrt((height - 2*q)/height)/2 |> simplify
p |> println

   3^(3/4)*q*sqrt(-12 + 10*sqrt(3))/3

手計算で簡約化すると，以下のように若干簡単になる。

p = q*sqrt(10 - 4*sqrt(Sym(3)))
p |> println

   q*sqrt(10 - 4*sqrt(3))

長半径 p は短半径 q の sqrt(10 - 4√3) 倍である。
短径が 1 寸のとき，長径は sqrt(10 - 4√3) = 1.7526542071168778 寸である。

function draw(q, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (q*(-1 + 3√3), q*(9 - √3)/3)
   p = q*sqrt(10 - 4√3)
   @printf("a = %g;  b = %g;  q = %g;  p = %g\n", a, b, q, p)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   segment(a, b, -a, -b, :magenta)
   segment(-a, b, a, -b, :magenta)
   ellipse(a - q, 0, q, p, color=:red)
   ellipse(q - a, 0, q, p, color=:red)
   ellipse(0, q - b, p, q, color=:green)
   ellipse(0, b - q, p, q, color=:green)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a - q, 0, "(a-q,0)", :red, :center, delta=-delta/2)
       point(0, b - q, "(0,b-q,0)", :green, :center, delta=-delta/2)
       point(a, 0, " a", :red, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(-0.3a, -0.5b, @sprintf(" a=%g, b=%g\np=%g, q=%g", a, b, p, q), :black, mark=false)
   end
end;

和算の心（その005）

和算の心（その005）

三角形に内接する楕円 -- 算法助術の公式97

1.　不等辺三角形の場合

不等辺三角形の場合について一般解が示されている。
和算の慣例として，楕円の長軸（長径），短軸（短径）を p, q としているが，ここでは長半径，短半径を使用することにする。

include("julia-source.txt");

@syms a, b, c, h, p, q
eq97 = -(b^2 - c^2)^2*h*q^2 + (b^2 - c^2)^2*q^3 + a^4*h*(2h - q)^2 - a^4*q*(2h - q)^2 - a^2*h*p^2*(2h - q)^2
eq97 = eq97(q => 2q, p => 2p)
eq_97 = eq97/4 |> simplify
eq_97 |> println

   a^4*h*(h - q)^2 - 2*a^4*q*(h - q)^2 - 4*a^2*h*p^2*(h - q)^2 - h*q^2*(b^2 - c^2)^2 + 2*q^3*(b^2 - c^2)^2

2.　二等辺三角形の場合

b = c の場合は，簡単になる。

eq97_2 = eq97(c => b) |> factor
eq97_2 |> println

   -4*a^2*(-h + q)^2*(-a^2*h + 2*a^2*q + 4*h*p^2)

更に a ≠ 0, h ≠ q なので非常に簡単な式になる。

eq97_3 = eq97_2/(4a^2*(q - h)^2)
eq97_3 |> println

   a^2*h - 2*a^2*q - 4*h*p^2

2.1.　a, h, q を与えて p を知るとき

res_p = solve(eq97_3, p)[2]
res_p |> println

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

a = 5; h = 4; q = 1 のとき

ans_p = res_p(a => 5, h => 4, q => 1)
ans_p |> println
ans_p.evalf() |> println

   5*sqrt(2)/4
   1.76776695296637

2.2.　a, h, p を与えて q を知るとき

res_q = solve(eq97_3, q)[1]
res_q |> println

   h/2 - 2*h*p^2/a^2

a = 5; h = 4; q = 2 のとき

ans_q = res_q(a => 5, h => 4, p => 2)
ans_q |> println
ans_q.evalf() |> println

   18/25
   0.720000000000000

3.　関数として定義する

function formula97(a, h; p=0, q = 0)
   p == 0 && q == 0 && return "error"
   p != 0 && return h/2 - 2h*p^2/a^2
   q != 0 && return a*sqrt((h - 2q)/h)/2
end;

formula97(5, 4, p=1)

   1.68

formula97(5, 4, q=1)

   1.7677669529663689

formula97(5, 4)

   "error"

4.　図を描いて確かめる

function draw(a, h, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot([0, a, a/2, 0], [0, 0, h, 0], color=:blue, lw=0.5)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   q = 1
   p = formula97(a, h, q=1)
   @printf("p = %g;  q = %g\n", p, q)
   point(a/2, q, @sprintf("長半径 = %g\n短半径 = %g", p, q), :green, :center, delta=-delta/2)
   ellipse(a/2, q, p, q, color=:green)

   p = 1
   q = formula97(a, h, p=1)
   @printf("p = %g;  q = %g\n", p, q)
   point(a/2, q, @sprintf("長半径 = %g\n短半径 = %g", p, q), :red, :center, delta=-delta/2)
   ellipse(a/2, q, p, q, color=:red)
   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;

和算の心（その004）

和算の心（その004）

円に引かれた水平の弦と円弧に囲まれた（狭い方の）隙間に一番大きい円は 1 個ある。二番目以降の円は 2 個ずつある。外円と一番大きい円の直径が与えられたとき，2番目以降の円の直径を求めよ。

「和算の心（その003）」では特別な場合として，弦が円の中心を通る場合について述べた。
今回は，一般の場合について述べる。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive,
     x::positive
y = R - 2r1
eq1 = x2^2 + (y + r2)^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = x/(2R - 2r1) - 2r1/x
eq3 = x^2 + y^2 - R^2
res2 = solve([eq1, eq2, eq3], (r2, x2, x))[1]

   (-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), 2*sqrt(r1)*sqrt(R - r1))

ans_r2 = res2[1](R=>3, r1=> 1/2)
ans_r2 |> println
ans_x2 = res2[2](R=>3, r1=> 1/2)
ans_x2 |> println

   0.416666666666667
   0.52704627669473*sqrt(3)

@syms r3::positive, x3::positive
eq4 = x3^2 + (y + r3)^2 - (R - r3)^2
eq5 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
res3 = solve([eq4, eq5], (r3, x3))[1]
res3 |> println
res3[1](R => 3, r1 => 1/2, r2 => ans_r2, x2 => ans_x2).evalf() |> println
res3[2](R => 3, r1 => 1/2, r2 => ans_r2, x2 => ans_x2).evalf() |> println

   (-(4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2)))/(4*(-R + r1 - r2)), -sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2))
   0.255102040816327
   1.56492159287190

@syms r4::positive, x4::positive
eq6 = x4^2 + (y + r4)^2 - (R - r4)^2
eq7 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq6, eq7], (r4, x4))[1]

   (-(4*R*r1 - 4*r1^2 + x3^2 - 2*x3*(-sqrt(r3)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r3 + x3^2))/(-R + r1 - r3) + x3*(-R + r1)/(-R + r1 - r3)))/(4*(-R + r1 - r3)), -sqrt(r3)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r3 + x3^2))/(-R + r1 - r3) + x3*(-R + r1)/(-R + r1 - r3))

r3, r4 の式を比較すれば，変数名の置き換えだけで次々と半径，と中心座標(x 座標)が求まることがわかる。
そこで，次のような関数を定義しておけば，簡単に計算を続けることができる。

R = 2, r1 = 1/2 のとき，円の半径は以下のようになる（分数は近似値）。

1: 0.5  1//2
2: 0.375  3//8
3: 0.18  9//50
4: 0.0688775510204081  27//392
5: 0.0240928019036288  79//3279
6: 0.00816312819134647  141//17273
7: 0.00273597297804469  2//731
8: 0.000913659014267574  2//2189
9: 0.00030473867949933  1//3281
10: 0.000101600202938299  1//9833

この場合もまた「算額の心（その003）」と同じく，簡単な一般項は求めがたい。

r1 = R/2 のときも，この式が使える。

nextcircle(r2, x2, R, r1) = (-(4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2)))/(4*(-R + r1 - r2)), -sqrt(r2)*sqrt((-R + r1)*(-4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2))/(-R + r1 - r2) + x2*(-R + r1)/(-R + r1 - r2))

function draw(r1, R, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2, x) = (-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), 2*sqrt(r1)*sqrt(R - r1))
   y = R - 2r1
   plot()
   circlef(0, 0, R, :aliceblue)
   circle(0, 0, R, :black)
   circlef(0, R - r1, r1, 1)
   @printf("%2d: %.15g  ", 1, r1)
   println(rationalize(r1, tol=1e-7))
   for i = 2:10
       @printf("%2d: %.15g  ", i, r2)
       println(rationalize(r2, tol=1e-7))
       circle2f(x2, y + r2, r2, i)
       (r2, x2) = nextcircle(r2, x2, R, r1)
   end
   segment(-x, y, x, y)
   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;

和算の心（その003）

和算の心（その003）

円に引かれた水平の弦と円弧に囲まれた（狭い方の）隙間に一番大きい円は 1 個ある。二番目以降の円は 2 個ずつある。一番大きい円の直径を 1 としたとき，2番目以降の円の直径を求めよ。

この問題を解く場合に，算法助法の公式29が役に立つ。

外円の中心を原点とする。
外円の半径を R とすれば，一番大きい円の半径 r1 は r1 = R/2 である。
弦の長さを length とすれば，length^2 = 4(2R)*(2r2) = 16R*r2 であるというのが公式29である。

ここでは特別な場合として，弦が円の中心を通る場合について述べる。この場合，length = 2R である。

include("julia-source.txt")

using SymPy

@syms R, r1, r2, length
formula29 = length^2 - 4(2R)*(2r2)
formula29 = formula29(length => 2R) |> simplify
formula29 |> println

   4*R*(R - 4*r2)

R ≠ 0 なので，R - 4r2 = 0 ということである。r2 について解けば， r2 = R/4 すなわち，2 番目に大きい円の半径は，外円の半径の 1/4 である。
r1 の半径は R/2 だったので，更にその 1/2 である。

公式29を直接使わなくても，以下の連立方程式を解けば，同じ結果が得られる。
このやり方だと，2番目に大きい円の中心座標（x 座標）も同時に得ることができる。また，3番目以降の円の半径を求めることもできる。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive
R = 2r1
eq1 = x2^2 + r2^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
res2 = solve([eq1, eq2], (r2, x2))[1]

   (r1/2, sqrt(2)*r1)

3 番目に大きい円の半径も同様に求めることができる。

@syms r3::positive, x3::positive
r2 = res2[1]
x2 = res2[2]
eq3 = x3^2 + r3^2 - (R - r3)^2
eq4 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
res3 = solve([eq3, eq4], (r3, x3))[1]

   (r1/9, 4*sqrt(2)*r1/3)

4 番目以降の円の半径は前式を初期値として，漸化式で表すことができる。

@syms r4::positive, x4::positive
r3 = res3[1]
x3 = res3[2]
eq5 = x4^2 + r4^2 - (R - r4)^2
eq6 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq5, eq6], (r4, x4))[1]

   (r1/50, 7*sqrt(2)*r1/5)

@syms r5::positive, x5::positive
r4 = res4[1]
x4 = res4[2]
eq7 = x5^2 + r5^2 - (R - r5)^2
eq8 = (x5 - x4)^2 + (r4 - r5)^2 - (r4 + r5)^2
res5 = solve([eq7, eq8], (r5, x5))[1]

   (r1/289, 24*sqrt(2)*r1/17)

@syms r6::positive, x6::positive
r5 = res5[1]
x5 = res5[2]
eq9 = x6^2 + r6^2 - (R - r6)^2
eq10 = (x6 - x5)^2 + (r5 - r6)^2 - (r5 + r6)^2
res6 = solve([eq9, eq10], (r6, x6))[1]

   (r1/1682, 41*sqrt(2)*r1/29)

一番大きい円の半径から順に列挙すると，r1 = R/2 で，
r1, r1/2, r1/9, r1/50, r1/289, r1/1682, ...
となる。一見何らかの規則性がありそうに見えるが，一筋縄ではいかないようだ。いずれにしろ，指数曲線よりは早く減衰するようだ。そのせいかどうか，この場合の累円を求めよという算額はまだ見たことがない。

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2) = (r1/2, sqrt(2)*r1)
   (r3, x3) = (r1/9, 4*sqrt(2)*r1/3)
   (r4, x4) = (r1/50, 7*sqrt(2)*r1/5)
   (r5, x5) = (r1/289, 24*sqrt(2)*r1/17)
   (r6, x6) = (r1/1682, 41*sqrt(2)*r1/29)
   R = 2r1
   @printf("""r1 = %.15g;
       r2 = %.15g;  x2 = %.15g;  
       r3 = %.15g;  x3 = %.15g;
       r4 = %.15g;  x4 = %.15g;
       r5 = %.15g;  x5 = %.15g;
       r6 = %.15g;  x6 = %.15g;\n""",
       r1, r2, x2, r3, x3, r4, x4, r5, x5, r6, x6)
   plot()
   circle(0, 0, R)
   circle22f(0, r1, r1, 1)
   circle4f(x2, r2, r2, 2)
   circle4f(x3, r3, r3, 3)
   circle4f(x4, r4, r4, 4)
   circle4f(x5, r5, r5, 5)
   segment(-R, 0, R, 0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;

算額（その1165）

算額（その1165）

六四 加須市不動岡 総願寺 慶応2年(1866)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，楕円，斜線2本

楕円の中に斜線を 2 本引き，甲円 1 個，乙円 2 個，丙円 2 個を容れる。甲円の直径が 1 寸のとき，丙円の直径はいかほどか。

『埼玉の算額』の図では乙円が楕円と 2 点で接しているように描かれているが，乙円は曲率円であり，楕円の端点と 1 点で接するものである。

楕円の中心を原点とする。
甲円の半径（楕円の短半径）を r1
楕円の長半径を a
乙円の半径と中心座標を r2, (a - r2, 0)
丙円の半径と中心座標を r3, (a - 2r2 - r3, 0)
斜線と楕円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。
(x0, y0) は図を描くためだけに必要なので，r3 を求めるだけなら eq1, eq2, eq3 の三元連立方程式を解くだけでよい。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, a::positive, x0::positive, y0::positive
eq1 = r2 - r1^2/a
eq2 = r2 - (a - r1)/2
eq3 = r3/(r1- r3) - r2/(a - r2)
eq4 = x0^2/a^2 + y0^2/r1^2 - 1
eq5 = dist2(0, 0, x0, y0, a - r2, 0, r2);

res = solve([eq1, eq2, eq3], (r3, r2, a))[1]

   (r1/4, r1/2, 2*r1)

res2 = solve([eq1, eq2, eq3, eq4, eq5], (r3, r2, a, x0, y0))[1]

   (r1/4, r1/2, 2*r1, 2*sqrt(6)*r1/3, sqrt(3)*r1/3)

丙円の半径 r3 は甲円の半径の 1/4 である。ちなみに，乙円の半径 r2 は甲円の半径の 1/2 である。
甲円の直径が 1 寸のとき，丙円の直径は 2 分 5 厘，乙円の直径は 5　分である。

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

  r1 = 0.5;  r3 = 0.125;  r2 = 0.25;  a = 1;; x0 = 0.816497;  y0 = 0.288675

function draw(r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, r2, a, x0, y0) = (r1/4, r1/2, 2*r1, 2*sqrt(6)*r1/3, sqrt(3)*r1/3)
   @printf("甲円の直径が %g のとき，丙円の直径は %g，乙円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  r3 = %g;  r2 = %g;  a = %g;; x0 = %g;  y0 = %g\n", r1, r3, r2, a, x0, y0)
   plot()
   ellipse(0, 0, a, r1, color=:red)
   circle(0, 0, r1, :orange)
   circle2(a - r2, 0, r2, :blue)
   circle2(r1 - r3, 0, r3, :green)
   segment(-x0, -y0, x0, y0, :magenta)
   segment(-x0, y0, x0, -y0, :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, r1, " r1", :red, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :red, :left, :bottom, delta=delta/2)
       point(a - r2, 0, "乙円:r2\n(a-r2,0)", :blue, :center, delta=-2delta)
       point(r1 - r3, 0, "丙円:r3\n(r1-r3,0)", :black, :center, delta=-2delta)
       point(x0, y0, "(x0,y0)  ", :magenta, :right, :vcenter)
   end
end;