算額（その2109）

算額（その2109）

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

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

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

include("julia-source.txt");

using SymPy

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

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

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

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

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

   6.31662479035540
   2.31662479035540

function draw(r2, r4, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (sqrt(r4)*sqrt(4*r2 + r4)/2 + 3*r4/2, sqrt(r4)*sqrt(4*r2 + r4)/2 - r4/2)
   r1 = (√r4*sqrt(4r2 + r4) + 3r4)/2
   r3 = (√r4*sqrt(4r2 + r4) - r4)/2
   @printf("小円の直径が %g，等円の直径が %g のとき，大円の直径は %g，甲円の直径は %g である。", 2r2, 2r4, 2r1, 2r3)
   plot()
   circle(0, 0, r1, :blue)
   circle2(r3 + r1, 0, r1, :blue)
   circle(0, r3 + r2, r2, :magenta)
   circle(0, 0, r3, :green)
   circle2(r3 + r4, 0, r4)
   circle(0, r3 + r4, r4)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(r3 + r1, 0, "大円:r1\n(r3+r1,0)", :blue, :center, delta=-delta/2)
       point(0, r3 + r2, "小円:r2\n(0,r3+r2)", :magenta, :center, :bottom, delta=delta/2)
       point(r3, 0, "r3 ", :green, :right, :vcenter)
       point(0, 0, "甲円", :green, :center, :vcenter, mark=false)
       point(r3 + r4, 0, "r3+r4", :red, :center, delta=-delta/2)
       point(r3 + r4, 0, "等円", :red, :center, :bottom, delta=delta/2, mark=false)
   end
end;

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

 

算額（その2108）

算額（その2108）

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

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

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

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

include("julia-source.txt");

using SymPy

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

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

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

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

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

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

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

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

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

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

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

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

   2.00000000000000

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

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

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

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

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

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

draw(5.4, 7.2, true)

#SymPy, #算額, #和算

算額（その2107）

算額（その2107）

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

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

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

include("julia-source.txt");

using SymPy

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

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

   2*r1 - r3

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

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

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

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

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

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

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

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

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

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

ans_x4 |> factor |> println

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

最後に r4 を求める。

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

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

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

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

   3.34783294256526

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

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

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

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

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

#SymPy #算額 #和算

算額（その2106）

算額（その2106）

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

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

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

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

include("julia-source.txt");

using SymPy

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

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

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

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

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

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

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

長弦を求める。

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

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

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

   1.44000000000000

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

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

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

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

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

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

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

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

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

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

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

   1.20000000000000

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

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

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

算額（その2105）

算額（その2105）

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

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

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

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

include("julia-source.txt");

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

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

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

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

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

2 * 0.5*(6 - 4√2)

   0.3431457505076194

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

draw(5, true)

算額（その2104）

算額（その2104）

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

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

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

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

   10.704662693192699

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

include("julia-source.txt");

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

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

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

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

draw(10/2, true)

算額（その2103）

算額（その2103）

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

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

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

include("julia-source.txt");

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

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

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

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

   15/2

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

draw(10/2, 8, true)

算額（その2102）

算額（その2102）

百八 群馬県邑楽郡板倉町板倉 雷電神社 慶応3年(1867)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，外円，二等辺三角形，弦3本，

外円の中に水平な弦と長さの同じ斜めの弦 2 本（倒立した二等辺三角形）と等円 2 個を容れる。外円の直径が 5 寸，斜線（二等辺三角形の斜辺）の長さが 3 寸のとき，等円の直径はいかほどか。

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

include("julia-source.txt");

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

   (-R + sqrt(2*R^2 - sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)) + sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)/(2*R), 斜*sqrt(4*R^2 - 斜^2)/(2*R))

等円の半径は -R + sqrt(2*R^2 - sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)) + sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4)/(2*R) である。
SymPy では sqrt(4*R^4 - 4*R^2*斜^2 + 斜^4) をうまく処理できないので，その部分を手作業で簡約化した後で simplify すると，斜*(2R - 斜)/(2R) になる。

問では「外円之径五寸□分」と 1 文字欠損のように見えるが，直径がちょうど 5 寸のとき（すなわち欠損文字は零か？），等円の直径が 2 寸 4 分になる。

