算額（その1313）
百三十八 群馬県利根郡月夜野町上津 八幡神社 明治22年(1889)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円39個，円環
#Julia, #SymPy, #算額, #和算
大円の内外にそれぞれ 39 個の等円（外等円，内等円と呼ぶ）が互いに接している。
外等円と内等円の直径の差が 8.454 寸，大円と内等円の直径の和が 171.305 寸，「外余積が 109.688535160125 歩」のとき，大円，外等円，内等円の直径はいかほどか。
算額（その672）と同じであるが，問題の提示法に難がある。
大円の半径と中心座標を R, (0, 0)
外等円の半径と中心座標を r1, ((R + r1)*cos(pi/19), (R + r1)*sin(pi/19)); (R + r1)*sin(pi/19) = r1
内等円の半径と中心座標を r2, ((R - r2)*cos(pi/19), (R - r2)*sin(pi/19)); (R - r2)*sin(pi/19) = r2
とする。
未知数は R, r1, r2 の 3 個なので，方程式も 3 本用意すれば解ける。
外等円と内等円は互いに接し合っているという 2 条件は eq1, eq2 で表されている。
残り1つの条件が必要であるが，「問」には 3 個の条件が提示されている。3 番目の条件に出てくる「外余積」は有効数字が多いので，この条件を使えば正確な解が得られると思われるが，この「外余積」が何であるかわからない。
残る 2 条件のいずれを使うかで解が異なる。それは「外等円と内等円の直径の差が 8.454 寸」も「大円と内等円の直径の和が 171.305 寸」も有理数ではなく誤差を含むのが原因である。
更に，後で分かるが，eq1, eq2 に使われている「円周率」は 3.16 である（円積率=0.79）。
1. 「外等円と内等円の直径の差が 8.454 寸」を使う場合
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive, r2::positive,
円周率::positive;
円周率 =3.16
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
eq3 = 2(r1 - r2) - 8.454;
res = solve([eq1, eq2, eq3], (r1, r2, R))
res |> println
Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(r2 => 10.6530275137650, R => 75.0022912276876, r1 => 14.8800275137650)
大円の直径 = 2*75.0022912276876 = 150.0045824553752
外等円の直径 = 2*14.8800275137650 = 29.76005502753
内等円の直径 = 2*10.6530275137650 = 21.30605502753
2. 「大円と内等円の直径の和が 171.305 寸」を使う場合
@syms R::positive, r1::positive, r2::positive,
円周率::positive;
円周率 =3.16
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
eq3 = 2(R + r2) - 171.305;
res = solve([eq1, eq2, eq3], (r1, r2, R))
res |> println
Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(r2 => 10.6526769444053, R => 74.9998230555947, r1 => 14.8795378424756)
大円の直径 = 2*74.9998230555947 = 149.9996461111894
外等円の直径 = 2*14.8795378424756 = 29.7590756849512
内等円の直径 = 2*10.6526769444053 = 21.3053538888106
3. 大円の直径を与える場合
r1, r2 は 大円 R の関数なので，条件は eq1, eq2 だけで十分である。求められた r1, r2 は R を含む式なので，R として任意の値を設定すれば，その場合の外等円，内等円の半径（直径）はすぐに得られる。
大円の直径 2R は，前 2 項から 150 寸と推定される（実際に「答」で 150 寸とされている）。
@syms R::positive, r1::positive, r2::positive,
円周率::positive;
円周率 = Sym(316)/100
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
ans_r1 = solve(eq1, r1)[1]
ans_r2 = solve(eq2, r2)[1]
(ans_r1, ans_r2) |> println
(2ans_r1(R => 150/2).evalf(), 2ans_r2(R => 150/2).evalf())
(R*sin(79/475)/(1 - sin(79/475)), R*sin(79/475)/(sin(79/475) + 1))
(29.7591458944761, 21.3054041537713)
外等円，内等円の半径はそれぞれ
R*sin(79/475)/(1 - sin(79/475))
R*sin(79/475)/(sin(79/475) + 1)
であり，R = 150/2 を代入し 2 倍すると，それぞれの直径は 29.7591458944761, 21.3054041537713 となり，小数点以下 3 桁までとると，29.759, 21.305 となり，「答」と一致する。
4. eq1, eq2 を，円周率 = π で計算する場合
近似計算ではなく正確な解を求める。
@syms R::positive, r1::positive, r2::positive,
円周率::positive;
円周率 = PI
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
ans_r1 = solve(eq1, r1)[1]
ans_r2 = solve(eq2, r2)[1]
(ans_r1, ans_r2) |> println
(2ans_r1(R => 150/2).evalf(), 2ans_r2(R => 150/2).evalf())
(R*sin(pi/19)/(1 - sin(pi/19)), R*sin(pi/19)/(sin(pi/19) + 1))
(29.5535416156890, 21.1998138649765)
大円の直径が 150 寸のとき，外等円の直径は 29.5535 寸, 内等円の直径は 21.1998 寸である。
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2) = (R*sin(pi/19)/(1 - sin(pi/19)), R*sin(pi/19)/(sin(pi/19) + 1))
@printf("大円の直径が %g のとき，外等円の直径は %g, 内等円の直径は %g である。\n", 2R, 2r1, 2r2)
@printf("R = %g; r1 = %g; r2 = %g\n", R, r1, r2)
plot()
circle(0, 0, R, :magenta)
circle(0, 0, R + r1, :gray90)
circle(0, 0, R - r2, :gray90)
rotate(0, -r1 - R, r1, :blue, angle=360/19)
rotate(0, r2 - R, r2, angle=360/19)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
end
end;
draw(150/2, true)
