算額（その1364）

算額（その1364）

三十四 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円15個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に，大円 5 個，小円 3 個，甲円 2 個，乙円 5 個を容れる。台形の上底が 14 寸，下底が 32 寸のとき，小円の直径はいかほどか。

注：「問」では台形の上底と下底が与えられているが，右側の 3 個の大円の中心を結ぶと正三角形になり，台形の斜辺と下底がなす角は 60° なので，下底が 32 寸のとき，上底は 14 寸にはならない。そこで，下底のみが与えられるとして，問を解く。算額（その1366）も参照のこと。

等脚台形の下底，上底，高さを 2a, 2b, h
大円の半径と中心座標を r1, (2r1, r1), (r1, h - r1)
小円の半径と中心座標を r2, (0, 2r1 + r2)
甲円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (0, h - r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive,
     r1::positive, r2::positive,
     r3::positive, x3::positive, r4::positive;
eq1 = r1^2 + (h - 2r1)^2 - 4r1^2
eq2 = r1/(a - 2r1) - 1/√Sym(3)
eq3 = r1/(b - r1) - √Sym(3);
# res1 = solve([eq1, eq2, eq3], (r1, b, h));

eq4 = r1^2 + (h - r1 - (2r1 + r2))^2 - (r1 + r2)^2
eq5 = r1^2 + (r1 - r4)^2 - (r1 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r4, b, h))

eq6 = r3/(a - x3) - 1/√Sym(3)
eq7 = x3 - 2r1 - 2sqrt(r1*r3)
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7],
   (r1, r2, r3, x3, r4, b, h))[2]  # 2 of 2

   (a*(2 - sqrt(3)), -1129438*a*sqrt(79315912984*sqrt(3) + 137379191137) - 2*a + sqrt(3)*a + 1956244*sqrt(3)*a*sqrt(79315912984*sqrt(3) + 137379191137)/3, a*(-5*sqrt(3) - 4*sqrt(7 - 4*sqrt(3)) + 10)/3, 2*a*(-5*sqrt(3) + 2*sqrt(21 - 12*sqrt(3)) + 9)/3, a*(2 - sqrt(3))/4, a*(3 - sqrt(3))/3, a*(-1694157*sqrt(79315912984*sqrt(3) + 137379191137) - 2*sqrt(3) + 4 + 978122*sqrt(412137573411 + 237947738952*sqrt(3))))

パラメータによっては二重根号を外して簡約化できる。

r1
res[1] |> println

   a*(2 - sqrt(3))

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

   a*(-12 + 7*sqrt(3))/3

r3
res[3] |> sympy.sqrtdenest |> simplify |> println

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

x3
res[4] |> sympy.sqrtdenest |> simplify |> println

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

r4
res[5] |> println

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

b
res[6] |> println

   a*(3 - sqrt(3))/3

h
res[7] |> sympy.sqrtdenest |> simplify |> println

   a

小円の半径は a*(7√3 - 12)/3 なので，下底が 32 寸のとき小円の直径は 32*(7√3 - 12)/3 = 1.3264602984761684 寸である。
このとき，上底は 32*(3 - sqrt(3))/3 = 13.524791385931977 である。

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

   a = 16;  r1 = 4.28719;  r2 = 0.66323;  r3 = 1.42906;  x3 = 13.5248;  r4 = 1.0718;  b = 6.7624;  h = 16

---

なお，小円の直径を知るだけならば，甲円，乙円は不要である（問題を複雑そうに見せる虚仮威しである）。
以下のような手順で小円の直径を（紙と鉛筆だけで）求めることができる。
このことからも，問で与えられた上底の条件は無用なものどころか，有害なものであることがわかる。

下底 2a と大円の半径 r1 の関係式は以下であり，大円の半径は r1 = a*(2 - √3) である。

@syms a, r1
eq01 = (2 + √Sym(3))r1 - a
eq01 |> println

   -a + r1*(sqrt(3) + 2)

r1 = solve(eq01, r1)[1] |> simplify
r1 |> println

   a*(2 - sqrt(3))

一辺の長さが 2r1 の正三角形の頂点を中心とする3個の大円が作る中央の隙間にピッタリハマる円の半径は以下のようになる（図中の灰色で描いた小さな直角三角形を参照）。

r1*2/√Sym(3) - r1 |> simplify |> println

   a*(-12 + 7*sqrt(3))/3

下底が 32 寸のとき，小円の直径は 1.3264602984761684 寸である。

a = 32/2
2*a*(-12 + 7*sqrt(3))/3

   1.3264602984761684

以下は図を描くプログラム