2res[1](R => 5.0/2, 斜 => 3) |> println

   2.40000000000000

2(斜*(2R - 斜)/(2R))(R => 5.0/2, 斜 => 3) |> println

   2.40000000000000

function draw(R, 斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 斜*(2R - 斜)/(2R)
   x = 斜*sqrt(4R^2 - 斜^2)/(2R)
   @printf("外円の直径が %g，等斜の長さが %g のとき，等円の直径は %g である。\n", 2R, 斜, 2r)
   @printf("R = %g;  斜 = %g;  r = %g;  x = %g\n", R, 斜, r, x)
   plot([x, -x, 0, x], [-sqrt(R^2 - x^2), -sqrt(R^2 - x^2), -R, -sqrt(R^2 - x^2)], color=:green, lw=0.5)
   circle(0, 0, R, :blue)
   circle2(r, r - sqrt(R^2 - x^2), r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(r, r - sqrt(R^2 - x^2), "等円:r,(r,r-sqrt(R^2-x^2))", :red, :center, delta=-delta/2)
       point(x, -sqrt(R^2 - x^2), "(x,-sqrt(R^2-x^2))   ", :green, :right, delta=-delta/2)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
   end
end;

draw(5/2, 3, true)

算額（その2101）

算額（その2101）

百八 群馬県邑楽郡板倉町板倉 雷電神社 慶応3年(1867)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，直角三角形，正方形

直角三角形の中に正方形と中鈎（直角の頂点から斜辺への垂線）を入れる。区画された領域に，甲円と乙円を 1 つずつ容れる。甲円と乙円の直径がそれぞれ 8 寸，15 寸のとき，正方形の一辺の長さはいかほどか。

正方形の一辺の長さを a とする。
直角三角形の鈎と股（直角を挟む 2 辺のうちの短い方と長い方）を，「鈎」，「股」とする。
股は問題を解く上では必要ないが，図を描くために求める。

また，図に示すように，正方形と中鈎の延長線を補助線として，交点座標を (股 - a, a + x)
甲円と乙円を含む直角三角形は相似で，相似比が a:(鈎 - a) である。
甲円の半径と中心座標を r1, (股 - r1, a + r1)
乙円の半径と中心座標を r2, (股 - a + r2, r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, x::positive,
     鈎::positive, 股::positive,
     r1::positive, r2::positive;
eq1 = a/(鈎 - a) - (a + x)/a
eq2 = (鈎 - a) + a - sqrt((鈎 - a)^2 + a^2) - 2r1
eq3 = a + (a + x) - sqrt(a^2 + ( a + x)^2) - 2r2
eq4 = (鈎 - a)/a - 鈎/股;

SymPy の能力的に，一度には解けないので逐次的に解いていく。

まず，eq1, eq2 を連立させて，鈎，股を求める。

res = solve([eq1, eq4], (鈎, 股))[1]

   (a*(2*a + x)/(a + x), 2*a + x)

鈎，股を eq2 に代入し，x を求める。

eq12 = eq2(鈎 => res[1], 股 => res[2]) |> simplify |> numerator |> expand
ans_x = solve(eq12, x)[1]
ans_x |> println

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

鈎，股，x を eq3 に代入し，a を求める。

eq13 = eq3(鈎 => res[1], 股 => res[2], x => ans_x)
ans_a = solve(eq13, a)[3]  # 3 of 3
ans_a |> println

   r1 + r2 + sqrt(r1^2 + r2^2)

正方形の一辺の長さは，r1 + r2 + sqrt(r1^2 + r2^2) である。
甲円と乙円の直径がそれぞれ 8 寸，15 寸のとき，正方形の一辺の長さは 20 寸である。

ans_a(r1 => 8/2, r2 => 15/2).evalf() |> println

   20.0000000000000

a が求まれば x が決まり，更に鈎，股も決まる。
a = r1 + r2 + sqrt(r1^2 + r2^2)
x = a*(a^2 - 4a*r1 + 2r1^2)/(2r1*(a - r1))
鈎 = a*(2a + x)/(a + x)
股 = 2a + x

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

   r1 = 4;  r2 = 7.5;  a = 20;  x = 17.5;  鈎 = 30.6667;  股 = 57.5

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = r1 + r2 + sqrt(r1^2 + r2^2)
   x = a*(a^2 - 4a*r1 + 2r1^2)/(2r1*(a - r1))
   鈎 = a*(2a + x)/(a + x)
   股 = 2a + x
   @printf("甲円と乙円の直径がそれぞれ %g 寸，%g 寸のとき，正方形の一辺の長さは %g 寸である。\n", 2r1, 2r2,a)
   @printf("r1 = %g;  r2 = %g;  a = %g;  x = %g;  鈎 = %g;  股 = %g\n", r1, r2, a, x, 鈎, 股)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   rect(股 - a, 0, 股, a, :red)
   circle(股 - r1, a + r1, r1, :blue)
   circle(股 - a + r2, r2, r2, :magenta)
   segment(股 - a, a + x, 股, 0)
   segment(股 - a, a + x, 股 - a, a)
   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(股, 鈎, "(股,鈎) ", :green, :right, :bottom, delta=delta/2)
       point(股 - a, 0, "股-a ", :red, :right, :bottom, delta=delta)
       point(股 - a, a + x, "(股-a,a+x) ", :black, :right, delta=-delta/2)
       point(股 - r1, a + r1, "甲円:r1\n(股-r1,a+r1)", :black, :center, :bottom, delta=delta/2, deltax=4delta)
       point(股 - a + r2, r2, "乙円:r2\n(股-a+r2,r2)", :magenta, :center, delta=-delta/2)
       xlims!(-2delta, 股 + 5delta)
   end
end;

draw(8/2, 15/2, true)

算額（その2100）

算額（その2100）

百三 群馬県高崎市八幡町 八幡宮 安政7年(1860)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円3個，長方形

長方形の中に相似な楕円 3 個を容れる。長方形の短辺が与えられたとき，長辺はいかほどか。

長方形の長辺を c
楕円の長半径，短半径，中心座標を a, b, (b, a), (c - a, b)
とおき，以下の連立方程式を解く。
長方形の短辺は 2a = 4b である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     x0::positive, y0::positive;
eq1 = (x0 - b)^2/b^2 + (y0 - a)^2/a^2 - 1
eq2 = (x0 -c + a)^2/a^2 + (y0 - b)^2/b^2 - 1
eq3 = a^2*(x0 - b)/(b^2*(y0 - a)) - b^2*(x0 - c + a)/(a^2*(y0 - 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)
   (c, x0, y0) = u
   return [
-1 + (-b + x0)^2/b^2 + (-a + y0)^2/a^2,  # eq1
-1 + (-b + y0)^2/b^2 + (a - c + x0)^2/a^2,  # eq2
a^2*(-b + x0)/(b^2*(-a + y0)) - b^2*(a - c + x0)/(a^2*(-b + y0)),  # eq3
   ]
end;
a = 1/2
b = a/2
iniv = BigFloat[2.6475964903602565, 0.8542484563784245, 0.5046106449806166] .*a
res = nls(H, ini=iniv)

   ([1.4708869390890313, 0.4745824757657914, 0.28033924721145365], true)

長方形の長辺は，楕円の長半径の 2.9417738781780627 倍である。
楕円の長半径が 1/2（すなわち，長方形の短辺が 1）のとき，長辺は 1.47089 である。

楕円の長半径，短半径は 0.5, 0.25，接点の座標は (0.474582, 0.280339) である。

function draw(短辺, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 短辺/2
   b = a/2
   (c, x0, y0) = a .*[2.9417738781780627, 0.9491649515315828, 0.5606784944229073]
   @printf("長方形の短辺が %g のとき，長辺は %g である。\n楕円の長半径，短半径は %g, %g，接点の座標は (%g, %g) である。", 短辺, c, a, b, x0, y0)
   plot([0, c, c, 0, 0], [0, 0, 2a, 2a, 0], color=:green, lw=0.5)
   ellipse(b, a, b, a, color=:red)
   ellipse(c - a, b, a, b, color=:red)
   ellipse(c - a, 3b, a, b, 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)
       point(b, a, "(b,a)")
       point(c - a, b, "(c-a,b)")
       point(x0, y0, "(x0,y0)")
       point(c, 4b, "(c,2a)", :green, :right, :bottom, delta=delta/2)
   end

end;

draw(1, true)

算額（その2099）

#和算

算額（その2099）

八十七 群馬県碓氷郡松井田町峠 熊野神社 安政4年(1857)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：直角三角形，鈎股弦

直角三角形（鈎股弦）において，「鈎*股」が 52440 歩，「股*弦」が 130065 歩のとき，鈎，股，弦を求めよ。

以下の連立方程式を解く。

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     鈎股相乗::positive, 股弦相乗::positive;
eq1 = 鈎*股 - 鈎股相乗
eq2 = 股*弦 - 股弦相乗
eq3 = 鈎^2 + 股^2 - 弦^2
eq3 = sqrt(鈎^2 + 股^2) - 弦
res = solve([eq1, eq2, eq3], (鈎, 股, 弦))[2]  # 2 of 2

   (鈎股相乗*(1/(股弦相乗^2 - 鈎股相乗^2))^(1/4), (1/(股弦相乗^2 - 鈎股相乗^2))^(-1/4), 股弦相乗*(1/(股弦相乗^2 - 鈎股相乗^2))^(1/4))

鈎 = res[1](鈎股相乗 => 52440, 股弦相乗 => 130065)
股 = res[2](鈎股相乗 => 52440, 股弦相乗 => 130065)
弦 = res[3](鈎股相乗 => 52440, 股弦相乗 => 130065)
(鈎, 股, 弦)

   (152, 345, 377)

「鈎*股」が 52440 歩，「股*弦」が 130065 歩のとき，鈎 = 152 寸，股 = 345 寸，弦 = 377 寸である。

算額（その2098）

#和算

算額（その2098）

七十六 群馬県桐生市天神町 天満宮 嘉永5年(1852)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：年利

銀 76 匁を，複利で 7 年預け，元利合計が 357.62786865234375 匁であった。年利はいかほどか。

「術」は 400 文字以上にわたり述べられている。

「元本」，「年利」，「年数」，「元利合計」の関係は以下の式になる。

元利合計 = 元本*(1 + 年利)^年数

using SymPy
@syms 元本::positive, 年利::positive, 年数::positive,
     元利合計::positive;
eq = 元利合計 ⩵ 元本*(1 + 年利)^年数
eq |> println

   Eq(元利合計, 元本*(年利 + 1)^年数)

方程式を解いて，年利を求める。
年利 = (元利合計/元本)^(1/年数) - 1

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

   元利合計^(1/年数)/元本^(1/年数) - 1

「元本」，「年数」，「元利合計」を式に代入し，「年利」を得る。
元本 76 匁を，複利で 7 年預け，元利合計が 357.62786865234375 匁のとき，年利は 0.25 = 2 割 5 分 である。

res(元本 => 75, 年数 => 7, 元利合計 => 357.62786865234375).evalf() |>  println

   0.250000000000000

 

算額（その2097）

#和算

算額（その2097）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，直角三角形，正五角形

直角三角形の中に正五角形と甲円，乙円を容れる。甲円の直径が 5 寸のとき，乙円の直径はいかほどか。

正五角形を内包する円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1 - R*cosd(18), R)
乙円の半径と中心座標を r2(x2, r2 -R*cosd(36))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     x2::positive;
eq1 = dist2(0, R, R*cosd(Sym(18)), R*sind(Sym(18)), r1 - R*cosd(Sym(18)), R, r1)
eq2 = dist2(0, R, R*cosd(Sym(18)), R*sind(Sym(18)), x2, r2 - R*cosd(Sym(36)), r2)
eq3 = dist2(R*sind(Sym(36)), -R*cosd(Sym(36)), R*cosd(Sym(18)), R*sind(Sym(18)), x2, r2 - R*cosd(Sym(36)), r2);

SymPy の能力的に，一度には解けないので，逐次解いていく。

まず，eq1 を解き，R を求める。

ans_R = solve(eq1, R)[1] |> simplify
ans_R |> println

   r1*((-sqrt(10) + 5*sqrt(2))*sqrt(sqrt(5) + 5) + 8*sqrt(5))/10

eq2 の R に ans_R を代入し，x2 を求める。

eq12 = eq2(R => ans_R) |> simplify
ans_x2 = solve(eq12, x2)[1]  # 1 of 2
@syms d
ans_x2 = apart(ans_x2, d) |> simplify |> sympy.sqrtdenest|> simplify
ans_x2 |> println

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

eq3 の R に ans_R，x2 に ans_x2 を代入し，r2 を求める。

eq13 = eq3(R => ans_R, x2 => ans_x2) |> simplify |> numerator
ans_r2 = solve(eq13, r2)[2]  # 2 of 2
ans_r2 = apart(ans_r2/r1, d) |> x -> x*r1
ans_r2 |> println

   r1*(-sqrt(2)*sqrt(sqrt(5) + 5) - sqrt(5) + 3 + sqrt(10)*sqrt(sqrt(5) + 5)/2)

甲円の直径が 5 寸のとき，各パラメータは以下のように計算される。
乙円の直径は 2*3.032399997699489 = 6.064799995398978 である。

r1 = 5/2
t = sqrt(5 + √5)
u = sqrt(5 - √5)
r2 = r1*(3 - √5 + t*(√10/2 - √2))
x2 = (r2*(5√2 - √10) -　r1*(4√2 + 2*u))/((√5 - 3)*u)
R = r1*((5√2 - √10)*t + 8√5)/10
(r2, x2, R)

   (3.032399997699489, 8.347481064940814, 7.1007915155952475)

function draw(r1, more)
    pyplot(size=(600, 600), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(5 + √5)
   u = sqrt(5 - √5)
   r2 = r1*(3 - √5 + t*(√10/2 - √2))
   x2 = (r2*(5√2 - √10) -　r1*(4√2 + 2*u))/((√5 - 3)*u)
   R = r1*((5√2 - √10)*t + 8√5)/10
   @printf("甲円の直径が %g のとき，乙円の直径は %.15g である。\n", 2r1, 2r2)
   plot([R*cosd(18), 0, -R*cosd(18), -R*sind(36), R*sind(36), R*cosd(18)],
        [R*sind(18), R, R*sind(18), -R*cosd(36), -R*cosd(36), R*sind(18)], color=:magenta, lw=0.5)
   circle(r1 - R*cosd(18), R, r1)
   circle(x2, r2 - R*cosd(36), r2, :blue)
   (x01, y01) = Float64.(intersection(0, R, R*cosd(18), R*sind(18), -R*cosd(18), 0, -R*cosd(18), -R*sind(36)))
   (x02, y02) = Float64.(intersection(0, R, R*cosd(18), R*sind(18), R*sind(36), -R*cosd(36), -R*cosd(18), -R*cosd(36)))
   plot!([x01, x02, x01, x01], [y02, y02, y01, y02], 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(r1 - R*cosd(18), R, "甲円:r1,(r1-R*cosd(18),R)", :red, :left, :bottom, delta=delta, deltax=-8delta)
       point(x2, r2 - R*cosd(36), "乙円:r2,(x2,r2-R*cosd(36))", :blue, :left, :bottom, delta=delta, deltax=-8delta)
       point(0, R, "(0,R)", :magenta, :center, delta=-2delta)
       point(R*cosd(18), R*sind(18), "(R*cosd(18),Rsind(18)) ", :magenta, :right, :vcenter)
       point(R*sind(36), -R*cosd(36), "(R*sind(36),-R*cosd(36))", :magenta, :right, :bottom, delta=delta)
   end  
end;

draw(5/2, true)

算額（その2096）

算額（その2096）

六十四 群馬県安中市板鼻 鷹巣神社 天保11年(1840)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円10個，外円

外円の中に甲乙丙丁戊の円を容れる。甲円，乙円の直径が 4 寸，9 寸のとき，外円の直径はいかほどか。
異版では，名前の付け方が違うが，本問に対応付けると「丙円の直径も求めよ」としている。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (0, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5. y5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, y1::positive,
     r2::positive, y2::negative, r3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive,
     r5::positive, x5::positive, y5::positive;
eq1 = r1^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = (x4 - r1)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq3 = (x5 - r1)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq4 = r2^2 + (y3 - y2)^2 - (r2 + r3)^2
eq5 = (r2 - x4)^2 + (y4 - y2)^2 - (r2 + r4)^2
eq6 = (x5 - r2)^2 + (y5 - y2)^2 - (r2 + r5)^2
eq7 = x4^2 + (y3 - y4)^2 - (r3 + r4)^2
eq8 = (x5 - x4)^2 + (y5 - y4)^2 - (r4 + r5)^2
eq9 = r1^2 + y1^2 - (R - r1)^2
eq10 = r2^2 + y2^2 - (R - r2)^2
eq11 = x5^2 + y5^2 - (R - r5)^2;

using NLsolve

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

function H(u)
   (R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
       r1^2 - (r1 + r3)^2 + (y1 - y3)^2,  # eq1
       (-r1 + x4)^2 - (r1 + r4)^2 + (y1 - y4)^2,  # eq2
       (-r1 + x5)^2 - (r1 + r5)^2 + (y1 - y5)^2,  # eq3
       r2^2 - (r2 + r3)^2 + (-y2 + y3)^2,  # eq4
       -(r2 + r4)^2 + (r2 - x4)^2 + (-y2 + y4)^2,  # eq5
       (-r2 + x5)^2 - (r2 + r5)^2 + (-y2 + y5)^2,  # eq6
       x4^2 - (r3 + r4)^2 + (y3 - y4)^2,  # eq7
       -(r4 + r5)^2 + (-x4 + x5)^2 + (-y4 + y5)^2,  # eq8
       r1^2 + y1^2 - (R - r1)^2,  # eq9
       r2^2 + y2^2 - (R - r2)^2,  # eq10
       x5^2 + y5^2 - (R - r5)^2,  # eq11
   ]
end;

(r1, r2) = (4/2, 9/2)
iniv = BigFloat[9.1, 6.9, -1.1, 1.9, 3.5, 0.94, 2.8, 4.0, 1.8, 5.3, 5.1]

res = nls(H, ini=iniv)

   ([9.142857142857142, 6.857142857142857, -1.1428571428571428, 1.9393939393939394, 3.463203463203463, 0.9411764705882353, 2.823529411764706, 4.033613445378151, 1.7777777777777777, 5.333333333333333, 5.079365079365079], true)

甲円，乙円の直径が 4 寸，9 寸のとき，外円の直径は 128/7 = 18 + 2/7 ≒ 18.2857 である。
丙円は 128/33 = 3 + 29/33 ≒ 3.87879 である。

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

   R = 9.14286;  y1 = 6.85714;  y2 = -1.14286;  r3 = 1.93939;  y3 = 3.4632;  r4 = 0.941176;  x4 = 2.82353;  y4 = 4.03361;  r5 = 1.77778;  x5 = 5.33333;  y5 = 5.07937

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5) = [37.4/2, 12.9, -3.5, 7.9/2, 5.6, 2.3, 6.2, 5.9, 7.9/2, 12.2, 8.1]
   (R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5) = res[1]
   @printf("甲円，乙円の直径が %g, %g のとき，外円の直径は %g, 丙円の直径は %g である。\n", 2r1, 2r2, 2R, 2r3)
   @printf("R = %g;  y1 = %g;  y2 = %g;  r3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5)
   plot()
   circle(0, 0, R)
   circle2(r1, y1, r1, :blue)
   circle2(r2, y2, r2, :green)
   circle(0, y3, r3, :magenta)
   circle2(x4, y4, r4, :orange)
   circle2(x5, y5, r5, :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(0, R, "R", :red, :center, :bottom, delta=delta/2)
       point(r1, y1, "甲円:r1\n(r1,y1)", :blue, :center, delta=-delta/2)
       point(r2, y2, "乙円:r2,(r2,y2)", :green, :center, delta=-delta/2)
       point(0, y3, "丙円:r3\n(0,y3)", :magenta, :center, delta=-delta/2)
       point(x4, y4, "丁円:r4,(x4,y4)", :black, :left, :bottom, delta=delta/2, deltax=-2delta)
       point(x5, y5, "戊円:r5\n(x5,y5)", :brown, :center, :bottom, delta=delta/2)
   end  
end;

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

 

#和算

算額（その2095）

算額（その2095）

二十七 群馬県太田市細谷 冠稲荷神社 文化11年(1814)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個，長方形

長方形の中に，甲乙丙丁戊己の 6 個の円を容れる。甲円，丁円の直径がそれぞれ 169 寸，36 寸のとき，己円の直径はいかほどか。

長方形の長辺と短辺を a, b; b = 2r1
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, b - r3)
丁円の半径と中心座標を r4, (x4, r4)
戊円の半径と中心座標を r5, (x5, r5)
己円の半径と中心座標を r6, (r6, r6)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive,
     x2::positive, r3::positive, x3::positive, r4::positive, x4::positive,
     r5::positive, x5::positive, r6::positive;
b = 2r1
eq1 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (b - r3 - r2)^2 - (r2 + r3)^2
eq4 = (x4 - x3)^2 + (b - r3 - r4)^2 - (r4 + r3)^2
eq5 = (x5 - x3)^2 + (b - r3 - r5)^2 - (r5 + r3)^2
eq6 = (r6 - x3)^2 + (b - r3 - r6)^2 - (r6 + r3)^2
eq7 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq8 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
eq9 = (x5 - r6)^2 + (r5 - r6)^2 - (r5 + r6)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (a, r2, x2, r3, x3, x4, r5, x5, r6))

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, r2, x2, r3, x3, x4, r5, x5, r6) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (a - r1 - x2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - x3)^2,  # eq2
       -(r2 + r3)^2 + (x2 - x3)^2 + (2*r1 - r2 - r3)^2,  # eq3
       -(r3 + r4)^2 + (-x3 + x4)^2 + (2*r1 - r3 - r4)^2,  # eq4
       -(r3 + r5)^2 + (-x3 + x5)^2 + (2*r1 - r3 - r5)^2,  # eq5
       -(r3 + r6)^2 + (r6 - x3)^2 + (2*r1 - r3 - r6)^2,  # eq6
       (r2 - r4)^2 - (r2 + r4)^2 + (x2 - x4)^2,  # eq7
       (r4 - r5)^2 - (r4 + r5)^2 + (x4 - x5)^2,  # eq8
       (r5 - r6)^2 - (r5 + r6)^2 + (-r6 + x5)^2,  # eq9
   ]
end;

(r1, r4) = (169/2, 36/2)
iniv = BigFloat[350, 27, 171, 67, 115, 127, 20, 89, 36]

res = nls(H, ini=iniv)

   ([350.3565608160057, 26.68421052631579, 170.88675952162194, 66.89583333333333, 115.48770876656474, 127.05454353959865, 19.6047261009667, 89.48407269786442, 36.202314049586775], true)

甲円の直径が 169 寸，丁円の直径が 36 寸のとき，己円の直径は 72.4046 寸である。

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

   a = 350.357;  b = 169;  r1 = 84.5;  r2 = 26.6842;  x2 = 170.887;  r3 = 66.8958;  x3 = 115.488;  r4 = 18;  x4 = 127.055;  r5 = 19.6047;  x5 = 89.4841;  r6 = 36.2023

function draw(r1, r4, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r2, x2, r3, x3, x4, r5, x5, r6) = res[1]
   b = 2r1
   @printf("甲円の直径が %g，丁円の直径が %g のとき，己円の直径は %g である。\n", 2r1, 2r4, 2r6)
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 = %g;  r5 = %g;  x5 = %g;  r6 = %g\n", a, b, r1, r2, x2, r3, x3, r4, x4, r5, x5, r6)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
   circle(a - r1, r1, r1)
   circle(x2, r2, r2, :magenta)
   circle(x3, b - r3, r3, :blue)
   circle(x4, r4, r4, :orange)
   circle(x5, r5, r5, :purple)
   circle(r6, r6, r6, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a, b, "(a,b)", :green, :right, :bottom, delta=2delta)
       point(a - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-2delta)
       point(x2, r2, "乙円:r2,(x2,r2)", :magenta, :left, :bottom, delta=2delta)
       point(x3, b - r3, "丙円:r3,(x3,b-r3)", :blue, :center, delta=-2delta)
       point(x4, r4, " 丁円:r4,(x4,r4)", :orange, :left, :vcenter)
       point(x5, r5, "戊円:r5,(x5,r5) ", :purple, :right, :vcenter)
       point(r6, r6, "己円:r6\n(r6,r6) ", :brown, :center, :bottom, delta=2delta)
   end  
end;

draw(169/2, 36/2, true)

#和算