function draw(a, more)
   pyplot(size=(800, 800), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x3, r4, b, h) = (a*(2 - sqrt(3)), -1129438*a*sqrt(79315912984*sqrt(3) + 137379191137) - 2*a + sqrt(3)*a + 1956244*sqrt(3)*a*sqrt(79315912984*sqrt(3) + 137379191137)/3, a*(-5*sqrt(3) - 4*sqrt(7 - 4*sqrt(3)) + 10)/3, 2*a*(-5*sqrt(3) + 2*sqrt(21 - 12*sqrt(3)) + 9)/3, a*(2 - sqrt(3))/4, a*(3 - sqrt(3))/3, a*(-1694157*sqrt(79315912984*sqrt(3) + 137379191137) - 2*sqrt(3) + 4 + 978122*sqrt(412137573411 + 237947738952*sqrt(3))))
   r1 = a*(2 - √3)
   r2 = a*(7√3 - 12)/3
   r3 = a*(2 - √3)/3
   x3 = 2*a*(3 - √3)/3
   r4 = a*(2 - √3)/4
   b = a*(3 - √3)/3
   h = a
   @printf("下底が %g のとき，小円の直径は %g である（上底は %g）。\n", 2a, 2r2, 2b)
   @printf("a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  b = %g;  h = %g\n", a, r1, r2, r3, x3, r4, b, h)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:brown, lw=0.5)
   circle(0, r1, r1)
   circle2(2r1, r1, r1)
   circle2(r1, h - r1, r1)

   circle(0, 2r1 + r2, r2, :blue)
   circle2(r1, h - 2r1 - r2, r2, :blue)

   circle(0, h - r4, r4, :green)
   circle2(r1, r4, r4, :green)

   circle2(x3, r3, r3, :orange)
   l = h - (2r1 + r2) - r4
   circle2(r1 + l*cosd(30), h - (2r1 + r2) + l*sind(30), 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(2r1, r1, "大円:r1,(2r1,r1)", :red, :center, delta=-delta)
       point(0, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta)
       point(r1, (1 + √3)r1, "大円:r1,(r1,(1+√3)r1)", :red, :center, delta=-delta)
       point(0, 2r1 + r2, "小円:r2,(0,2r1+r2)", :blue, :left, :vcenter, deltax=4delta)
       point(0, h - r4, "乙円:(0,h-r4)", :green, :left, :vcenter, deltax=7delta)
       point(x3, r3, "甲円:(x3,r3)", :orange, :right, :vcenter, deltax=-9delta)
       point(b, h, "(b,h)", :brown, :right, :bottom, delta=delta/2)
       point(a, 0, "a", :brown, :left, :bottom, delta=delta/2)
       plot!([r1, 2r1, r1, r1], [r1, r1, h - 2r1 - r2, r1], color=:gray70, lw=0.5)
   end  
end;

draw(32/2, true)

 

算額（その1363）

算額（その1363）

三十三 群馬県佐波郡境町伊与久 雷電神社 文化14年(1817)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円14個，外円，正三角形
#Julia, #SymPy, #算額, #和算

外円の中に正三角形と円 13 個を容れる。甲円の直径が 39.9 寸のとき，丙円の直径はいかほどか。

図形を反時計回りに 90° 回転させたものを考える。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 0)
乙円の半径と中心座標を r2, (r1 + r2, r2)
丙円の半径と中心座標を r3, (r1 + r3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r3::positive, y3::positive;
R = 2r1
eq1 = (r1 + r2)^2 + r2^2 - (R - r2)^2
eq2 = (r1 + r3)^2 + y3^2 - (R - r3)^2
eq3 = y3 - r2 - 2sqrt(r2*r3)
res = solve([eq1, eq2, eq3], (r2, r3, y3))[1];

丙円の半径は r1*(-267 - 56*sqrt(21 - 12*sqrt(3)) + 36*sqrt(7 - 4*sqrt(3)) + 190*sqrt(3))/169 となるが，これは二重根号を解消することにより，実にスッキリと簡約化される。

res[2] |> println

   r1*(-267 - 56*sqrt(21 - 12*sqrt(3)) + 36*sqrt(7 - 4*sqrt(3)) + 190*sqrt(3))/169

丙円の半径は，甲円の半径の 3(14√3 - 9)/169 倍である。

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

   3*r1*(-9 + 14*sqrt(3))/169

甲円の直径が 39.9 寸のとき，丙円の直径は 39.9 * 3(14√3 - 9)/169 = 10.800418599549849 寸である。

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

  R = 39.9;  r1 = 19.95;  r2 = 9.25883;  r3 = 5.40021;  y3 = 23.4009

function draw(r1, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, y3) = (r1*(-3 + 2*sqrt(3)), r1*(-267 - 56*sqrt(21 - 12*sqrt(3)) + 36*sqrt(7 - 4*sqrt(3)) + 190*sqrt(3))/169, r1*(-6*sqrt(3) + 12*sqrt(7 - 4*sqrt(3)) + 16*sqrt(21 - 12*sqrt(3)) + 15)/13)
   R = 2r1
   @printf("甲円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r3)
   @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  y3 = %g\n", R, r1, r2, r3, y3)
   plot([r1, -R, r1, r1], [√3r1, 0, -√3r1, √3r1], color=:brown, lw=0.5)
   circle(0, 0, R, :orange)
   circle(0, 0, r1, :magenta)
   rotate(r1 + r2, r2, r2, :blue)
   rotate(r1 + r2, -r2, r2, :blue)
   rotate(r1 + r3, y3, r3, :green)
   rotate(r1 + r3, -y3, 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(r1 + r2, r2, "乙円:r2\n(r1+r2,r2)", :blue, :center, delta=-delta)
       point(r1 + r3, y3, "丙円:r3,(r1+r3,y3) ", :green, :right, :vcenter)
       point(r1, 0, "r1 ", :magenta, :right, :bottom, delta=delta/2)
       point(R, 0, " R", :orange, :left, :bottom, delta=delta/2)
       point(r1, √3r1, "(r1,√3r1)", :brown, :left, :vcenter)
   end  
end;

draw(39.9/2, true)

算額（その1362）

算額（その1362）

三十一 群馬県藤岡市立石 立石神社 文化13年(1816)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に全円 1 個，等円 4 個を容れる。全円の直径が与えられたとき，等円の直径を求める術をのべよ。

等脚台形の下底，上底，高さを 2a, 2b, h
全円の半径と中心座標を r1, (0, r1)
等円の半径と中心座標を r2, (x2, r2), (r2, h - r2)
とおき，以下の連立方程式の数値解を求める。
SymPy では，性能的に解析解を求めることができない。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive,
     r1::positive, r2::positive, x2::positive;
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = r2^2 + (h - r2 - r1)^2 - (r1 + r2)^2
eq3 = dist2(b, h, a, 0, x2, r2, r2)
eq4 = dist2(b, h, a, 0, r2, h - r2, r2)
eq5 = dist2(b, h, a, 0, 0, r1, r1);

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)
   (a, b, h, r2, x2) = u
   return [
       x2^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       r2^2 - (r1 + r2)^2 + (h - r1 - r2)^2,  # eq2
       h*(a^2*h - 2*a^2*r2 + 2*a*b*r2 - 2*a*h*x2 + 2*a*r2*x2 - 2*b*r2*x2 - h*r2^2 + h*x2^2),  # eq3
       h*(2*a*b*r2 - 2*a*r2^2 + b^2*h - 2*b^2*r2 - 2*b*h*r2 + 2*b*r2^2),  # eq4
       h*(a^2*h - 2*a^2*r1 + 2*a*b*r1 - h*r1^2),  # eq5
   ]
end;

r1 = 1
iniv = BigFloat[1.5, 0.5, 2.5, 0.28, 1.1]
res = nls(H, ini=iniv)

   ([1.4706834198711607, 0.4706834198711606, 2.5293165801288393, 0.28017604199929, 1.0586331602576788], true)

全円の直径が 2 のとき，等円の直径は 0.56035208399858 である。

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

    r1 = 1;  a = 1.47068;  b = 0.470683;  h = 2.52932;  r2 = 0.280176;  x2 = 1.05863

function draw(r1, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, h, r2, x2) = res[1]
   @printf("全円の直径が %g のとき，等円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  a = %g;  b = %g;  h = %g;  r2 = %g;  x2 = %g\n", r1, a, b, h, r2, x2)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:green, lw=0.5)
   circle(0, r1, r1, :blue)
   circle2(x2, r2, r2)
   circle2(r2, h - r2, r2)
   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,(0,r1)", :blue, :center, delta=-delta/2)
       point(x2, r2, "等円:r2\n(x2,r2)", :red, :center, delta=-delta/2)
       point(r2, h - r2, "等円:r2\n(r2,h-r2)", :red, :center, :bottom, delta=delta/2)
       point(b, h, "(b,h)", :green, :right, :bottom, delta=delta/2)
       point(a, 0, "a", :green, :left, :bottom, delta=delta/2)
   end  
end;

draw(1, true)

算額（その1361）

算額（その1361）

二十六 群馬県安中市簗瀬 稲荷神社 文化11年(1814)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：正方形2個，長方形，斜線2本
#Julia, #SymPy, #算額, #和算

長方形の中に，斜線 2 本と，大小の正方形 2 個を容れる。2 つの正方形の一辺の長さの積が 15 平方寸，和が 8 寸のとき，長方形の面積が最小になるときの長方形の短辺はいかほどか。

長方形の長辺と短辺の長さを a, b
大きな正方形と小さな正方形の一辺の長さをそれぞれ d, e; d ≦ e
2本の斜線の長方形の長辺上での交点座標を (c, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     d::positive, e::positive;

まず，「2 つの正方形の一辺の長さの積が 15 平方寸，和が 8 寸」ということから，正方形の一辺の長さを求める。

eq1 = d*e - 15
eq2 = d + e - 8
res1 = solve([eq1, eq2], (d, e))[1]  # 1 of 2

   (3, 5)

長方形の一辺の長さは d = 3 寸と e = 5 寸である。

これが既知として，残りの条件を列挙するが，独立なのは 2 つだけである。

(d, e) = (3, 5)
eq3 = eq1 = d/(c - d) - b/c
eq4 = (b - d)/d - b/c
eq5 = e/(a - e - c) - b/(a - c)
eq6 = (b - e)/e - b/(a - c);

eq3, eq5 を連立させて，a, b を求める。

res2 = solve([eq3, eq5], (a, b))[1]

   (2*c*(c - 15)/(2*c - 15), 3*c/(c - 3))

a, b は c の関数として表すことができる。
a, b に基づき，長方形の面積 a * b を計算する。

S = res2[1] * res2[2]
S |> println

   6*c^2*(c - 15)/((c - 3)*(2*c - 15))

S のグラフを描いてみると，（少なくとも）c が 3 〜 6 の間で最小値を取りそうであるということがわかる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(S,  xlims=(3.1, 7.4), xlabel="c", ylabel="面積")

面積が最小になるときの c は，面積を表す関数の導関数を取り，導関数 = 0 の方程式を解けばよい。

導関数
diff_S = diff(S, c) |> simplify
diff_S |> println

   12*c*(c^3 - 21*c^2 + 225*c - 675)/(4*c^4 - 84*c^3 + 621*c^2 - 1890*c + 2025)

res3 = solve(diff_S, c)[3];  # 3 of 3
res3 |> println

   -(107 + 15*sqrt(129))^(1/3) + 26/(107 + 15*sqrt(129))^(1/3) + 7

ans_c = res3.evalf()
ans_c |> println

   4.46521053650337

c = 26/(107 + 15√129)^(1/3) - (107 + 15√129)^(1/3) + 7 = 4.46521053650337 のとき，面積は最小になる。

c = 4.46521053650337 のとき，a, b は 15.5002689570601, 9.14246197101333 である。
すなわち，長方形の長辺は 15.5002689570601，短辺は 9.14246197101333 である。

また，面積は 15.5002689570601 * 9.14246197101333 = 141.710619480400 である。

(res2[1](c => ans_c), res2[2](c => ans_c)) |> println
(res2[1] * res2[2])(c => ans_c) |> println

   (15.5002689570601, 9.14246197101333)
   141.710619480400

function draw(K1, K2, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = 4.46521053650337
   (d, e) = [3, 5]
   (a, b) = (2*c*(c - 15)/(2*c - 15), 3*c/(c - 3))
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
   plot!([0, c, a], [b, 0, b], color=:blue, lw=0.5)
   rect(0, 0, d, d, :red)
   rect(a - e, 0, a, e, :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)
       point(a, b, "(a,b)", :green, :right, :bottom, delta=delta/2)
       point(d, d, " (d,d)", :red, :left, :vcenter)
       point(a - e, e, "(a-e,e)  ", :red, :right, :vcenter)
       point(c, 0, "c", :blue, :center, delta=-delta/2)
       ylims!(-5delta, b + 5delta)
   end  
end;

draw(15, 8, true)

算額（その1360）

算額（その1360）

二十六 群馬県安中市簗瀬 稲荷神社 文化11年(1814)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，正方形，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に正方形 1 個と等円 4 個を容れる。等円の直径が 3 寸のとき，正方形の一辺の長さはいかほどか。



等脚台形の下底，上底を 2a, 2b，高さを h
等円の半径と中心座標を r, (x, r), (r, h - r)
等脚台形の斜辺上の正方形の頂点座標を (c, c)
とおき，以下の連立方程式を解く。
正方形の一辺の長さは √2c である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, h::positive,
     r::positive, x::positive
eq1 = dist2(a, 0, b, h, x, r, r)
eq2 = dist2(a, 0, b, h, r, h - r, r)
eq3 = dist2(0, 0, c, c, x, r, r)
eq4 = dist2(0, 2c, c, c, r, h - r, r)
eq5 = c/(a - c) - h/(a - 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 Float64.(v), r.f_converged
end;

function H(u)
   (a, b, c, h, x) = u
   return [
       h*(a^2*h - 2*a^2*r + 2*a*b*r - 2*a*h*x + 2*a*r*x - 2*b*r*x - h*r^2 + h*x^2),  # eq1
       h*(2*a*b*r - 2*a*r^2 + b^2*h - 2*b^2*r - 2*b*h*r + 2*b*r^2),  # eq2
       -r^2/2 - r*x + x^2/2,  # eq3
       2*c^2 - 2*c*h + h^2/2 - r^2,  # eq4
       c/(a - c) - h/(a - b),  # eq5
   ]
end;

r = 3/2
iniv = BigFloat[36.8/2, 18.7/2, 29.1/2, 40, 12]
res = nls(H, ini=iniv)

   ([5.576158402821451, 2.6509904819684404, 4.397777478867205, 10.916875301294052, 3.6213203435596424], true)

√2*4.397777478867205

   6.219396554912959

等円の直径が 3 寸のとき，正方形の一辺の長さは 6.219396554912959 である。

「答」では 6 寸としているが，「術」は 等円の直径の (√0.5 + 0.5 + √0.75) 倍としており，正しく 3(√0.5 + 0.5 + √0.75) = 6.219396554912959 となる。

function draw(r, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, h, x) = [36.8/2, 18.7/2, 29.1/2, 36, 12]
   (a, b, c, h, x) = res[1]
   length = √2c
   @printf("等円の直径が %g のとき，正方形の一辺の長さは %g である。\n", 2r, length)
   @printf("r = %g;  a = %g;  b = %g;  c = %g;  h = %g;  x = %g\n", r, a, b, c, h, x)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:green, lw=0.5)
   plot!([0, c, 0, -c, 0], [0, c, 2c, c, 0], color=:magenta, lw=0.5)
   circle2(x, r, r)
   circle2(r, h - 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, 2c, "(0,2c)", :magenta, :left, :vcenter, deltax=delta)
       point(c, c, "(c,c)", :magenta, :left, :vcenter, deltax=delta)
       point(a, 0, "a", :green, :left, :bottom, delta=delta/2, deltax=delta/3)
       point(b, h, "(b,h)", :green, :right, :bottom, delta=delta/2)
       point(x, r, "等円:r,(x,r)", :red, :center, delta=-delta/2)
       point(r, h - r, "等円:r,(r,h-r)", :red, :center, :bottom, delta=delta/2)
   end  
end;

draw(3/2, true)

算額（その1359）

算額（その1359）

二十 群馬県群馬郡群馬町 旧福寿庵 文化8年(1811)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，長方形，斜線2本
#Julia, #SymPy, #算額, #和算

長方形の中に斜線を 2 本引き，上円 1 個，下円 2 個を容れる。（長方形の）縦（長辺）横（短辺）の積が 336 平方寸，上円と下円の直径の積が 63 平方寸のとき，上円の直径はいかほどか。

縦横の積を K1，上円と下円の直径の積を K2
横を 2a，縦を b
上円の半径と中心座標を r1, (0, b - r1)
下円の半径と中心座標を r2, (a - r2, r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive,
     K1::positive, K2::positive
eq1 = dist2(0, 0, a, b, 0, b - r1, r1)
eq2 = dist2(0, 0, a, b, a - r2, r2, r2)
eq2 = a + b - sqrt(a^2 + b^2) - 2r2
eq2 = (a + b + sqrt(a^2 + b^2))*r2 - a*b
eq3 = 2a*b - K1
eq4 = 2r1*2r2 - K2
res = solve([eq1, eq2, eq3, eq4], (a, b, r1, r2));  # 3 of 4

a, b, r1, r2 は K1, K2 を含む長い式になり，　SymPy では簡約化できないが，実際の値を代入すると数値が求まる。

a = 横/2
res[3][1] |> println
2res[3][1](K1 => 336, K2 => 63).evalf() |> println

   sqrt(2)*sqrt(-sqrt(K1)*sqrt(K1 - 4*K2) + 3*K1)*(-48*K1^(11/2)*sqrt(K1 - 4*K2) - 40*K1^(9/2)*K2*sqrt(K1 - 4*K2) + 4*K1^(7/2)*K2^2*sqrt(K1 - 4*K2) + 6*K1^(5/2)*K2^3*sqrt(K1 - 4*K2) - 2*K1^(3/2)*K2^4*sqrt(K1 - 4*K2) - sqrt(K1)*K2^5*sqrt(K1 - 4*K2) + 48*K1^6 + 136*K1^5*K2 + 76*K1^4*K2^2 - 14*K1^3*K2^3 - 10*K1^2*K2^4 + 5*K1*K2^5 + 2*K2^6)/(8*sqrt(K2)*sqrt(2*K1 + K2)*(48*K1^5 + 40*K1^4*K2 - 4*K1^3*K2^2 - 6*K1^2*K2^3 + 2*K1*K2^4 + K2^5))
   14.0000000000000

b　= 縦
res[3][2] |> println
2res[3][2](K1 => 336, K2 => 63).evalf() |> println

   sqrt(-2*sqrt(K1)*sqrt(K1 - 4*K2) + 6*K1)*(12*K1^(11/2)*sqrt(K1 - 4*K2) + 10*K1^(9/2)*K2*sqrt(K1 - 4*K2) - K1^(7/2)*K2^2*sqrt(K1 - 4*K2) - 3*K1^(5/2)*K2^3*sqrt(K1 - 4*K2)/2 + K1^(3/2)*K2^4*sqrt(K1 - 4*K2)/2 + sqrt(K1)*K2^5*sqrt(K1 - 4*K2)/4 + 12*K1^6 + 10*K1^5*K2 - K1^4*K2^2 - 3*K1^3*K2^3/2 + K1^2*K2^4/2 + K1*K2^5/4)/(sqrt(K2)*sqrt(2*K1 + K2)*(48*K1^5 + 40*K1^4*K2 - 4*K1^3*K2^2 - 6*K1^2*K2^3 + 2*K1*K2^4 + K2^5))
   48.0000000000000

r1 = 上円の半径
res[3][3] |> println
2res[3][3](K1 => 336, K2 => 63).evalf() |> println

   sqrt(2)*sqrt(K2)*sqrt(-sqrt(K1)*sqrt(K1 - 4*K2) + 3*K1)*(8*K1^(11/2)*sqrt(K1 - 4*K2) + 12*K1^(9/2)*K2*sqrt(K1 - 4*K2) + 6*K1^(7/2)*K2^2*sqrt(K1 - 4*K2) + K1^(5/2)*K2^3*sqrt(K1 - 4*K2) + 24*K1^6 + 36*K1^5*K2 + 18*K1^4*K2^2 + 3*K1^3*K2^3)/(8*K1^3*sqrt(2*K1 + K2)*(8*K1^3 + 12*K1^2*K2 + 6*K1*K2^2 + K2^3))
   10.5000000000000

r2 = 下円の半径
res[3][4] |> println
2res[3][4](K1 => 336, K2 => 63).evalf() |> println

   sqrt(2)*sqrt(K2)*sqrt(-sqrt(K1)*sqrt(K1 - 4*K2) + 3*K1)/(4*sqrt(2*K1 + K2))
   6.00000000000000

縦横の積が 336，上円下円の直径の積が 63 のとき，横は 14，縦は 24，大円の直径は 10.5，小円の直径は 6 である。

function draw(K1, K2, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r1, r2) = (sqrt(2)*sqrt(-sqrt(K1)*sqrt(K1 - 4*K2) + 3*K1)*(-48*K1^(11/2)*sqrt(K1 - 4*K2) - 40*K1^(9/2)*K2*sqrt(K1 - 4*K2) + 4*K1^(7/2)*K2^2*sqrt(K1 - 4*K2) + 6*K1^(5/2)*K2^3*sqrt(K1 - 4*K2) - 2*K1^(3/2)*K2^4*sqrt(K1 - 4*K2) - sqrt(K1)*K2^5*sqrt(K1 - 4*K2) + 48*K1^6 + 136*K1^5*K2 + 76*K1^4*K2^2 - 14*K1^3*K2^3 - 10*K1^2*K2^4 + 5*K1*K2^5 + 2*K2^6)/(8*sqrt(K2)*sqrt(2*K1 + K2)*(48*K1^5 + 40*K1^4*K2 - 4*K1^3*K2^2 - 6*K1^2*K2^3 + 2*K1*K2^4 + K2^5)), sqrt(-2*sqrt(K1)*sqrt(K1 - 4*K2) + 6*K1)*(12*K1^(11/2)*sqrt(K1 - 4*K2) + 10*K1^(9/2)*K2*sqrt(K1 - 4*K2) - K1^(7/2)*K2^2*sqrt(K1 - 4*K2) - 3*K1^(5/2)*K2^3*sqrt(K1 - 4*K2)/2 + K1^(3/2)*K2^4*sqrt(K1 - 4*K2)/2 + sqrt(K1)*K2^5*sqrt(K1 - 4*K2)/4 + 12*K1^6 + 10*K1^5*K2 - K1^4*K2^2 - 3*K1^3*K2^3/2 + K1^2*K2^4/2 + K1*K2^5/4)/(sqrt(K2)*sqrt(2*K1 + K2)*(48*K1^5 + 40*K1^4*K2 - 4*K1^3*K2^2 - 6*K1^2*K2^3 + 2*K1*K2^4 + K2^5)), sqrt(2)*sqrt(K2)*sqrt(-sqrt(K1)*sqrt(K1 - 4*K2) + 3*K1)*(8*K1^(11/2)*sqrt(K1 - 4*K2) + 12*K1^(9/2)*K2*sqrt(K1 - 4*K2) + 6*K1^(7/2)*K2^2*sqrt(K1 - 4*K2) + K1^(5/2)*K2^3*sqrt(K1 - 4*K2) + 24*K1^6 + 36*K1^5*K2 + 18*K1^4*K2^2 + 3*K1^3*K2^3)/(8*K1^3*sqrt(2*K1 + K2)*(8*K1^3 + 12*K1^2*K2 + 6*K1*K2^2 + K2^3)), sqrt(2)*sqrt(K2)*sqrt(-sqrt(K1)*sqrt(K1 - 4*K2) + 3*K1)/(4*sqrt(2*K1 + K2)))
   @printf("縦横の積が %g，上円下円の直径の積が %g のとき，横は %g，縦は %g，大円の直径は %g，小円の直径は %g である。\n", K1, K2, 2a, b, 2r1, 2r2)
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:green, lw=0.5)
   circle(0, b - r1, r1)
   circle(a - r2, r2, r2, :blue)
   circle(r2 - a, r2, r2, :blue)
   plot!([-a, 0, a], [b, 0, b], color=:magenta, 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, b - r1, "上円:r1,(0,b-r1)", :red, :center, delta=-delta/2)
       point(a - r2, r2, "下円:re2\n(a-r2,r2)", :blue, :center, delta=-delta/2)
       point(0, b, "b", :green, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
   end  
end;

draw(336, 63, true)

算額（その1358）

算額（その1358）

十九 群馬県群馬郡榛名町榛名山 榛名神社 文化8年(1811)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，長方形，斜線
#Julia, #SymPy, #算額, #和算

長方形の中に斜線を隔てて大円 1 個，中円 2 個，小円 1 個を容れる。大円の直径が 10 寸のとき，小円の直径はいかほどか。

長方形の短辺と長辺の長さを a, b
斜線と長方形の辺の交点座標を (c, 0), (a, d)
大円の半径と中心座標を r1, (r1, b - r1)
中円の半径と中心座標を r2, (r2, r2), (a - r2, r2)
小円の半径と中心座標を r3, (a - r3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, d::positive,
     r1::positive, r2::positive, r3::positive, y3::positive
eq1 = dist2(c, 0, a, d, r1, b - r1, r1)
eq2 = dist2(c, 0, a, d, r2, r2, r2)
eq3 = dist2(c, 0, a, d, a - r2, r2, r2)
eq4 = dist2(c, 0, a, d, a - r3, y3, r3)
eq5 = (r1 - r2)^2 + (b - r1 - r2)^2 - (r1 + r2)^2
eq6 = (r2 - r3)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq7 = (a - c) + d - sqrt((a - c)^2 + d^2) - 2r2;
# eq8 = ((a - c) + d + sqrt((a - c)^2 + d^2))*r2 - (a - c)*d;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (a, b, c, d, r2, r3, y3))

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)
   (a, b, c, d, r2, r3, y3) = u
   return [
       a^2*b^2 - 2*a^2*b*r1 - 2*a*b^2*c + 2*a*b*c*d + 4*a*b*c*r1 - 2*a*b*d*r1 - 2*a*c*d*r1 + 2*a*d*r1^2 + b^2*c^2 - 2*b*c^2*d - 2*b*c^2*r1 + 2*b*c*d*r1 + c^2*d^2 + 2*c^2*d*r1 - 2*c*d^2*r1 - 2*c*d*r1^2,  # eq1
       d*(2*a*c*r2 - 2*a*r2^2 + c^2*d - 2*c^2*r2 - 2*c*d*r2 + 2*c*r2^2),  # eq2
       d*(a^2*d - 2*a^2*r2 - 2*a*c*d + 4*a*c*r2 - 2*a*d*r2 + 2*a*r2^2 + c^2*d - 2*c^2*r2 + 2*c*d*r2 - 2*c*r2^2),  # eq3
       a^2*d^2 - 2*a^2*d*y3 - a^2*r3^2 + a^2*y3^2 - 2*a*c*d^2 + 4*a*c*d*y3 + 2*a*c*r3^2 - 2*a*c*y3^2 - 2*a*d^2*r3 + 2*a*d*r3*y3 + c^2*d^2 - 2*c^2*d*y3 - c^2*r3^2 + c^2*y3^2 + 2*c*d^2*r3 - 2*c*d*r3*y3,  # eq4
       (r1 - r2)^2 - (r1 + r2)^2 + (b - r1 - r2)^2,  # eq5
       (-r2 + y3)^2 + (r2 - r3)^2 - (r2 + r3)^2,  # eq6
       a - c + d - 2*r2 - sqrt(d^2 + (a - c)^2),  # eq7
   ]
end;
r1 = 10/2
#iniv = BigFloat[13.2, 16, 5, 16, 3, 2, 8]
iniv = BigFloat[12.2, 15.2, 4.3, 12.2, 2.8, 1.5, 6.9]
res = nls(H, ini=iniv)

   ([12.21558138691064, 15.247885379719106, 4.29489761141387, 12.16723318645839, 2.784847053174034, 1.5510746219144202, 6.941531111454818], true)

数値解は初期値の与え方で，見た目の図が変わる。

例えば，先に掲げた図は，以下のパラメータで描画したものである。

大円の直径が 10 のとき，小円の直径は 3.10215 である。中円の直径は 5.56969 である。
r1 = 5;  a = 12.2156;  b = 15.2479;  c = 4.2949;  d = 12.1672;  r2 = 2.78485;  r3 = 1.55107;  y3 = 6.94153

答曰では，小円の直径を 3.8497有奇としているが，その解に近い図は下のようになる。

パラメータは，以下のとおりである。
大円の直径が 10 のとき，小円の直径は 3.86769 である。中円の直径は 6.21908 である。
r1 = 5;  a = 13.1965;  b = 15.9957;  c = 5.01652;  d = 16.0811;  r2 = 3.10954;  r3 = 1.93385;  y3 = 8.01398

解が一意にならない原因は，この図形の条件として，大円の直径一つだけでは解は無数にあるということであろう。制約条件がもう一つ必要である。

function draw(r1, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, d, r2, r3, y3) = res[1]
   @printf("大円の直径が %g のとき，小円の直径は %g である。中円の直径は %g である。\n", 2r1, 2r3, 2r2)
   @printf("r1 = %g;  a = %g;  b = %g;  c = %g;  d = %g;  r2 = %g;  r3 = %g;  y3 = %g\n", r1, a, b, c, d, r2, r3, y3)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
   circle(r1, b - r1, r1)
   circle(r2, r2, r2, :blue)
   circle(a - r2, r2, r2, :blue)
   circle(a - r3, y3, r3, :magenta)
   segment(c, 0, a, d, :orange)
   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, b - r1, "大円:r1,(r1,b-r1)", :red, :center, delta=-delta)
       point(r2, r2, "中円:r2,(r2,r2)", :blue, :center, delta=-delta)
       point(a - r2, r2, "中円:r2,(a-r2,r2)", :blue, :center, delta=-delta)
       point(a - r3, y3, "小円:r3\n(a-r3,y3)", :magenta, :center, :bottom, delta=delta)
       point(a, b, "(a,b)", :green, :right, :bottom, delta=delta/2)
       point(a, d, "(a,d)", :green, :right, :bottom, delta=delta/2)
       point(c, 0, "(c,0)", :green, :left, :bottom, delta=delta/2, deltax=2delta)
   end  
end;

draw(10/2, true)

算額（その1357）

算額（その1357）

十八 群馬県高崎市八幡町 八幡宮 文化7年(1810)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，正方形
#Julia, #SymPy, #算額, #和算

正方形の中に，大円 2 個，小円 3 個を容れる。正方形の一辺の長さが 5 寸のとき，小円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive
eq1 = 2(r1 - a)^2 - (r1 - r2)^2
eq2 = (a - r1 - (a - r2))^2 + (r2 - a - (a - r1))^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r2, r1))[4];  # 4 of 4

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

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

小円の半径は，a*(2sqrt(√2 + 2) - √2)*(2sqrt(√2 + 2) - 2 - √2)/2
a = 5 のとき，小円の直径は 1.60435348784316

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

   a*(-2*sqrt(sqrt(2) + 2) + sqrt(2) + 3)

大円の半径は，a*(√2 + 3 - 2sqrt(√2 + 2))
a = 5 のとき，大円の直径は 3.59347716163974

正方形の一辺の長さが 5 のとき，小円の直径は 1.60435348784316，大円の直径は 3.59347716163974 である。

2*res[1](a => 5/2).evalf() |> println
2*res[2](a => 5/2).evalf() |> println

   1.60435348784316
   3.59347716163974

function draw(a, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (-(-10*a^3 - 20*a^3*(-sqrt(28*sqrt(2) + 40) + 3*sqrt(2) + 5)^2 + a^3*(-sqrt(28*sqrt(2) + 40) + 3*sqrt(2) + 5)^3 + 33*a^3*(-sqrt(28*sqrt(2) + 40) + 3*sqrt(2) + 5))/(2*a^2), a*(-sqrt(28*sqrt(2) + 40) + 3*sqrt(2) + 5))
   @printf("正方形の一辺の長さが %g のとき，小円の直径は %g，大円の直径は %g である。\n", 2a, 2r2, 2r1)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:green, lw=0.5)
   circle(a - r1, a - r1, r1)
   circle(r1 - a, r1 - a, r1)
   circle(0, 0, r2, :blue)
   circle(a - r2, r2 - a, r2, :blue)
   circle(r2 - a, a - 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 - r1, a - r1, "大円:r1,(a-r1,a-r1)", :red, :left, :bottom, delta=delta, deltax=-4delta)
       point(a - r2, r2 - a, "小円:r2,(a-r2,r2-a)", :blue, :center, delta=-delta/2)
       point(a, a, "(a,a)", :green, :right, :bottom, delta=delta)
   end  
end;

draw(5/2, true)

 

算額（その1356）

算額（その1356）

十七 群馬県高崎市八幡町 八幡宮 文化7年(1810)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形（梯形）の中に大円 1 個，小円 2 個を容れる。大円の直径が 7 寸，大頭（下底）が  18 寸と 2/3 のとき，小円の直径はいかほどか。

等脚台形の下底，上底を 2a, 2b，高さを h
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (r2, h - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

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

SymPy の性能上，一度には解けないので，eq1, eq2 を解いて，h, b を求め，その解を eq3 に代入して r2 について解く。

res = solve([eq1, eq2], (h, b))[1]  # 1 of 2

   (2*r1*(a^2 - a*r2 - r1*r2)/(a^2 - r1^2), r2*(a + r1)/a)

eq13 = eq3(h => res[1])
eq13 |> println

   r2^2 - (r1 + r2)^2 + (-r1 + 2*r1*(a^2 - a*r2 - r1*r2)/(a^2 - r1^2) - r2)^2

ans_r2 = solve(eq13, r2)[1]  # 1 of 2
ans_r2 |> println

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

小円の直径は，大円の直径と下底の関数として得られる。

r1 = 14/2, a = (18 +2/3)/2 を代入して 2 倍すれば，小円の直径 32/7 = 4 + 4/7 が得られる。

2ans_r2(r1 => 14/2, a => (18 + 2/3)/2) |> x -> rationalize(float(x), tol=1e-3)

   32//7

大円の直径が 14，下底が 18.6667(56/3) のとき，小円の直径は 4.57143(32/7) である。
上底は 8, 高さは 18.2857(128/7) である。

str(x) = replace(string(rationalize(x, tol=1e-3)), "//" => "/")

function draw(r1, a, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2r1*(a^2 - a*r1 - a*sqrt(a^2 + r1^2) + r1^2 + r1*sqrt(a^2 + r1^2))/(a^2 + 2a*r1 + r1^2)
   (h, b) = (2r1*(a^2 - a*r2 - r1*r2)/(a^2 - r1^2), r2*(a + r1)/a)
   str1 = str(2a)
   str2 = str(2r2)
   str3 = str(h)
   @printf("大円の直径が %g，下底が %g(%s) のとき，小円の直径は %g(%s) である。\n上底は %g, 高さは %g(%s) である。\n", 2r1, 2a, str1, 2r2, str2, 2b, h, str3)
   println()
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:blue, lw=0.5)
   circle2(r2, h - r2, r2)
   circle(0, r1, r1, :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, 0, "a", :blue, :left, :bottom, delta=delta/2)
       point(b, h, "(b,h)", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, "大円:r1,(0,r1)", :green, :center, delta=-delta/2)
       point(r2, h - r2, "小円:r2\n(r2,h-r2)", :red, :center, delta=-delta/2)
   end  
end;

draw(14/2, (18 + 2/3)/2, true)

 

算額（その1355）

算額（その1355）

十三 群馬県前橋市元総社町 総社神社 享和年中(1801〜1804)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個，三角形
#Julia, #SymPy, #算額, #和算

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

三角形の底辺の長さを a，頂点の座標を (b, c)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x31, r3), (x32, r3)
丁円の半径と中心座標を r4, (x41, r4), (x42, y42)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, x1::positive,
     r2::positive, x2::positive,
     r3::positive, x31::positive, x32::positive,
     r4::positive, x41::positive, x42::positive, y42::positive
eq1 = x2 - x1 - 2sqrt(r1*r2)
eq2 = x32 - x2 - 2sqrt(r2*r3)
eq3 = x41 - x32 - 2sqrt(r3*r4)
eq4 = x1 - x31 - 2sqrt(r1*r3)
eq5 = (x1 - x42)^2 + (y42 - r1)^2 - (r1 + r4)^2
eq6 = r4/(a - x41) - r1/(a - x1)
eq7 = r3/(a - x32) - r1/(a - x1)
eq8 = r2/(a - x2) - r1/(a - x1)
eq9 = r3/x31 - r1/x1
eq10 = dist2(b, c, a, 0, x2, r2, r2)
eq11 = dist2(b, c, a, 0, x32, r3, r3)
eq12 = dist2(b, c, a, 0, x42, y42, r4)
eq13 = dist2(b, c, 0, 0, x31, r3, r3);
# eq01 = dist2(b, c, a, 0, x1, r1, r1)
# eq02 = dist2(b, c, 0, 0, x1, r1, r1)
# eq03 = dist2(b, c, 0, 0, x42, y42, r4)
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13], (a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42))

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)
   (a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42) = u
   return [
       -2*sqrt(r1)*sqrt(r2) - x1 + x2,  # eq1
       -2*sqrt(r2)*sqrt(r3) - x2 + x32,  # eq2
       -2*sqrt(r3)*sqrt(r4) - x32 + x41,  # eq3
       -2*sqrt(r1)*sqrt(r3) + x1 - x31,  # eq4
       (-r1 + y42)^2 - (r1 + r4)^2 + (x1 - x42)^2,  # eq5
       -r1/(a - x1) + r4/(a - x41),  # eq6
       -r1/(a - x1) + r3/(a - x32),  # eq7
       -r1/(a - x1) + r2/(a - x2),  # eq8
       -r1/x1 + r3/x31,  # eq9
       c*(a^2*c - 2*a^2*r2 + 2*a*b*r2 - 2*a*c*x2 + 2*a*r2*x2 - 2*b*r2*x2 - c*r2^2 + c*x2^2),  # eq10
       c*(a^2*c - 2*a^2*r3 + 2*a*b*r3 - 2*a*c*x32 + 2*a*r3*x32 - 2*b*r3*x32 - c*r3^2 + c*x32^2),  # eq11
       a^2*c^2 - 2*a^2*c*y42 - a^2*r4^2 + a^2*y42^2 + 2*a*b*c*y42 + 2*a*b*r4^2 - 2*a*b*y42^2 - 2*a*c^2*x42 + 2*a*c*x42*y42 - b^2*r4^2 + b^2*y42^2 - 2*b*c*x42*y42 - c^2*r4^2 + c^2*x42^2,  # eq12
       c*(-2*b*r3*x31 - c*r3^2 + c*x31^2),  # eq13
   ]
end;

r1 = 8/2
iniv = BigFloat[20.818, 5.28282, 9.60552, 7.7588, 2.27159, 12.81312, 1.27967, 2.1886, 16.22124, 0.74556, 18.20508, 5.53907, 8.57916]
res = nls(H, ini=iniv)

   ([20.807871178546222, 5.278637230682365, 9.621452184679763, 6.7569787601756, 2.2804168100957263, 12.79739836162148, 1.3000752069417922, 2.1961451399841465, 16.24106696141371, 0.7411783390746858, 18.204316902254963, 5.545974845788723, 8.58391116430616], true)

甲円の直径 = 8 のとき，乙円直径 = 4.56083，丙円直径 = 2.60015，丁円直径 = 1.48236 である。

function draw(r1, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42) = (20.818, 5.28282, 9.60552, 6.7588, 2.27159, 12.81312, 1.27967, 2.1886, 16.22124, 0.74556, 18.20508, 5.53907, 8.57916)
   (a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42) = res[1]
   @printf("甲円の直径 = %g，乙円直径 = %g，丙円直径 = %g，丁円直径 = %g\n", 2r1, 2r2, 2r3, 2r4)
   plot([0, a, b, 0], [0, 0, c, 0], color=:blue, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :green)
   circle(x31, r3, r3, :magenta)
   circle(x32, r3, r3, :magenta)
   circle(x41, r4, r4, :brown)
   circle(x42, y42, 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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
       point(x31, r3, " 丙円:r3,(x31,r2)", :magenta, :center, delta=-10delta, deltax=8delta)
       point(x32, r3, " 丙円:r3,(x32,r2)", :magenta, :left, :bottom, delta=12delta, deltax=-5delta)
       point(x41, r4, " 丁円:r4,(x41,r4)", :brown, :center, delta=-6delta, deltax=-5delta)
       point(x42, y42, " 丁円:r4,(x42,y42)", :brown, :left, :bottom, delta=2delta, deltax=5delta)
       point(b, c, "(b,c)", :blue, :center, :bottom, delta=delta)
       point(a, 0, "a", :blue, :left, :bottom, delta=delta)
       ylims!(-1, 10)
   end  
end;

draw(8/2, true)

算額（その1354）

算額（その1354）

十二 群馬県高崎市石原町 清水寺 安永3年(1774)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：三角形，面積最大化
#Julia, #SymPy, #算額, #和算

三辺が，大斜，中斜，小斜の三角形がある。中斜と中鈎（頂点から大斜へおろした垂線の長さ。三角形の高さ）の和が 126 寸，小斜と中鈎の和が 121 寸のとき，三角形の面積が最大になるときの中鈎はいかほどか。

大斜，中斜，小斜，中鈎，面積をそのまま変数名として，以下の連立方程式を解き，中鈎，大斜，小斜を求める。解は中斜の関数になる。

include("julia-source.txt");

using SymPy
@syms 小斜::positive, 中斜::positive, 大斜::positive,
     中鈎::positive, 面積::positive
eq1 = 中斜 + 中鈎 - 126
eq2 = 小斜 + 中鈎 - 121
eq3 = sqrt(中斜^2 - 中鈎^2) + sqrt(小斜^2 - 中鈎^2) - 大斜
res = solve([eq1, eq2, eq3], (中鈎, 大斜, 小斜))[1]

   (126 - 中斜, 11*sqrt(2*中斜 - 131) + 6*sqrt(7*中斜 - 441), 中斜 - 5)

中鈎と大斜は以下のようになるので，面積は 中鈎*大斜/2 である。

中鈎 = 126 - 中斜
大斜 = 11*sqrt(2*中斜 - 131) + 6*sqrt(7*中斜 - 441)

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

   126 - 中斜
   11*sqrt(2*中斜 - 131) + 6*sqrt(7*中斜 - 441)

面積 = (126 - 中斜)*(11*sqrt(2*中斜 - 131) + 6*sqrt(7*中斜 - 441))/2;

面積も中斜の関数になり，下図のように中斜が 80〜90 の間に最大値 3000 程を取る。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(面積, xlims=(60, 126), xlabel="中斜", ylabel="面積")

面積が最大になるときの中斜は，面積の導関数を取り，それが 0 になるときの値を求めればよい。

導関数 = diff(面積, 中斜);

ans_中斜 = solve(導関数, 中斜)[1].evalf()
ans_中斜 |> println

   84.8502093388565

中斜が 84.8502093388565 のときの面積は，2934.69516643686 である。

面積(中斜 => 84.8502093388565) |> println

   2934.69516643686

中鈎は「126 - 中斜」なので，126 - 84.8502093388565 = 41.1497906611435 寸である。

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   中斜 = 84.8502093388565
   (中鈎, 大斜, 小斜) = (126 - 中斜, 11*sqrt(2*中斜 - 131) + 6*sqrt(7*中斜 - 441), 中斜 - 5)
   大斜1 = sqrt(中斜^2 - 中鈎^2)
   plot([0, 大斜, 大斜1, 0], [0, 0, 中鈎, 0], color=:blue, lw=0.5)
   segment(大斜1, 0, 大斜1, 中鈎, :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)
       dimension_line(0, -3, 大斜, -3, "大斜", :blue, :center, delta=-4delta)
       point(大斜1, 中鈎/2, "中鈎", :red, :left, :vcenter, deltax=2delta, mark=false)
       point(大斜1/2, 中鈎/2, "中斜", :blue, :right, :vcenter, mark=false, deltax=-7delta)
       point((大斜1 + 大斜)/2, 中鈎/2, "小斜", :blue, :left, :vcenter, mark=false, deltax=7delta)
       ylims!(-15delta, 中鈎 + 2delta)
   end  
end;

draw(true)

算額（その1353）

算額（その1353）

中山道碓氷峠熊野社
藤田貞資門人 上州安中驛 角田左源司親信 享和元年(1801)
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

十一 群馬県碓氷郡松井田町峠 熊野神社 享和元年(1801)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

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

菱形の中に等円を数個描く（仮に11個描く）。菱形の対角線の長い方（菱長），短い方（菱平）の和と積，等円の直径が与えられたとき，等円の個数を求める術を述べよ。

菱形の中心を原点とする
菱長，菱平を 2a, 2b
菱長と菱平の和を K1, 菱長と菱平の積を K2
等円の半径を r
等円の中心座標を (x, y) として， x ≧ 0 かつ y = 0 の等円の個数を m，x = 0 かつ y ≧ 0 の等円の個数を n とする。

使う使わないはともかくとして，以下の 4 本の方程式を立てる。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r::positive,
     m, n,
     K1::positive, K2::positive
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2)

eq3 = 2a + 2b - K1
eq4 = (2a)*(2b) - K2;

1.　r, m, n を与えて a, b を求める場合

算額問題としては最も多いパターンである。 eq1, eq2 を連立させて a, b を求める。SymPy の能力上，r, m, n を未知数のままとして解くことができないので，r, m, n が既知として解く。
図の場合を例として， r = 1/2, m = 4, n = 3 として解く。

using SymPy
@syms a::positive, b::positive, r::positive,
     m, n,
     K1::positive, K2::positive
(r, m, n) = (1/2, 4, 3)
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2)
solve([eq1, eq2], (a, b))[1]

   (3.90138781886600, 2.60092521257733)

2 .* (3.90138781886600, 2.60092521257733)

   (7.802775637732, 5.20185042515466)

7.802775637732 + 5.20185042515466 |> println
7.802775637732 * 5.20185042515466 |> println

   13.004626062886661
   40.58887176852263

菱長 = 7.802775637732，菱平 = 5.20185042515466 である。
これを基礎として，図を書いたものが最初に示したものである。
後で使うことになる，菱長と菱平の和と積も 13.004626062886661, 40.58887176852263 である。

2.　r, a, b を与えられて m, n を求める場合

eq1, eq2 は独立なので，それぞれを m，n について解けばよい。

using SymPy
@syms a::positive, b::positive, r::positive,
     m, n,
     K1::positive, K2::positive
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2);

まずは m を求める。

ans_m = solve(eq1, m)[1]
ans_m |> println

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

前項のようにして別途計算しておいた結果，r = 1/2, a = 3.90138781886600, b = 2.60092521257733 を与えると，m = 4 が得られる。

ans_m(r => 1/2, a => 3.90138781886600, b => 2.60092521257733) |> println

   4.00000000000000

次に n を求め，r，a, b を代入し実際の値を求める。
n = 3 が得られる。

ans_n = solve(eq2, n)[1]
ans_n |> println
ans_n(r => 1/2, a => 3.90138781886600, b => 2.60092521257733) |> println

   b/(2*r) + 1 - sqrt(a^2 + b^2)/(2*a)
   3.00000000000000

a, b に近似値を与えると，m, n は整数から外れるが，四捨五入すればよい。

ans_m(r => 1/2, a => 3.9, b => 2.6) |> println
ans_n(r => 1/2, a => 3.9, b => 2.6) |> println

   3.99861218113400
   2.99907478742267

3.　r, K1, K2 を与えられて m, n を求める場合

eq1, eq2, eq3, eq4 を連立方程式として解けばよい。

using SymPy
@syms a::positive, b::positive, r::positive,
     m, n,
     K1::positive, K2::positive
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2)

eq3 = 2a + 2b - K1
eq4 = (2a)*(2b) - K2

res = solve([eq1, eq2, eq3, eq4], (m, n, a, b))[2]  # 2 of 2

   ((K1*r + K2/2 - r*sqrt(K1^2 - 4*K2) - r*sqrt(K1^2 - 2*K2))/(r*(K1 - sqrt(K1^2 - 4*K2))), (K1*r + K2/2 + r*sqrt(K1^2 - 4*K2) - r*sqrt(K1^2 - 2*K2))/(r*(K1 + sqrt(K1^2 - 4*K2))), K1/4 + sqrt(K1^2 - 4*K2)/4, K1/4 - sqrt(K1^2 - 4*K2)/4)

m, n はちょっと複雑な式になる。また，問では「長平の和と積」，この解答例では 長，平を 2a, 2b としている点に注意が必要である。
長平の和と積は 13.004626062886661, 40.58887176852263 である。

r, K1, K2 を代入して m，n は以下のようになる。

res[1](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println
res[2](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println

   4.00000000000000
   3.00000000000000

第 3，第 4 要素として a, b も求まる。
最初の項で求めた (3.90138781886600, 2.60092521257733) と一致する（当たり前だが）。

res[3](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println
res[4](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println

   3.90138781886600
   2.60092521257733

二数の和 u と積 v が与えられたとき，元の二数 x, y は以下の方程式を解けばわかる。

@syms u, v, x, y
eq01 = x + y ⩵ u
eq02 = x * y ⩵ v
res2 = solve((eq01, eq02), (x, y))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (u/2 - sqrt(u^2 - 4*v)/2, u/2 + sqrt(u^2 - 4*v)/2)
    (u/2 + sqrt(u^2 - 4*v)/2, u/2 - sqrt(u^2 - 4*v)/2)

長平の和と積を与えて，長，平を得る。それを 1/2 すれば，a, b である。
事前に a, b を求めて，前項に従えば同じ結果になる。

res2[2][1](u => 13.004626062886661, v => 40.58887176852263) |> println
res2[2][2](u => 13.004626062886661, v => 40.58887176852263) |> println

   7.80277563773200
   5.20185042515466

res2[2][1](u => 13.004626062886661, v => 40.58887176852263)/2 |> println
res2[2][2](u => 13.004626062886661, v => 40.58887176852263)/2 |> println

   3.90138781886600
   2.60092521257733

4.　閑話休題

問題で要求されているのは，等円の個数である。

求めた m, n から，2(m - 1) + 2(n - 1) + 1 = 2m + 2n - 3 で求めることができる。

m = 4, n = 3 のとき，個数は 11 個である。

include("julia-source.txt");

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, m, n) = (1/2, 4, 3)
   (a, b) = (3.90138781886600, 2.60092521257733)
   @printf("菱長 = %g，菱平 = %g\n", 2a, 2b)
   @printf("長平和 = %g，長平積 = %g\n", 2a + 2b, 2a*2b)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:blue, lw=0.5)
   circle(0, 0, r)
   for i = 1:m - 1
       circle2(i, 0, r)
   end
   for j = 1:n - 1
       circle22(0, j, r)
   end
   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", :blue, :left, :bottom, delta=delta)
       point(0, b, "b", :blue, :center, :bottom, delta=delta)
   end  
end;

draw(true)

 

算額（その1352）

算額（その1352）

九 群馬県群馬郡群馬町引間 妙見寺 寛政9年(1797)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：正方形，4等分割
#Julia, #SymPy, #算額, #和算

正方形の田んぼがある（45°回転）。図のように 3 本の割線で区切って，面積を 4 等分する。甲の長さはいかほどか。

正方形の一辺の長さを x
a, b, c を図のように定め，連立方程式を解く。
甲の長さは c - a + b である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive
c = x/√Sym(2)
eq1 = (c - b)^2 - x^2/4
eq2 = ((c - a) + (c - b - a)) * b/2 + (c - a)^2/2 - x^2/4
res = solve([eq1, eq2], (a, b))[1]  # 1 of 4

   (-x/2 + sqrt(2)*x - sqrt(-x^2 + (-x + 2*sqrt(2)*x)^2)/2, -x/2 + sqrt(2)*x/2)

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

甲の長さは正方形の一辺の長さの (x/sqrt(Sym(2)) - a + b) 倍である。
a, b を代入して，sqrt(2 - √2) 倍である。

甲 = (x/sqrt(Sym(2)) - a + b) 
甲 = 甲(a => res[1], b => res[2]) |> simplify
甲 |> println

   x*sqrt(2 - sqrt(2))

甲の長さは，正方形の一辺の長さの sqrt(2 - √2) 倍である。

正方形の一辺の長さが 10 寸のとき，甲は 10*sqrt(2 - √2) = 7.653668647301794 である。

x = 10
x*sqrt(2 - √2)

   7.653668647301794

function draw(x, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = x/√2
   (a, b) = (-x/2 + sqrt(2)*x - sqrt(-x^2 + (-x + 2*sqrt(2)*x)^2)/2, -x/2 + sqrt(2)*x/2)
   甲 = sqrt(2 - √2)*x
   @printf("正方形の一辺の長さが %g のとき，甲は %g\n", x, b + c - a)
   println((a, b))
   @printf("x = %g;  a = %g;  b = %g;  c = %g\n", x, a, b, c)
   plot([c, 0, -c, 0, c], [0, c, 0, -c, 0], color=:blue, lw=1)
   segment(c - b, -b, b - c, -b, :red)
   segment(a, -b, a, c - a, :red, lw=1)
   segment(-a, -b, -a, c - a, :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)
       point(c - b, -b, " (c-b,-b)")
       point(a, -b, "(a,-b)")
       point(a, c - a, " (a,c-a)", :green, :left, :vcenter)
       point(c, 0, " c", :green, :left, :vcenter)
       point(0, c, "c", :green, :center, :bottom, delta=delta/2)
       point(0, 0, " 0", :green, :left, :vcenter)
       point(a, (c - a - b)/2, " 甲", :green, :left, :vcenter, mark=false)
   end  
end;

draw(10, true)

 

算額（その1351）

算額（その1351）

九 群馬県群馬郡群馬町引間 妙見寺 寛政9年(1797)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：直角三角形2個
#Julia, #SymPy, #算額, #和算

2 つの直角三角形が共通の底辺を持つ。長股（高さの長い方）が 30 寸，短股（高さの短い方）が 20 寸のとき，矢（2 つの斜辺の交点から底辺へおろした垂線の長さ）はいかほどか。

注：底辺の長さにかかわらず矢は一定の長さである。

短股，長股，矢をそのまま変数名とする
斜辺の交点座標を (a, 矢) として，以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, 矢::positive,
     短股::positive, 長股::positive
eq = 矢/短股 - (長股 - 矢)/長股

res = solve(eq, 矢)[1]
res |> println

   短股*長股/(短股 + 長股)

矢は，短股*長股/(短股 + 長股) である。
長股，短股が 30 寸，20 寸のとき，矢は 12 寸である。

function draw(長股, 短股, a, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   矢 =  短股*長股/(短股 + 長股)
   b = a*長股/短股
   @printf("短股 = %g;  長股 = %g;  矢 = %g;  a = %g;  b = %g;  a + b = %g\n", 短股, 長股, 矢, a, b, a + b)
   plot([0, a + b, 0, 0, a + b, a + b, 0], [0, 0, 短股, 0, 長股, 0, 0], color=:blue, lw=1)
   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, 矢, "(a,矢)", :red, :center, delta=-3delta)
       point(0, 短股, "短股", :blue, :left, :bottom, deltax=2delta)
       point(a*(長股 + 短股)/短股, 長股, "(a*(長股+短股)/短股,長股)", :blue, :right, :vcenter, deltax=-4delta)
       # xlims!(-12delta, a + b + 10delta)
   end  
end;

draw(30, 20, 18, true)

 

 

算額（その1350）

算額（その1350）

二十九 群馬県高崎市山名町 八幡宮 文化11年(1814)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：直角三角形，正三角形，正方形，円
#Julia, #SymPy, #算額, #和算

直角三角形の中に，正三角形，正方形，円を容れる。直角三角形の面積を，直角三角形の三辺，正三角形の一辺，正方形の一辺，円の直径の 6 個の和で表す術を述べよ。

正三角形の一辺の長さを a
正方形の一辺の長さを c
直角三角形の 3 辺をそのまま「鈎」，「股」，「弦」
股 = 2a + b, b = a - c
円の半径を r
「直角三角形の三辺，正三角形の一辺，正方形の一辺，円の直径の 6 個の和」を K
とおく。

図中の（隠れているものも含む）直角三角形は，辺の比が 1:√3:2 であり，相似（合同）である（証明，説明省略）。よって，弦 = 2鈎，股 = √3鈎である。
また，正三角形の右にある三角形は頂角 120° の二等辺三角形である。
これらの長さの関係について以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎, 股, 弦, a, b, c, r, K
弦 = 2鈎
股 = √Sym(3)鈎
b = a - c
eq1 = 鈎 + 股 + 弦 + a + c + 2r - K
eq2 = 2a + b - 股
eq3 = (a + 2b)/√Sym(3) + √Sym(3)a - 弦
eq4 = (a - c + 2b) + c - (鈎 - 2c/√Sym(3) + √Sym(3)b) - 2r
eq5 = a - c + 2b - √Sym(3)c;

res = solve([eq1, eq2, eq4, eq5], (鈎, a, c, r));

鈎は長精度整数を含む分数式で表される。

res[鈎] |> println

   294669784542913660789*K/(1846975238059608604496 + 1066351650880288157387*sqrt(3)) + 170127679427900227272*sqrt(3)*K/(1846975238059608604496 + 1066351650880288157387*sqrt(3))

しかし，以下のようにあっけなく簡約化される。

ans_鈎 = res[鈎] |> simplify |> factor
ans_鈎 |> println

   K*(sqrt(3) + 8)/61

直角三角形の面積は 鈎*股/2 = √3*鈎^2 で，簡約化すると K^2*(48 + 67*sqrt(3))/7442 になる。

すなわち，「直角三角形の面積は，直角三角形の三辺，正三角形の一辺，正方形の一辺，円の直径の 6 個の和の平方の，(48 + 67√3)/7442 倍である。」

ans_鈎 * √Sym(3)ans_鈎 / 2 |> simplify |> println

   K^2*(48 + 67*sqrt(3))/7442

ちなみに，a, c, r も K の関数として，それぞれ以下のように表される。

res[a] |> simplify |> factor |> println
res[c] |> simplify |> factor |> println
res[r] |> simplify |> factor |> println

   K*(-5 + 7*sqrt(3))/61
   K*(-18 + 13*sqrt(3))/61
   -K*(-57 + 31*sqrt(3))/122

たとえば，K = 100 のとき，パラメータは以下のようになる。

K = 100;  鈎 = 15.9541816517523;  股 = 27.6334532140181;  弦 = 31.9083633035045;  a = 11.6792715622658;  b = 4.27491008948646;  c = 7.40436147277935;  r = 2.71018439784;

6 和が 100 になることが確かめられる。

鈎 + 股 + 弦 + a + c + 2r

   100.00000000000004

function draw(K, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   鈎 = K*(√3 + 8)/61
   a = K*(7√3 - 5)/61
   c = K*(13√3 - 18)/61
   r = K*(57 - 31√3)/122
   弦 = 2鈎
   股 = √3鈎
   b = a - c
   @printf("K = %.15g;  鈎 = %.15g;  股 = %.15g;  弦 = %.15g;  a = %.15g;  b = %.15g;  c = %.15g;  r = %.15g\n ", K, 鈎, 股, 弦, a, b, c, r)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r, r, r)
   plot!([b, b + a, b + a/2, b], [0, 0, √3a/2, 0], color=:blue, lw=0.5)
   plot!([0, √3c/2, √3c/2 + c/2, c/2, 0], [鈎 - 2c/√3, 鈎 - 2c/√3 - c/2, 鈎 - 2c/√3 - c/2 + √3c/2, 鈎 - 2c/√3 + √3c/2, 鈎 - 2c/√3], color=:green, 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(b, 0, "b", :blue, :center, delta=-delta)
       point(b + a, 0, "b+a", :blue, :center, delta=-delta)
       point(b + 2a, 0, "b+2a = 股", :blue, :center, delta=-delta)
       point(0, 鈎, "鈎 ", :black, :right, :vcenter)
       point(0, 鈎 - 2c/√3, "鈎-2c/√3 ", :black, :right, :vcenter)
       point(r, r, "(r,r)", :red, :center, delta=-delta)
       dimension_line(√3c/2 + c/2 + 1.1delta, 鈎 - 2c/√3 - c/2 + √3c/2 + 2delta, c/2 +  1.1delta, 鈎 - 2c/√3 + √3c/2 +  2delta, "c", :green, :center, :bottom, delta=2delta, deltax=2delta)
       dimension_line(b + 2delta, -1.1delta, b + (a - c)/2 + 2delta, √3(a - c)/2 - 1.1delta, "b", :blue, :left, deltax=3delta)
       plot!(xlims=(-28delta, 股 + 15delta), ylims=(-8delta, 鈎 + 5delta))
   end  
end;

draw(100, true)